ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

mscohere

振幅二乗コヒーレンス

説明

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

  • xy の両方がベクトルの場合、それらは同じ長さでなければなりません。

  • 信号の 1 つが行列でもう 1 つがベクトルの場合は、ベクトルの長さが行列の行数と等しくなければなりません。関数はベクトルを展開し、列ごとの振幅二乗コヒーレンス推定の行列を返します。

  • xy が同じ行数だが異なる列数をもつ行列の場合、mscohere は多重コヒーレンス行列を返します。cxy の m 番目の列には、すべての入力信号と m 番目の出力信号間の相関度の推定が含まれます。詳細は、振幅二乗コヒーレンスを参照してください。

  • xy が等しいサイズの行列の場合、mscohere は列方向に動作します (cxy(:,n) = mscohere(x(:,n),y(:,n)))。多重コヒーレンス行列を取得するには、引数リストに 'mimo' を追加します。

cxy = mscohere(x,y,window) は、window を使用して xy をセグメントに分割し、ウィンドウ処理を実行します。少なくとも 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 を返し、これにより、振幅二乗コヒーレンスが推定されます。fsmscohere に対する 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 つのカラード ノイズ シーケンス間のコヒーレンス推定を計算して、プロットします。

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

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);

xy の振幅二乗コヒーレンスを推定します。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

ランダムな 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);

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

noverlap = 512;
nfft = 1024;

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

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

[H,f] = freqz(h);

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

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

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

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

入力信号と出力信号を切り替え、多重コヒーレンス関数を計算します。同じハミング ウィンドウを使用します。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

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

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

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

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

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

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)

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

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

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

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

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

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

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

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

入力引数

すべて折りたたむ

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

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

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

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

  • window が整数の場合、mscoherexy を長さ window のセグメントに分割し、各セグメントにその長さのハミング ウィンドウを適用します。

  • window がベクトルの場合、mscoherexy をベクトルと同じ長さのセグメントに分割し、window を使用して各セグメントにウィンドウを適用します。

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

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

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

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

データ型: single | double

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

  • window がスカラーの場合、noverlapwindow より小さくなければなりません。

  • window がベクトルの場合、noverlapwindow の長さより小さくなければなりません。

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

データ型: double | single

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

データ型: single | double

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

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

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

データ型: double

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

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

データ型: double

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

  • 'onesided' — 2 つの実数値の入力信号 xy の間の振幅二乗コヒーレンス推定の片側推定を返します。nfft が偶数の場合、cxynfft/2 + 1 行で、計算区間は [0,π] ラジアン/サンプルです。nfft が奇数の場合、cxy は (nfft + 1)/2 行で、計算区間は [0,π) ラジアン/サンプルです。fs を指定すると、nfft が偶数の場合は、対応する計算区間は [0,fs/2] サイクル/単位時間で、nfft が奇数の場合は、[0,fs/2) サイクル/単位時間です。

  • 'twosided' — 実数値または複素数値の 2 つの入力信号 xy の間の振幅二乗コヒーレンス推定の両側推定を返します。この場合、cxynfft 行で、計算区間は [0,2π) ラジアン/サンプルです。fs を指定した場合、計算区間は [0,fs) サイクル/単位時間となります。

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

出力引数

すべて折りたたむ

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

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

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

詳細

すべて折りたたむ

振幅二乗コヒーレンス

振幅二乗コヒーレンス推定は 0 ~ 1 の値をもつ周波数の関数です。これらの値は、各周波数において xy の一致度合い示します。振幅二乗コヒーレンスは、パワー スペクトル密度 Pxx(f)Pyy(f) の関数で、xy:

Cxy(f)=|Pxy(f)|2Pxx(f)Pyy(f).

のクロス パワー スペクトル密度 Pxy(f) の関数です。

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

CXyi(f)=PXyi(f)PXX1(f)PXyi(f)Pyiyi(f)=[Px1yi*(f)Pxmyi*(f)][Px1x1(f)Px1x2(f)Px1xm(f)Px2x1(f)Px2x2(f)Px2xm(f)Pxmx1(f)Pxmx2(f)Pxmxm(f)]1[Px1yi(f)Pxmyi(f)]1Pyiyi(f)

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

  • X は m 個の入力の配列に対応します。

  • PXyi は、入力と yi の間のクロス パワー スペクトル密度の m 次元ベクトルです。

  • PXX は、入力のパワー スペクトル密度とクロス パワー スペクトル密度の m 行 m 列の行列です。

  • Pyiyi は出力のパワー スペクトル密度です。

  • ダガー (†) は複素共役転置を意味します。

アルゴリズム

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++ コードを生成します。

R2006a より前に導入