スペクトル
*フーリエスペクトルの計算 [#bc0c7995]
関数fftでは,データ長が値にかかってくる.
f: データ
fの振幅とスペクトルの振幅を合わせたいときは
P = 2*abs(fft(f))/length(f)
Powersectrum
P = c.*conj(c)/length(c);
とする.
*** 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が今後サポートされる予定.
csdはそのうち無くなり,cpsdが今後サポートされる予定.
しかし,cpsdはcpdとちょっと仕様が異なるので移行に注意が必要.
cpsdはパワースペクトルが周波数で割られており,単位周波数辺りのパワースペクトル密度となっている(こっちが一般的).
csdはそうでなく,fftベースで求められたパワー量をベースとしている.
信号のパワーの平均を得るためにはcsdの方は,fftベースと同様に周波数点数で割るが,cpsdの方は更に周波数Fsを掛けることで求めることができる.