mscohere

振幅二乗コヒーレンス

ページ内をすべて折りたたむ

構文

cxy = mscohere(x,y)

cxy = mscohere(x,y,window)

cxy = mscohere(x,y,window,noverlap)

cxy = mscohere(x,y,window,noverlap,nfft)

cxy = mscohere(___,'mimo')

[cxy,w] = mscohere(___)

[cxy,f] = mscohere(___,fs)

[cxy,w] = mscohere(x,y,window,noverlap,w)

[cxy,f] = mscohere(x,y,window,noverlap,f,fs)

[___] = mscohere(x,y,___,freqrange)

mscohere(___)

説明

cxy = mscohere(x,y) は、入力信号 x と y の振幅二乗コヒーレンス推定 cxy を求めます。

x と y の両方がベクトルの場合、それらは同じ長さでなければなりません。
信号の 1 つが行列でもう 1 つがベクトルの場合は、ベクトルの長さが行列の行数と等しくなければなりません。関数はベクトルを展開し、列ごとの振幅二乗コヒーレンス推定の行列を返します。
x と y が同じ行数だが異なる列数をもつ行列の場合、mscohere は多重コヒーレンス行列を返します。cxy の m 番目の列には、すべての入力信号と m 番目の出力信号間の相関度の推定が含まれます。詳細については、振幅二乗コヒーレンスを参照してください。
x と y が等しいサイズの行列の場合、mscohere は列方向に動作します (cxy(:,n) = mscohere(x(:,n),y(:,n)))。多重コヒーレンス行列を取得するには、引数リストに 'mimo' を追加します。

cxy = mscohere(x,y,window) は、window を使用して x と y をセグメントに分割し、ウィンドウ処理を実行します。少なくとも 2 つのセグメントを使用しなければなりません。そうでない場合、すべての周波数の振幅二乗コヒーレンスが 1 になります。MIMO では、セグメント数が入力チャネル数より大きくなければなりません。

cxy = mscohere(x,y,window,noverlap) は、隣り合ったセグメント間で noverlap 個のサンプルのオーバーラップを使用します。

例

cxy = mscohere(x,y,window,noverlap,nfft) は、nfft サンプリング点を使用して離散フーリエ変換を計算します。

例

cxy = mscohere(___,'mimo') は、行列入力の多重コヒーレンス行列を計算します。この構文には、前の構文の入力引数を任意に組み合わせて含めることができます。

[cxy,w] = mscohere(___) は、正規化周波数のベクトル w を返し、これにより振幅二乗コヒーレンスが推定されます。

例

[cxy,f] = mscohere(___,fs) は、サンプルレート fs で表される周波数のベクトル f を返し、これにより、振幅二乗コヒーレンスが推定されます。fs は mscohere に対する 6 番目の数値入力でなければなりません。サンプルレートを入力した場合でも、前のオプション引数の既定値を使用するには、これらの引数を空 [] として指定します。

[cxy,w] = mscohere(x,y,window,noverlap,w) は、w で指定される正規化周波数における振幅二乗コヒーレンス推定を返します。

[cxy,f] = mscohere(x,y,window,noverlap,f,fs) は、f で指定される周波数における振幅二乗コヒーレンス推定を返します。

[___] = mscohere(x,y,___,freqrange) は、freqrange で指定される周波数範囲での振幅二乗コヒーレンス推定を返します。freqrange の有効なオプションは、'onesided'、'twosided' および 'centered' です。

例

出力引数を設定せずに mscohere(___) を使用すると、現在の Figure ウィンドウに振幅二乗コヒーレンス推定がプロットされます。

例

すべて折りたたむ

2 つのシーケンスのコヒーレンス推定

ライブスクリプトを開く

2 つのカラードノイズシーケンス間のコヒーレンス推定を計算して、プロットします。

ホワイトガウスノイズから成る信号を生成します。

r = randn(16384,1);

最初のシーケンスを作成するために、バンドパスフィルターで信号をフィルター処理します。0.2π ～ 0.4π ラジアン/サンプルの正規化周波数を通過させる 16 次フィルターを設計します。60 dB の阻止帯域の減衰量を指定します。元の信号をフィルター処理します。

