Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

クロス スペクトルと振幅二乗コヒーレンス

この例は、クロス スペクトルを使用して、二変量時系列における正弦波成分の間の位相遅れを求める方法を示しています。この例では、また、正弦波の周波数における著しい周波数領域の相関を特定するために振幅二乗コヒーレンスを使用します。

二変量時系列を作成します。個々の系列は、周波数が 100 Hz および 200 Hz の 2 つの正弦波から構成されます。系列は加法性ホワイト ガウス ノイズに組み込まれ、1 kHz でサンプリングされます。系列 x の正弦波はともに振幅 1 です。系列 y の 100 Hz の正弦波は振幅 0.5 で、系列 y の 200 Hz の正弦波は振幅 0.35 です。系列 y の正弦波には、それぞれ π/4 ラジアン (100 Hz) および π/2 ラジアン (200 Hz) の位相遅れがあります。系列 y は、ノイズで破損された、入力 x を持つ線形システムの出力と見なすことができます。再現性のある結果を得るために、乱数発生器を既定の状態に設定します。

rng default

Fs = 1000;
t = 0:1/Fs:1-1/Fs;

x = cos(2*pi*100*t) + sin(2*pi*200*t) + 0.5*randn(size(t));
y = 0.5*cos(2*pi*100*t - pi/4) + 0.35*sin(2*pi*200*t - pi/2) + 0.5*randn(size(t));

二変量時系列について振幅二乗コヒーレンス推定を求めます。振幅二乗コヒーレンスによって、2 つの時系列間の著しい周波数領域の相関を特定することができます。クロス スペクトルの位相推定は、著しい周波数領域相関がある場合にのみ有用です。

求められた振幅二乗コヒーレンス推定がすべての周波数に対して等しく 1 となることを避けるため、平均化されたコヒーレンス推定子を使用しなければなりません。ウェルチのオーバーラップ セグメント平均 (WOSA) とマルチテーパー手法のどちらも適していますが、mscohere は WOSA 推定子を実装します。

ウィンドウの長さを 100 サンプルに設定します。このウィンドウの長さに 100 Hz の正弦波 10 周期および 200 Hz の正弦波 20 周期が含まれます。既定のハミング ウィンドウで 80 サンプルのオーバーラップを使用します。出力周波数が Hz 単位となるように明示的にサンプル レートを入力します。振幅二乗コヒーレンスをプロットします。振幅二乗コヒーレンスは 100 Hz および 200 Hz において 0.8 より大きくなります。

[Cxy,F] = mscohere(x,y,hamming(100),80,100,Fs);

plot(F,Cxy)
title('Magnitude-Squared Coherence')
xlabel('Frequency (Hz)')
grid

Figure contains an axes object. The axes object with title Magnitude-Squared Coherence, xlabel Frequency (Hz) contains an object of type line.

xy のクロス スペクトルを cpsd を使用して求めます。コヒーレンス推定で使用したのと同じパラメーターを使用してクロス スペクトルを求めます。コヒーレンスが小さい場合、クロス スペクトルを無視します。クロス スペクトルの位相をプロットして、2 つの時系列間で著しいコヒーレンスがある周波数を示します。正弦波成分の間の既知の位相遅れをマークします。100 Hz および 200 Hz において、クロス スペクトルから推定された位相遅れが真の値に近くなります。

[Pxy,F] = cpsd(x,y,hamming(100),80,100,Fs);

Pxy(Cxy < 0.2) = 0;

plot(F,angle(Pxy)/pi)
title('Cross Spectrum Phase')
xlabel('Frequency (Hz)')
ylabel('Lag (\times\pi rad)')
grid

Figure contains an axes object. The axes object with title Cross Spectrum Phase, xlabel Frequency (Hz), ylabel Lag ( times pi blank rad) contains an object of type line.

参考

| |