| Signal Processing Toolbox | |
表示
Pxy=csd(x,y) Pxy=csd(x,y,nfft) [Pxy,f]=csd(x,y,nfft,fs) Pxy=csd(x,y,nfft,fs,window) Pxy=csd(x,y,nfft,fs,window,numoverlap) Pxy=csd(x,y,...,'dflag') [Pxy,Pxyc,f]=csd(x,y,nfft,fs,window,numoverlap,p) csd(x,y,...)
詳細
Pxy は、スペクトル推定にWelch法を用いて、長さ = csd(x,y)
nのデータ列xとyのクロススペクトル密度を推定します。Pxy = csd(x,y)は、つぎのデフォルト値を使用します。
nfft = min(256,length(x))fs = 2window は、長さnfftの周期性をもつHann(Hanning)ウィンドウnumoverlap = 0nfftは、csdが使用するFFTの長さを設定します。この値により、クロススペクトルを推定する周波数が決定されます。fsは、サンプリング周波数を設定するスカラ値です。windowは、csdが分割したxおよびyベクトルに使用するウィンドウ関数とウインドウ長を設定します。numoverlapは、分割したデータ列をオーバラップさせるウインドウ長です。入力引数からどの引数を省略しても、上記のデフォルト値が使用されます。
xとyが実数の場合、csdは正の周波数のクロススペクトル密度のみを推定します。この場合、出力Pxyは、nfftが偶数の場合、長さ(nfft/2 + 1)の列ベクトル、nが奇数の場合、長さ(nfft + 1)/2の列ベクトルとなります。xまたはyが複素数の場合、csdは正と負の周波数のクロススペクトル密度を共に推定し、Pxyの長さはnfftとなります。
Pxy は、 = csd(x,y,nfft)
xとyのクロススペクトル密度の推定に、長さnfftのFFTを使用します。実行速度を最も速くするには、nfftを2のベキ乗として設定します。
[Pxy,f] は、CSDを評価する周波数のベクトル = csd(x,y,nfft,fs)
fを出力します。fは、Pxyと同じ大きさであるため、plot(f,Pxy)は適切に換算された周波数のスペクトルをプロットします。fsは周波数を求めるだけのもので、出力Pxyには影響しません。
Pxy は、ベクトル = csd(x,y,nfft,fs,window)
xおよびyの分割したものを対象にウィンドウ関数およびウインドウ長を設定します。windowにスカラを与えた場合、csdは、その値を長さとするHannウィンドウを使用します。ウィンドウの長さは、nfft以下でなければなりません。ウィンドウの長さがnfftより小さい場合、csdは、ゼロを付加します。ウィンドウの長さがnfftを越えている場合、csdはエラーを出力します。
Pxy は、 = csd(x,y,nfft,fs,window,numoverlap)
numoverlapサンプルでxおよびyの区分を重複させます。
入力引数xまたはy以外のどの入力引数に対しても、空行列を使って、デフォルト値を設定することができます。たとえば、
csd(x,y,[],10000)
は、つぎのものと等価ですが、デフォルトの2 Hzの代わりに10,000 Hzのサンプリング周波数を使用します。
csd(x)
Pxy では、トレンド除去オプションを設定します。ここで、 = csd(x,y,...,'dflag')
'dflag'には、つぎの文字列を設定できます。
'linear'：ウィンドウを適用する前に、x、yから、最適近似直線を除去します。'mean'：ウィンドウを適用する前に、x、yから各々の平均値を除去します。'none'：トレンド除去を行ないません(デフォルト) 引数dflagは、入力引数のリストの最後に付けなければなりません。csdは、省略した入力引数の数に関係なく、dflag文字列を認識します。
[Pxy,Pxyc,f] は、 = csd(x,y,nfft,fs,window,numoverlap,p)
Pxyの
p*100パーセントの信頼区間の推定値を含んだベクトルPxycを出力します。ここで、pは0と1の間の正のスカラです。Pxycは、Pxyと同じ長さの2列の行列となります。区間[Pxyc(:,1), Pxyc(:,2)]は、確率pで真のCSDを示します。plot(f,[Pxy Pxyc])は、p*100パーセントの信頼区間内のクロススペクトルをプロットします。未設定の場合、pは0.95にデフォルト設定されます。
csd(x,y,...)
は、周波数に対するCSDを、カレントのFigureウィンドウにプロットします。p パラメータを設定した場合、プロットには信頼区間も含まれます。
例題
2つの有色ノイズ信号を生成し、95%の信頼区間と共に、そのCSDをプロットします。長さ1024のFFT、500点の非重複三角ウィンドウ、および10 Hzのサンプリング周波数を設定します。
randn('state',0);
h = fir1(30,0.2,boxcar(31));
h1 = ones(1,10)/sqrt(10);
r = randn(16384,1);
x = filter(h1,1,r);
y = filter(h,1,x);
csd(x,y,1024,10000,triang(500),0,[])
アルゴリズム
csdは、つぎのようにスペクトル密度推定のWelch法を使っています(参考文献[1]および[2]を参照)。
windowwベクトルで設定したウィンドウを適用します。
nfft点のFFTを使って変換します。
xの変換したものとyの変換したものの共役との積をスケーリングすることで、各分割のピリオドグラムを作成します。
xとyのクロススペクトル密度Pxyを作成します。
csdが平均化に使う分割数はkです。ここで、kは、つぎのように表されます。
fix((length(x)-numoverlap)/(length(window)-numoverlap))
診断
csdに対して、誤った引数を使用した場合、つぎのような警告メッセージが表示されます。
Requires window's length to be no greater than the FFT length.Requires NOVERLAP to be strictly less than the window length.Requires positive integer values for NFFT and NOVERLAP.Requires vector (either row or column) input.Requires inputs X and Y to have the same length.Requires confidence parameter to be a scalar between 0 and 1.
参考
|
2つの信号間の2乗コヒーレンス関数の推定 |
|
Burg法を使ったパワースペクトルの推定 |
|
Multitaper法(MTM)を使ったパワースペクトル推定 |
|
MUSIC アルゴリズムを使った擬似スペクトルの推定 |
|
Welch法を使って信号のパワースペクトルの推定 |
|
Yule-Walker AR法を使って、パワースペクトルの推定 |
|
入力および出力からの伝達関数の推定 |
参考文献
[1] Rabiner, L.R., and B. Gold. Theory and Application of Digital Signal Processing. Englewood Cliffs, NJ: Prentice-Hall, 1975. Pgs. 414-419.
[2] Welch, P.D. "The Use of Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging Over Short, Modified Periodograms." IEEE Trans. Audio Electroacoust. Vol. AU-15 (June 1967). Pgs. 70-73.
[3] Oppenheim, A.V., and R.W. Schafer. Discrete-Time Signal Processing. Englewood Cliffs, NJ: Prentice-Hall, 1989. Pg. 737.
| | cremez | czt | |