dx = designfilt('bandpassiir','FilterOrder',16, ...
    'StopbandFrequency1',0.2,'StopbandFrequency2',0.4, ...
    'StopbandAttenuation',60);
x = filter(dx,r);

2 番目のシーケンスを作成するために、0.6π ～ 0.8π ラジアン/サンプルの正規化周波数を阻止する 16 次フィルターを設計します。0.1 dB の通過帯域リップルを指定します。元の信号をフィルター処理します。

dy = designfilt('bandstopiir','FilterOrder',16, ...
    'PassbandFrequency1',0.6,'PassbandFrequency2',0.8, ...
    'PassbandRipple',0.1);
y = filter(dy,r);

x と y の振幅二乗コヒーレンスを推定します。512 サンプルのハミングウィンドウを使用します。隣り合ったセグメント間のオーバーラップを 500 サンプル、DFT 点を 2048 に指定します。

[cxy,fc] = mscohere(x,y,hamming(512),500,2048);

コヒーレンス関数をプロットし、フィルターの周波数応答を重ね合わせます。

[qx,f] = freqz(dx);
qy = freqz(dy);

plot(fc/pi,cxy)
hold on
plot(f/pi,abs(qx),f/pi,abs(qy))
hold off

Figure contains an axes object. The axes object contains 3 objects of type line.

多重コヒーレンスと通常のコヒーレンス

ライブスクリプトを開く

ランダムな 2 チャネル信号 x を生成します。もう 1 つの信号 y を生成するために、2 チャネルをローパスフィルター処理してこれらを加算します。カットオフ周波数が 0.3π で、箱型ウィンドウを使用して設計された 30 次の FIR フィルターを指定します。

h = fir1(30,0.3,rectwin(31));
x = randn(16384,2);
y = sum(filter(h,1,x),2);

x と y の多重コヒーレンス推定を計算します。1024 サンプルのハンウィンドウを使用して信号にウィンドウを適用します。隣り合ったセグメント間のオーバーラップを 512 サンプル、DFT 点を 1024 に指定します。推定をプロットします。

noverlap = 512;
nfft = 1024;

mscohere(x,y,hann(nfft),noverlap,nfft,'mimo')

Figure contains an axes object. The axes object with title Coherence Estimate via Welch contains an object of type line.

コヒーレンス推定をフィルターの周波数応答と比較します。コヒーレンスの低下が周波数応答のゼロと対応します。

[H,f] = freqz(h);

hold on
yyaxis right
plot(f/pi,20*log10(abs(H)))
hold off

Figure contains an axes object. The axes object with title Coherence Estimate via Welch contains 2 objects of type line.

x と y の通常の振幅二乗コヒーレンス推定を計算してプロットします。どのチャネルでも推定は 1 に達しません。

figure
mscohere(x,y,hann(nfft),noverlap,nfft)

Figure contains an axes object. The axes object with title Coherence Estimate via Welch contains 2 objects of type line.

MIMO システムのコヒーレンス

ライブスクリプトを開く

1 kHz でそれぞれ 2 秒間サンプリングされた 2 つのマルチチャネル信号を生成します。最初の信号は入力で、周波数が 120 Hz、360 Hz、480 Hz の 3 つの正弦波から構成されます。2 番目の信号は出力で、周波数が 120 Hz と 360 Hz の 2 つの正弦波から構成されます。一方の正弦波は最初の信号より π/2 遅延します。他方の正弦波の遅延は π/4 です。両方の信号がホワイトガウスノイズに組み込まれます。

fs = 1000;
f = 120;
t = (0:1/fs:2-1/fs)';

inpt = sin(2*pi*f*[1 3 4].*t);
inpt = inpt+randn(size(inpt));
oupt = sin(2*pi*f*[1 3].*t-[pi/2 pi/4]);
oupt = oupt+randn(size(oupt));

