このページの翻訳は最新ではありません。ここをクリックして、英語の最新版を参照してください。

cpsd

クロス パワー スペクトル密度

構文

pxy = cpsd(x,y)

pxy = cpsd(x,y,window)

pxy = cpsd(x,y,window,noverlap)

pxy = cpsd(x,y,window,noverlap,nfft)

pxy = cpsd(___,'mimo')

[pxy,w] = cpsd(___)

[pxy,f] = cpsd(___,fs)

[pxy,w] = cpsd(x,y,window,noverlap,w)

[pxy,f] = cpsd(x,y,window,noverlap,f,fs)

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

cpsd(___)

説明

pxy = cpsd(x,y) では、スペクトル推定の平均修正ピリオドグラム法であるウェルチ法を使用して、2 つの離散時間信号 x と y のクロス パワー スペクトル密度 (CPSD) が推定されます。

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

• 信号の 1 つが行列でもう 1 つがベクトルの場合は、ベクトルの長さが行列の行数と等しくなければなりません。関数はベクトルを展開し、列ごとのクロス パワー スペクトル密度推定の行列を返します。

• x と y が同じ行数だが異なる列数をもつ行列の場合、cpsd は、入力列のすべての組み合わせに対するクロス パワー スペクトル密度推定を含む 3 次元配列 pxy を返します。pxy の各列は x の列に対応し、各ページは y の列に対応します (pxy(:,m,n) = cpsd(x(:,m),y(:,n)))。

• x と y が等しいサイズの行列の場合、cpsd は列方向に動作します (pxy(:,n) = cpsd(x(:,n),y(:,n)))。多入力/多出力配列を取得するには、引数リストに 'mimo' を追加します。

x と y が実数の場合、cpsd は片側 CPSD を返します。x または y が複素数の場合は、cpsd は両側 CPSD を返します。

pxy = cpsd(x,y,window) は、window を使用して x と y をセグメントに分割し、ウィンドウ処理を実行します。

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

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

pxy = cpsd(___,'mimo') は、クロス パワー スペクトル密度推定の多入力/多出力配列を計算します。この構文には、前の構文の入力引数を任意に組み合わせて含めることができます。

[pxy,w] = cpsd(___) は、正規化周波数のベクトル w を返し、これによりクロス パワー スペクトル密度が推定されます。

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

[pxy,w] = cpsd(x,y,window,noverlap,w) は、w で指定した正規化周波数におけるクロス パワー スペクトル密度推定を返します。

[pxy,f] = cpsd(x,y,window,noverlap,f,fs) は、f で指定した周波数におけるクロス パワー スペクトル密度推定を返します。

[___] = cpsd(x,y,___,freqrange) は、freqrange で指定した周波数範囲におけるクロス パワー スペクトル密度推定を返します。freqrange の有効なオプションは、'onesided'、'twosided' および 'centered' です。

出力引数を設定せずに cpsd(___) を使用すると、現在の Figure ウィンドウにクロス パワー スペクトル密度推定がプロットされます。

例

カラード ノイズ信号のクロス パワー スペクトル密度

2 つのカラード ノイズ信号を生成し、そのクロス パワー スペクトル密度をプロットします。長さ 1024 の FFT と、500 の点からなるオーバーラップのない三角ウィンドウを指定します。

r = randn(16384,1);

hx = fir1(30,0.2,rectwin(31));
x = filter(hx,1,r);

hy = ones(1,10)/sqrt(10);
y = filter(hy,1,r);

cpsd(x,y,triang(500),250,1024)

SISO および MIMO クロス パワー スペクトル密度

1 kHz で 1 秒間サンプリングされた 2 つの 2 チャネルの正弦波を生成します。最初の信号のチャネルは 200 Hz と 300 Hz の周波数をもちます。2 番目の信号のチャネルは 300 Hz と 400 Hz の周波数をもちます。両方の信号を単位分散のホワイト ガウス ノイズに組み込みます。

