- 追加された行はこの色です。
- 削除された行はこの色です。
*フーリエスペクトルの計算
*フーリエスペクトルの計算 [#bc0c7995]
関数fftでは，データ長が値にかかってくる．
f: データ
fの振幅とスペクトルの振幅を合わせたいときは
P = 2*abs(fft(f))/length(f)
Powersectrum
P = c.*conj(c)/length(c);
とする．
Original routine is
*** Sample of cpsd [#eb73fa00]
Fs = 1/dt;
filt = hanning(nfft/16);
n_overlap = 15;
[P_cpsd, f] = cpsd(x,x,filt, n_overlap, nfft, Fs);
*** Sample of csd [#eb73fa00]
Fs = 1/dt;
filt = hanning(nfft/16);
n_overlap = 15;
[S,freq] = csd( eta2, eta2, nfft, Fs, filt, n_overlap, 'mean' );
* パワースペクトルの計算方法 [#l26be0c9]
パワースペクトルは，複素フーリエ係数cにより定義するため，fftの出力から
c*conj(c)で計算できる．
しかし，パワースペクトルは確率密度関数であり，適切なフィルターを用いてスムージングを掛ける必要がある．
MATLABでは，フィルター付パワースペクトル推定法としてcsd, cpsdが用意されている．
*** データの用意 [#w563a5ef]
dt = 0.001;
t = 0:dt:0.6;
tmax = max(t);
df = 1/tmax;
Fs = 1/dt;
x = 5*sin(2*pi*50*t)+sin(2*pi*120*t);
nfft = length(x);
plot(t,x)
disp('mean energy of x')
var(x)
*** fftで直接計算 [#naf72ac4]
c = fft(x,nfft);
S_fft = c.*conj(c)/length(x);
P_fft = S_fft;
disp('mean energy of x by fft')
sum(P_fft)/length(x)
*** cpsdで計算 [#k9fe556e]
[P_cpsd, f] = cpsd(x,x,[],[],nfft,Fs);
disp('mean energy of x by cpsd')
sum(P_cpsd)/length(x)*Fs
*** csdで計算 [#p6586c74]
[P_csd, f] = csd(x,x,nfft,Fs);
disp('mean energy of x by csd')
sum(P_csd)/length(f)
* ETC [#nf9204d4]
** cpsdとcsd [#p6a18c8e]
csdはそのうち無くなり，cpsdが今後サポートされる予定．
しかし，cpsdはcpdとちょっと仕様が異なるので移行に注意が必要．
cpsdはパワースペクトルが周波数で割られており，単位周波数辺りのパワースペクトル密度となっている（こっちが一般的）．
csdはそうでなく，fftベースで求められたパワー量をベースとしている．
信号のパワーの平均を得るためにはcsdの方は，fftベースと同様に周波数点数で割るが，cpsdの方は更に周波数Fsを掛けることで求めることができる．