すべての入力信号と各出力チャネル間の相関度を推定します。長さ 100 のハミングウィンドウを使用してデータにウィンドウを適用します。mscohere は、出力チャネルごとに 1 つのコヒーレンス関数を返します。コヒーレンス関数は、入力と出力で共有される周波数において最大値に達します。

[Cxy,f] = mscohere(inpt,oupt,hamming(100),[],[],fs,'mimo');

for k = 1:size(oupt,2)
    subplot(size(oupt,2),1,k)
    plot(f,Cxy(:,k))
    title(['Output ' int2str(k) ', All Inputs'])
end

Figure contains 2 axes objects. Axes object 1 with title Output 1, All Inputs contains an object of type line. Axes object 2 with title Output 2, All Inputs contains an object of type line.

入力信号と出力信号を切り替え、多重コヒーレンス関数を計算します。同じハミングウィンドウを使用します。480 Hz において入力と出力は相関しません。したがって、3 番目の相関関数にピークはありません。

[Cxy,f] = mscohere(oupt,inpt,hamming(100),[],[],fs,'mimo');

for k = 1:size(inpt,2)
    subplot(size(inpt,2),1,k)
    plot(f,Cxy(:,k))
    title(['Input ' int2str(k) ', All Outputs'])
end

Figure contains 3 axes objects. Axes object 1 with title Input 1, All Outputs contains an object of type line. Axes object 2 with title Input 2, All Outputs contains an object of type line. Axes object 3 with title Input 3, All Outputs contains an object of type line.

mscohere のプロット機能を使用して計算を繰り返します。

clf
mscohere(oupt,inpt,hamming(100),[],[],fs,'mimo')

Figure contains an axes object. The axes object with title Coherence Estimate via Welch contains 3 objects of type line.

2 番目の信号および最初の信号の最初の 2 つのチャネルの通常のコヒーレンス関数を計算します。オフピーク値は多重コヒーレンス関数とは異なります。

[Cxy,f] = mscohere(oupt,inpt(:,[1 2]),hamming(100),[],[],fs);
plot(f,Cxy)

Figure contains an axes object. The axes object contains 2 objects of type line.

位相差を求めるために、最大コヒーレンス点におけるクロススペクトルの角度を計算します。

Pxy = cpsd(oupt,inpt(:,[1 2]),hamming(100),[],[],fs);
[~,mxx] = max(Cxy);
for k = 1:2
    fprintf('Phase lag %d = %5.2f*pi\n',k,angle(Pxy(mxx(k),k))/pi)
end

Phase lag 1 = -0.51*pi
Phase lag 2 = -0.22*pi

振幅二乗コヒーレンスプロットの変更

ライブスクリプトを開く

1 kHz でそれぞれ 1 秒間サンプリングされた 2 つの正弦波信号を生成します。各正弦波の周波数は 250 Hz です。一方の信号の位相を他方より π/3 ラジアン遅らせます。両方の信号を単位分散のホワイトガウスノイズに組み込みます。

fs = 1000;
f = 250;
t = 0:1/fs:1-1/fs;
um = sin(2*pi*f*t)+rand(size(t));
un = sin(2*pi*f*t-pi/3)+rand(size(t));

mscohere を使用し、信号の振幅二乗コヒーレンスを計算してプロットします。

mscohere(um,un,[],[],[],fs)

Figure contains an axes object. The axes object with title Coherence Estimate via Welch contains an object of type line.

プロットのタイトル、x 軸のラベルおよび y 軸の範囲を変更します。

title('Magnitude-Squared Coherence')
xlabel('f (Hz)')
ylim([0 1.1])

Figure contains an axes object. The axes object with title Magnitude-Squared Coherence contains an object of type line.

gca を使用して、現在の座標軸へのハンドルを取得します。目盛りの位置を変更します。y 軸のラベルを削除します。

ax = gca;
ax.XTick = 0:250:500;
ax.YTick = 0:0.25:1;
ax.YLabel.String = [];

Figure contains an axes object. The axes object with title Magnitude-Squared Coherence contains an object of type line.

ハンドルの Children プロパティを呼び出して、プロットされたラインの色と幅を変更します。