fs = 1e3;
t = (0:1/fs:1-1/fs)';

q = 2*sin(2*pi*[200 300].*t);
q = q+randn(size(q));

r = 2*sin(2*pi*[300 400].*t);
r = r+randn(size(r));

2 つの信号のクロス パワー スペクトル密度を計算します。256 サンプルのバートレット ウィンドウを使用し、信号をセグメントに分割してセグメントにウィンドウを適用します。隣り合ったセグメント間のオーバーラップを 128 サンプル、DFT 点を 2048 に指定します。cpsd の組み込み機能を使用して、結果をプロットします。

cpsd(q,r,bartlett(256),128,2048,fs)

既定では、cpsd は同じサイズの行列入力に対して列ごとに処理を行います。各チャネルのピークは元の正弦波の周波数において出現します。

計算を繰り返しますが、今回は引数のリストに 'mimo' を追加します。

cpsd(q,r,bartlett(256),128,2048,fs,'mimo')

'mimo' オプションを指定して呼び出す場合、cpsd は、入力列のすべての組み合わせに対するクロス パワー スペクトル密度推定を含む 3 次元配列を返します。2 番目のチャネル q および最初のチャネル r の推定には、共通の周波数 300 Hz に強調されたピークが示されます。

遅れのある正弦波のクロス スペクトル位相

1 kHz で 296 ms 間サンプリングされた 100 Hz 正弦波信号を 2 つ生成します。一方の正弦波を他方より 2.5 ms 遅らせます。これは、π/2 の位相遅れと等価です。両方の信号を、分散 1/4² のホワイト ガウス ノイズに組み込みます。

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

x = cos(2*pi*t*100)+0.25*randn(size(t));
tau = 1/400;
y = cos(2*pi*100*(t-tau))+0.25*randn(size(t));

クロス パワー スペクトル密度の大きさを計算してプロットします。cpsd の既定の設定を使用します。大きさは、信号間で著しいコヒーレンスがある周波数でピークになります。

cpsd(x,y,[],[],[],Fs)

振幅二乗コヒーレンス関数およびクロス スペクトルの位相をプロットします。コヒーレンスの高い周波数の縦座標は正弦波間の位相遅れに対応しています。

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

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

subplot(2,1,1)
plot(F,Cxy)
title('Magnitude-Squared Coherence')

subplot(2,1,2)
plot(F,angle(Pxy))

hold on
plot(F,2*pi*100*tau*ones(size(F)),'--')
hold off

xlabel('Hz')
ylabel('\Theta(f)')
title('Cross Spectrum Phase')

指数シーケンスのクロス パワー スペクトル密度

2 つの $$N$$ サンプルの指数シーケンス、$$x_a=a^n$$ および $$x_b=b^n$$ ($$n\ge 0$$) を生成します。$$a=0.8$$、$$b=0.9$$、および小さい $$N$$ を指定して有限サイズの効果を確認します。

N = 10;
n = 0:N-1;

a = 0.8;
b = 0.9;

xa = a.^n;
xb = b.^n;

シーケンスのクロス パワー スペクトル密度を正規化周波数の完全な区間 $$[-\pi ,\pi ]$$ にわたって計算してプロットします。長さが $$N$$ で、セグメント間のオーバーラップがない箱型ウィンドウを指定します。

w = -pi:1/1000:pi;
wind = rectwin(N);
nove = 0;

[pxx,f] = cpsd(xa,xb,wind,nove,w);

2 つのシーケンスのクロス パワー スペクトルの解析式は、$$N$$ が大きい場合は以下のようになります。

$$R(\omega )=\frac{1}{1-a{e}^{-j\omega }}\frac{1}{1-b{e}^{j\omega }}.$$

この式をクロス パワー スペクトル密度に変換するには、$$2\pi N$$ で除算します。結果を比較します。cpsd の結果に含まれるリップルは、ウィンドウ処理の結果です。

