Navigation : Top / MATLAB / スペクトル

スペクトル

フーリエスペクトルの計算

関数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を掛けることで求めることができる.