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 の点からなり、セグメント間のオーバーラップが 50% である三角ウィンドウを指定します。

rng default
r = randn(2e4,1);

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

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

win = triang(500);
nov = 250;

cpsd(x,y,win,nov,1024)

Figure contains an axes object. The axes object with title Welch Cross Power Spectral Density Estimate, xlabel Normalized Frequency ( times pi blank rad/sample), ylabel Power/frequency (dB/(rad/sample)) contains an object of type line.

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)

Figure contains an axes object. The axes object with title Welch Cross Power Spectral Density Estimate, xlabel Frequency (Hz), ylabel Power/frequency (dB/Hz) contains 2 objects of type line.

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

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

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

Figure contains an axes object. The axes object with title Welch Cross Power Spectral Density Estimate, xlabel Frequency (Hz), ylabel Power/frequency (dB/Hz) contains 4 objects of type line.

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

Figure contains an axes object. The axes object with title Welch Cross Power Spectral Density Estimate, xlabel Frequency (Hz), ylabel Power/frequency (dB/Hz) contains an object of type line.

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

[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')

$Figure contains 2 axes objects. Axes object 1 with title Magnitude-Squared Coherence contains an object of type line. Axes object 2 with title Cross Spectrum Phase, xlabel Hz, ylabel \Theta(f) contains 2 objects of type line.$

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

2 つの $N$ サンプルの指数シーケンス、 $x_{a} = a^{n}$ および $x_{b} = b^{n}$ ( $n \geq 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;

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

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

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

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

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

この式をクロスパワースペクトル密度に変換するには、 $2 π 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')

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent 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')

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent 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 | single

`f` — 周波数
ベクトル

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

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

データ型: double | single

`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

詳細