スペクトル
フーリエスペクトルの計算 †
関数fftでは,データ長が値にかかってくる.
f: データ fの振幅とスペクトルの振幅を合わせたいときは P = 2*abs(fft(f))/length(f) Powersectrum P = c.*conj(c)/length(c);
とする.
Sample of cpsd †
Fs = 1/dt; filt = hanning(nfft/16); n_overlap = 15; [P_cpsd, f] = cpsd(x,x,filt, n_overlap, nfft, Fs);
Sample of csd †
Fs = 1/dt; filt = hanning(nfft/16); n_overlap = 15; [S,freq] = csd( eta2, eta2, nfft, Fs, filt, n_overlap, 'mean' );
パワースペクトルの計算方法 †
パワースペクトルは,複素フーリエ係数cにより定義するため,fftの出力から c*conj(c)で計算できる. しかし,パワースペクトルは確率密度関数であり,適切なフィルターを用いてスムージングを掛ける必要がある. MATLABでは,フィルター付パワースペクトル推定法としてcsd, cpsdが用意されている.
データの用意 †
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で直接計算 †
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で計算 †
[P_cpsd, f] = cpsd(x,x,[],[],nfft,Fs); disp('mean energy of x by cpsd') sum(P_cpsd)/length(x)*Fs
csdで計算 †
[P_csd, f] = csd(x,x,nfft,Fs); disp('mean energy of x by csd') sum(P_csd)/length(f)
ETC †
cpsdとcsd †
csdはそのうち無くなり,cpsdが今後サポートされる予定. しかし,cpsdはcpdとちょっと仕様が異なるので移行に注意が必要. cpsdはパワースペクトルが周波数で割られており,単位周波数辺りのパワースペクトル密度となっている(こっちが一般的). csdはそうでなく,fftベースで求められたパワー量をベースとしている.
信号のパワーの平均を得るためにはcsdの方は,fftベースと同様に周波数点数で割るが,cpsdの方は更に周波数Fsを掛けることで求めることができる.