nfac = 2*pi*N;

X = 1./(1-a*exp(-1j*w));
Y = 1./(1-b*exp( 1j*w));
R = X.*Y/nfac;

semilogy(f/pi,abs(pxx))
hold on
semilogy(w/pi,abs(R))
hold off
legend('cpsd','Analytic')

$$N=25$$ を使用して計算を繰り返します。$$N$$ を 100 まで小さくすると、曲線は 6 つの Figure と一致します。

N = 25;
n = 0:N-1;
xa = a.^n;
xb = b.^n;

wind = rectwin(N);

[pxx,f] = cpsd(xa,xb,wind,nove,w);
R = X.*Y/(2*pi*N);

semilogy(f/pi,abs(pxx))
hold on
semilogy(w/pi,abs(R))
hold off
legend('cpsd','Analytic')

ダイヤル トーン認識

クロス パワー スペクトル密度を使用して、大きく破損したトーンを識別します。

デジタル電話の番号またはシンボルをダイヤルしたときに生成される音声信号は、2 つの異なるグループから取得された周波数をもつ正弦波の和で構成されます。トーンの各ペアは、低群 (697 Hz、770 Hz、852 Hz または 941 Hz) の 1 つの周波数と高群 (1209 Hz、1336 Hz または 1477 Hz) の 1 つの周波数で構成されます。

すべてのシンボルに対応する信号を生成します。各トーンを 4 kHz で 0.5 秒間サンプリングします。参照表を準備します。

fs = 4e3;
t = 0:1/fs:0.5-1/fs;

nms = ['1';'2';'3';'4';'5';'6';'7';'8';'9';'*';'0';'#'];

ver = [697 770 852 941];
hor = [1209 1336 1477];

v = length(ver);
h = length(hor);

for k = 1:v
    for l = 1:h
        idx = h*(k-1)+l;
        tone = sum(sin(2*pi*[ver(k);hor(l)].*t))';
        tones(:,idx) = tone;
    end
end

各信号のウェルチ ピリオドグラムをプロットし、周波数成分に注釈を付けます。200 サンプルのハミング ウィンドウを使用し、信号をオーバーラップしないセグメントに分割してセグメントにウィンドウを適用します。

[pxx,f] = pwelch(tones,hamming(200),0,[],fs);

for k = 1:v
    for l = 1:h
        idx = h*(k-1)+l;
        ax = subplot(v,h,idx);
        plot(f,10*log10(pxx(:,idx)))
        ylim([-80 0])
        title(nms(idx))
        tx = [ver(k);hor(l)];
        ax.XTick = tx;
        ax.XTickLabel = int2str(tx);
    end
end

番号 8 をダイヤルすると生成される信号は、ノイズの多いチャネルを通して送信されます。受信信号は大きく破損しているため、調べても番号を識別できません。

mys = sum(sin(2*pi*[ver(3);hor(2)].*t))'+5*randn(size(t'));

% To hear, type soundsc(mys,fs)

破損した信号と基準信号のクロス パワー スペクトル密度を計算します。512 サンプルのカイザー ウィンドウを使用し、形状係数を β = 5 として、信号にウィンドウを適用します。各スペクトルの振幅をプロットします。

[pxy,f] = cpsd(mys,tones,kaiser(512,5),100,[],fs);

for k = 1:v
    for l = 1:h
        idx = h*(k-1)+l;
        ax = subplot(v,h,idx);
        plot(f,10*log10(abs(pxy(:,idx))))
        ylim([-80 0])
        title(nms(idx))
        tx = [ver(k);hor(l)];
        ax.XTick = tx;
        ax.XTickLabel = int2str(tx);
    end
end

破損した信号の数字は、最も高いピークと最も高い RMS 値のスペクトルをもちます。

[~,loc] = max(rms(abs(pxy)));

digit = nms(loc)

