2つのアプローチを使用して信号のパワースペクトルを見つける

Question

私は信号のパワースペクトルを見つけようとしています。信号の長さは100000、サンプル周波数は1000Hz、ポイント数は100000です。私は2つのアプローチを使ってパワースペクトルを見つけました。最初の1つは、すべての長さを1つの部分として取り、そのパワースペクトルを見つけ、2番目の方法は、信号を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); 

コードを更新するだけです。私は分かりませんが、ノイズを加えて2つのプロットがシフトしたように見えます。

Answer 1

主な問題はreshapeであり、各行が別々のシーケンスで作業しています。ただし、最初の列を塗りつぶしてから2番目の列に移動します。

代わりに以下を使用できます。

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

正規化も別の問題です。デフォルトで正規化されているので、fftの代わりにifftを使用することができます（理由はわかりません）。または、代わりに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)が同じになるようにベクトルの長さを得るためのハック方法です以下を使用プロットする終わり

。