2016-11-08 19 views
0

我正试图找到信号的功率谱。信号长度为100000,采样频率为1000Hz,点数为100000。我使用两种方法找到功率谱。第一种方法是将所有长度作为一个部分并找到功率谱,而第二种方法是将信号分成100*1000并找出每一行的频谱,然后得到所有行的平均值。我的问题是,我必须在两种方法中得到相同的答案,但我得到了不同的答案。我不知道我的代码中有什么错误。使用两种方法找到信号的功率谱

N=100000; 
SF=1000;  
a=0.1; 
b=0.3; 
amplitude1=1; 
amplitude2=0.5; 
t=0:1/SF:100; 
f1=SF*a; 
f2=SF*b; 
A=amplitude1*sin(2*pi*f1*t)+amplitude2*sin(2*pi*f2*t); 
Y=2*randn(1,length(A))+A; 
bin=[0 :N/2]; 
fax_Hz=(bin*SF)/N; 
FFT=fft(Y); 
spectra=2/(SF*length(Y))*(FFT.*conj(FFT)); 
plot(fax_Hz,spectra(1,1:50001)); 
D=reshape(Y(1,1:100000),[100,1000]); 
M=length(D(1,:)); 
for i=1:100 
    FFT_1(i,:)=fft(D(i,:)); 
    S(i,:)=(2/(SF*M))*(FFT_1(i,:).*conj(FFT_1(i,:))); 
end 
S_f=mean(S); 
figure 
plot (S_f); 

我只是更新了代码。我不知道,但是当我添加噪音来表明这两个地块看起来转移。

回答

1

主要问题是与reshape你正在使用每一行是一个单独的序列。然而,在转移到第二个专栏之前,Reshape会填满第一个专栏。

您可以改用以下代码。

D=reshape(A(1,1:100000),[1000,100]).'; 

标准化是另一个问题。您可以使用ifft而不是fft,因为它默认为标准化(不知道为什么)。或者,也可以保持规范化,而不是使用mean,则应该使用sum,这可能是由于您可能犯的错误。在幅度上似乎仍然存在一个小的差异,但不确定这些差异是从哪里来的。

在最后绘制使用以下内容:

bin=[0 :N]; 
fax_Hz=(bin*SF)/N; 
FFT=ifft(A); 
spectra=FFT.*conj(FFT); 
plot(fax_Hz,spectra); hold on 
D=reshape(A(1,1:100000),[1000,100]).'; 
M=length(D(1,:)); 
for i=1:100 
    FFT_1(i,:)=ifft(D(i,:)); 
    S(i,:)=FFT_1(i,:).*conj(FFT_1(i,:)); 
end 
S_f=mean(S); 
plot(fax_Hz(1:100:end-1), S_f); 

注:fax_Hz(1:100:end-1)是获得向量的长度是相同的哈克的方式。

+0

谢谢你的帮助。请,我不知道,但是当我给信号添加噪音时,这两个地块看起来转移了。我只是更新代码来看看。 – user6052232

+0

我无法重现该问题。从你的代码中删除下面'bin'中的所有内容,并将其替换为我发布的代码。 (显然将'A'重命名为'Y') – mpaskov