ln = ax.Children;
ln.Color = [0.8 0 0];
ln.LineWidth = 1.5;

Figure contains an axes object. The axes object with title Magnitude-Squared Coherence contains an object of type line.

または、set と get を使用してラインプロパティを変更することもできます。

set(get(gca,'Children'),'Color',[0 0.4 0],'LineStyle','--','LineWidth',1)

Figure contains an axes object. The axes object with title Magnitude-Squared Coherence contains an object of type line.

入力引数

すべて折りたたむ

`x`, `y` — 入力信号
ベクトル | 行列

ベクトルまたは行列として指定される入力信号。

例: cos(pi/4*(0:159))+randn(1,160) は、ホワイトガウスノイズに含まれる正弦波を指定します。

データ型: single | double
複素数のサポート: あり

`window` — ウィンドウ
整数 | ベクトル | `[]`

ウィンドウ。整数、あるいは行ベクトルまたは列ベクトルとして指定します。window は信号をセグメントに分割するために使用します。

window が整数の場合、mscohere は x と y を長さ window のセグメントに分割し、各セグメントにその長さのハミングウィンドウを適用します。
window がベクトルの場合、mscohere は x と y をベクトルと同じ長さのセグメントに分割し、window を使用して各セグメントにウィンドウを適用します。

x と y の長さが noverlap 個のオーバーラップサンプルをもつ整数個のセグメントに厳密に分割できない場合、信号はそれに応じた長さで打ち切られます。

window を空として指定した場合、mscohere はハミングウィンドウを使用して、x と y を noverlap 個のオーバーラップサンプルをもつ 8 個のセグメントに分割します。

利用可能なウィンドウのリストについては、ウィンドウを参照してください。

