2016-11-08 66 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