digit = 
'8'

入力引数

x, y — 入力信号

ベクトル | 行列

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

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

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

window — ウィンドウ

整数 | ベクトル | []

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

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

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

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

window を空として指定した場合、cpsd はハミング ウィンドウを使用して、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 を空として指定した場合、cpsd はセグメント間で 50% のオーバーラップが発生する数を使用します。セグメントの長さを指定していない場合、関数により noverlap が ⌊N/4.5⌋ に設定されます。ここで、N は入力および出力信号の長さです。

データ型: double | single

nfft — DFT 点の数

正の整数 | []

正の整数として指定する DFT 点の数。nfft を空として指定した場合、cpsd によりパラメーターが 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 が偶数の場合、pxy は nfft/2 + 1 行で、計算区間は [0,π] ラジアン/サンプルです。nfft が奇数の場合、pxy は (nfft + 1)/2 行で、計算区間は [0,π) ラジアン/サンプルです。fs を指定すると、nfft が偶数の場合は、対応する計算区間は [0,fs/2] サイクル/単位時間で、nfft が奇数の場合は、[0,fs/2) サイクル/単位時間です。

• 'twosided' — 2 つの実数値または複素数値の入力信号 x と y のクロス パワー スペクトル密度の両側推定を返します。この場合、pxy は nfft 行で、計算区間は [0,2π) ラジアン/サンプルです。fs を指定した場合、計算区間は [0,fs) サイクル/単位時間となります。

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

出力引数

pxy — クロス パワー スペクトル密度

ベクトル | 行列 | 3 次元配列

クロス パワー スペクトル密度。ベクトル、行列または 3 次元配列として返されます。

w — 正規化周波数

ベクトル

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

• pxy が片側の場合、nfft が偶数のときは、w は区間 [0,π] をカバーし、nfft が奇数のときは、[0,π) をカバーします。

• pxy が両側の場合、w は区間 [0,2π) をカバーします。

• pxy が DC を中央に揃えたものの場合、nfft が偶数のときは、w は区間 (–π,π] をカバーし、nfft が奇数のときは、(–π,π) をカバーします。

データ型: double | single

f — 周波数

ベクトル

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

データ型: double | single

詳細

クロス パワー スペクトル密度

クロス パワー スペクトル密度は、単位周波数あたりのパワーの分布で、以下のように定義されます。

$$P_{xy}(\omega )={\displaystyle \sum _{m=-\infty }^{\infty }{R}_{xy}(m)e^{-j\omega m}.}$$

相互相関列は、以下のように定義されます。

$$R_{xy}(m)=E\left\{x_{n+m}{y}_n^{\ast }\right\}=E\left\{x_n{y}_{n-m}^{\ast }\right\},$$

ここで、x_n と y_n は共に定常的なランダム過程で、–∞ < n < ∞、$$-\infty <n<\infty$$ であり、E {· } は期待値演算子です。

アルゴリズム

cpsd では、スペクトル推定にウェルチの平均修正ピリオドグラム法が使用されます。

参照

[1] Oppenheim, Alan V., Ronald W. Schafer, and John R. Buck. Discrete-Time Signal Processing. 2nd Ed. Upper Saddle River, NJ: Prentice Hall, 1999.

[2] Rabiner, Lawrence R., and B. Gold. Theory and Application of Digital Signal Processing. Englewood Cliffs, NJ: Prentice-Hall, 1975, pp. 414–419.

[3] Welch, Peter D. “The Use of the 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, June 1967, pp. 70–73.

拡張機能

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

GPU 配列

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

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

バージョン履歴

R2006a より前に導入

参考

mscohere | pburg | pcov | peig | periodogram | pmcov | pmtm | pmusic | pwelch | pyulear | tfestimate

トピック

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

×

例を開く

この例の変更されたバージョンがあります。編集された方の例を開きますか?

いいえ、変更されたバージョンを上書きします。 はい