例: hann(N+1) と (1-cos(2*pi*(0:N)'/N))/2 は、いずれも長さ N + 1 のハンウィンドウを指定します。

データ型: single | double

`noverlap` — オーバーラップするサンプル数
正の整数 | `[]`

オーバーラップするサンプル数。正の整数で指定します。

window がスカラーの場合、noverlap は window より小さくなければなりません。
window がベクトルの場合、noverlap は window の長さより小さくなければなりません。

noverlap を空として指定した場合、mscohere はセグメント間で 50% のオーバーラップが発生する数を使用します。セグメントの長さを指定していない場合、関数により noverlap が ⌊N/4.5⌋ に設定されます。ここで、N は入力および出力信号の長さです。

データ型: double | single

`nfft` — DFT 点の数
正の整数 | `[]`

正の整数として指定する DFT 点の数。nfft を空として指定した場合、mscohere によりこの引数が max(256,2^p) に設定されます。ここで p = ⌈log₂ N⌉ は入力信号の長さ N についての値です。

データ型: single | double

`fs` — サンプルレート
正のスカラー

サンプルレート。正のスカラーで指定します。サンプルレートは単位時間あたりのサンプル数です。時間の単位が秒の場合、サンプルレートの単位は Hz です。

`w` — 正規化周波数
ベクトル

正規化周波数。少なくとも 2 つの要素をもつ行ベクトルまたは列ベクトルとして指定します。正規化周波数の単位はラジアン/サンプルです。

例: w = [pi/4 pi/2]

データ型: double

`f` — 周波数
ベクトル

周波数。少なくとも 2 つの要素をもつ行ベクトルまたは列ベクトルとして指定します。周波数の単位は単位時間あたりのサイクルです。単位時間はサンプルレート fs で指定されます。fs の単位がサンプル/秒であれば、f の単位は Hz です。

例: fs = 1000; f = [100 200]

データ型: double

`freqrange` — 振幅二乗コヒーレンス推定の周波数範囲
`'onesided'` | `'twosided'` | `'centered'`

振幅二乗コヒーレンス推定の周波数範囲。'onesided'、'twosided' または 'centered' で指定します。既定値は、実数値信号の場合は 'onesided'、複素数値信号の場合は 'twosided' です。

'onesided' — 2 つの実数値の入力信号 x と y の間の振幅二乗コヒーレンス推定の片側推定を返します。nfft が偶数の場合、cxy は nfft/2 + 1 行で、計算区間は [0,π] ラジアン/サンプルです。nfft が奇数の場合、cxy は (nfft + 1)/2 行で、計算区間は [0,π) ラジアン/サンプルです。fs を指定すると、nfft が偶数の場合は、対応する計算区間は [0,fs/2] サイクル/単位時間で、nfft が奇数の場合は、[0,fs/2) サイクル/単位時間です。
'twosided' — 実数値または複素数値の 2 つの入力信号 x と y の間の振幅二乗コヒーレンス推定の両側推定を返します。この場合、cxy は nfft 行で、計算区間は [0,2π) ラジアン/サンプルです。fs を指定した場合、計算区間は [0,fs) サイクル/単位時間となります。
'centered' — 実数値または複素数値の 2 つの入力信号 x と y の間の振幅二乗コヒーレンス推定の中央に揃えた両側推定を返します。この場合、cxy は nfft 行で、偶数の nfft については区間 (–π,π] ラジアン/サンプルで、奇数の nfft については (–π,π) ラジアン/サンプルで計算されます。fs を指定すると、nfft が偶数の場合は、対応する計算区間は (–fs/2, fs/2] サイクル/単位時間で、nfft が奇数の場合は、(–fs/2, fs/2) サイクル/単位時間です。

出力引数

すべて折りたたむ

`cxy` — 振幅二乗コヒーレンス推定
ベクトル | 行列 | 3 次元配列

振幅二乗コヒーレンス推定。ベクトル、行列または 3 次元配列として返されます。

`w` — 正規化周波数
ベクトル

正規化周波数。実数値の列ベクトルとして返されます。

`f` — 周波数
ベクトル

周波数。実数値の列ベクトルとして返されます。

詳細

すべて折りたたむ

振幅二乗コヒーレンス

振幅二乗コヒーレンス推定は 0 ～ 1 の値をもつ周波数の関数です。これらの値は、各周波数において x と y の一致度合い示します。振幅二乗コヒーレンスは、パワースペクトル密度 P_xx(f) と P_yy(f)、および x と y のクロスパワースペクトル密度 P_xy(f) の関数です。

$C_{x y} (f) = \frac{{| P_{x y} (f) |}^{2}}{P_{x x} (f) P_{y y} (f)} .$

多入力/多出力システムの場合、多重コヒーレンス関数は以下のようになります。

$\begin{matrix} C_{X y_{i}} (f) = \frac{P_{X y_{i}}^{†} (f) P_{X X}^{- 1} (f) P_{X y_{i}} (f)}{P_{y_{i} y_{i}} (f)} \\ = [\begin{matrix} P_{x_{1} y_{i}}^{*} (f) & \dots & P_{x_{m} y_{i}}^{*} (f) \end{matrix}] {[\begin{matrix} P_{x_{1} x_{1}} (f) & P_{x_{1} x_{2}} (f) & \dots & P_{x_{1} x_{m}} (f) \\ P_{x_{2} x_{1}} (f) & P_{x_{2} x_{2}} (f) & \dots & P_{x_{2} x_{m}} (f) \\ ⋮ & ⋮ & ⋱ & ⋮ \\ P_{x_{m} x_{1}} (f) & P_{x_{m} x_{2}} (f) & \dots & P_{x_{m} x_{m}} (f) \end{matrix}]}^{- 1} [\begin{matrix} P_{x_{1} y_{i}} (f) \\ ⋮ \\ P_{x_{m} y_{i}} (f) \end{matrix}] \frac{1}{P_{y_{i} y_{i}} (f)} \end{matrix}$

i 番目の出力信号について、以下のようになります。

X は m 個の入力の配列に対応します。
P_{Xy_i} は、入力と y_i の間のクロスパワースペクトル密度の m 次元ベクトルです。
P_XX は、入力のパワースペクトル密度とクロスパワースペクトル密度の m 行 m 列の行列です。
P_{y_i}_{y_i} は出力のパワースペクトル密度です。
ダガー (†) は複素共役転置を意味します。

アルゴリズム

mscohere では、ウェルチのオーバーラップ平均ピリオドグラム法[3]、[5]を使用して、振幅二乗コヒーレンス関数[2]が推定されます。

参照

[1] Gómez González, A., J. Rodríguez, X. Sagartzazu, A. Schumacher, and I. Isasa. “Multiple Coherence Method in Time Domain for the Analysis of the Transmission Paths of Noise and Vibrations with Non-Stationary Signals.” Proceedings of the 2010 International Conference of Noise and Vibration Engineering, ISMA2010-USD2010. pp. 3927–3941.

[2] Kay, Steven M. Modern Spectral Estimation. Englewood Cliffs, NJ: Prentice-Hall, 1988.

[3] Rabiner, Lawrence R., and Bernard Gold. Theory and Application of Digital Signal Processing. Englewood Cliffs, NJ: Prentice-Hall, 1975.

[4] Stoica, Petre, and Randolph Moses. Spectral Analysis of Signals. Upper Saddle River, NJ: Prentice Hall, 2005.

[5] Welch, Peter D. “The Use of Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging Over Short, Modified Periodograms.” IEEE^® Transactions on Audio and Electroacoustics. Vol. AU-15, 1967, pp. 70–73.

拡張機能

C/C++ コード生成
MATLAB® Coder™ を使用して C および C++ コードを生成します。

GPU 配列
Parallel Computing Toolbox™ を使用してグラフィックス処理装置 (GPU) 上で実行することにより、コードを高速化します。

この関数は、GPU 配列を完全にサポートします。詳細については、GPU での MATLAB 関数の実行 (Parallel Computing Toolbox)を参照してください。

バージョン履歴

R2006a より前に導入

参考

cpsd | periodogram | pwelch | tfestimate

mscohere

構文

説明

例

2 つのシーケンスのコヒーレンス推定

多重コヒーレンスと通常のコヒーレンス

MIMO システムのコヒーレンス

振幅二乗コヒーレンス プロットの変更

入力引数

x, y — 入力信号 ベクトル | 行列

window — ウィンドウ 整数 | ベクトル | []

noverlap — オーバーラップするサンプル数 正の整数 | []

nfft — DFT 点の数 正の整数 | []

fs — サンプルレート 正のスカラー

w — 正規化周波数 ベクトル

f — 周波数 ベクトル

freqrange — 振幅二乗コヒーレンス推定の周波数範囲 'onesided' | 'twosided' | 'centered'

出力引数

cxy — 振幅二乗コヒーレンス推定 ベクトル | 行列 | 3 次元配列

w — 正規化周波数 ベクトル

f — 周波数 ベクトル

詳細

振幅二乗コヒーレンス

アルゴリズム

参照

拡張機能

C/C++ コード生成 MATLAB® Coder™ を使用して C および C++ コードを生成します。

GPU 配列 Parallel Computing Toolbox™ を使用してグラフィックス処理装置 (GPU) 上で実行することにより、コードを高速化します。

バージョン履歴

参考

トピック

振幅二乗コヒーレンスプロットの変更

`x`, `y` — 入力信号
ベクトル | 行列

`window` — ウィンドウ
整数 | ベクトル | `[]`

`noverlap` — オーバーラップするサンプル数
正の整数 | `[]`

`nfft` — DFT 点の数
正の整数 | `[]`

`fs` — サンプルレート
正のスカラー

`w` — 正規化周波数
ベクトル

`f` — 周波数
ベクトル

`freqrange` — 振幅二乗コヒーレンス推定の周波数範囲
`'onesided'` | `'twosided'` | `'centered'`

`cxy` — 振幅二乗コヒーレンス推定
ベクトル | 行列 | 3 次元配列

`w` — 正規化周波数
ベクトル

`f` — 周波数
ベクトル

C/C++ コード生成
MATLAB® Coder™ を使用して C および C++ コードを生成します。

GPU 配列
Parallel Computing Toolbox™ を使用してグラフィックス処理装置 (GPU) 上で実行することにより、コードを高速化します。