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

xspectrogram

短時間フーリエ変換を使用したクロス スペクトログラム

構文

s = xspectrogram(x,y)

s = xspectrogram(x,y,window)

s = xspectrogram(x,y,window,noverlap)

s = xspectrogram(x,y,window,noverlap,nfft)

[s,w,t] = xspectrogram(___)

[s,f,t] = xspectrogram(___,fs)

[s,w,t] = xspectrogram(x,y,window,noverlap,w)

[s,f,t] = xspectrogram(x,y,window,noverlap,f,fs)

[___,c] = xspectrogram(___)

[___] = xspectrogram(___,freqrange)

[___] = xspectrogram(___,spectrumtype)

[___] = xspectrogram(___,'MinThreshold',thresh)

xspectrogram(___)

xspectrogram(___,freqloc)

説明

s = xspectrogram(x,y) は、x と y で指定された信号のクロス スペクトログラムを返します。入力信号は同数の要素をもつベクトルでなければなりません。s の各列は、x と y に共通の短時間の、時間が局所化された周波数成分の推定を含みます。

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

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

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

[s,w,t] = xspectrogram(___) は、正規化周波数 w のベクトルとクロス スペクトログラムを計算した時点 t のベクトルを返します。この構文には、前の構文の入力引数を任意に組み合わせて含めることができます。

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

[s,w,t] = xspectrogram(x,y,window,noverlap,w) は、w で指定した正規化周波数におけるクロス スペクトログラムを返します。

[s,f,t] = xspectrogram(x,y,window,noverlap,f,fs) は、f で指定した周波数におけるクロス スペクトログラムを返します。

[___,c] = xspectrogram(___) は、入力信号の時変複素クロス スペクトルの推定を含む行列 c も返します。クロス スペクトログラム s は c の振幅です。

[___] = xspectrogram(___,freqrange) は、freqrange で指定される周波数範囲にわたるクロス スペクトログラムを返します。freqrange の有効なオプションは、'onesided'、'twosided' および 'centered' です。

[___] = xspectrogram(___,spectrumtype) は、spectrumtype を 'psd' に指定した場合は、短時間のクロス パワー スペクトル密度推定を返し、spectrumtype を 'power' に指定した場合は、短時間のクロス パワー スペクトル推定を返します。

[___] = xspectrogram(___,'MinThreshold',thresh) は、10 log₁₀(s) ≤ thresh となる s の要素をゼロに設定します。thresh をデシベル単位で指定します。

出力引数なしで xspectrogram(___) を使用すると、現在の Figure ウィンドウにクロス スペクトログラムがプロットされます。

xspectrogram(___,freqloc) では、周波数をプロットする軸を指定します。freqloc は 'xaxis' または 'yaxis' のいずれかに指定します。

例

線形チャープのクロス スペクトログラム

1 MHz で 10 ms 間サンプリングした 2 つの線形チャープを生成します。

• 最初のチャープの初期周波数は 150 kHz で、測定の最後には 350 kHz に増加します。

• 2 番目のチャープの初期周波数は 200 kHz で、測定の最後には 300 kHz に増加します。

S/N 比が 40 dB となるホワイト ガウス ノイズを付加します。

nSamp = 10000;
Fs = 1000e3;
SNR = 40;
t = (0:nSamp-1)'/Fs;

x1 = chirp(t,150e3,t(end),350e3);
x1 = x1+randn(size(x1))*std(x1)/db2mag(SNR);

x2 = chirp(t,200e3,t(end),300e3);
x2 = x2+randn(size(x2))*std(x2)/db2mag(SNR);

2 つのチャープのクロス スペクトログラムを計算し、プロットします。信号を 200 サンプルのセグメントに分割し、各セグメントにハミング ウィンドウを適用します。隣接するセグメント間の 80 サンプルのオーバーラップ、および 1024 サンプルの DFT 長を指定します。

xspectrogram(x1,x2,hamming(200),80,1024,Fs,'yaxis')

2 番目のチャープを変更し、測定中に周波数が 50 kHz から 350 kHz に増加するようにします。形状係数 で 500 サンプルのカイザー ウィンドウを使用して、セグメントにウィンドウを適用します。450 サンプルのオーバーラップおよび 256 の DFT 長を指定します。クロス スペクトログラムを計算してプロットします。

x2 = chirp(t,50e3,t(end),350e3);
x2 = x2+randn(size(x2))*std(x2)/db2mag(SNR);

xspectrogram(x1,x2,kaiser(500,5),450,256,Fs,'yaxis')

両方の場合について、関数は 2 つの信号が共通にもつ周波数成分を強調表示します。

音声信号のクロス スペクトログラム

44,100 Hz でサンプリングされた 2 つの音声信号を含むファイルを読み込みます。

• 最初の信号は "transform function" という女性の声の録音です。

• 2 番目の信号は "reform justice" という同じ女性の声の録音です。

2 つの信号をプロットします。2 番目の信号を垂直方向にオフセットして、両方が表示されるようにします。

load(fullfile(matlabroot,'examples','signal','voice.mat'))

% To hear, type soundsc(transform,fs),pause(2),soundsc(reform,fs)

t = (0:length(reform)-1)/fs;

plot(t,transform,t,reform+0.3)
legend('"Transform function"','"Reform justice"')

2 つの信号のクロス スペクトログラムを計算します。信号を 1000 サンプルのセグメントに分割し、ハミング ウィンドウを適用します。隣接するセグメント間に 800 個のサンプルのオーバーラップを指定します。4 kHz までの周波数のみを含めます。

nwin = 1000;
nvlp = 800;
fint = 0:4000;

[s,f,t] = xspectrogram(transform,reform,hamming(nwin),nvlp,fint,fs);

mesh(t,f,20*log10(s))
view(2)
axis tight

クロス スペクトログラムは、信号がより多くの周波数成分を共通してもつ時間間隔を強調表示します。"form" という音節が特に顕著になります。

2 つの 2 次チャープ間の位相シフト

それぞれ 1 kHz で 2 秒間サンプリングされた 2 つの 2 次チャープを生成します。両方のチャープの初期周波数は 100 Hz で、測定の半ばには 200 Hz に増加します。2 番目のチャープの最初のチャープに対する位相差は 23° です。

fs = 1e3;
t = 0:1/fs:2;

y1 = chirp(t,100,1,200,'quadratic',0);
y2 = chirp(t,100,1,200,'quadratic',23);

チャープの複素クロス スペクトログラムを計算し、その間の位相シフトを抽出します。信号を 128 サンプルのセグメントに分割します。隣接するセグメント間に 120 個のサンプルのオーバーラップを指定します。形状係数が β = 18 のカイザー ウィンドウを使用して各セグメントをウィンドウ処理し、128 サンプルの DFT 長を指定します。xspectrogram のプロット機能を使用して、クロス スペクトログラムを表示します。

[~,f,xt,c] = xspectrogram(y1,y2,kaiser(128,18),120,128,fs);

xspectrogram(y1,y2,kaiser(128,18),120,128,fs,'yaxis')

クロス スペクトログラムの最大エネルギーの時間-周波数リッジを抽出して表示します。

[tfr,~,lridge] = tfridge(c,f);

hold on
plot(xt,tfr,'k','linewidth',2)
hold off

位相シフトは、リッジに沿った時変クロス スペクトルの虚数部と実数部の比です。位相シフトを計算して、度単位で表します。平均値を表示します。

pshft = angle(c(lridge))*180/pi;

mean(pshft)

ans = -23.0000

複素信号のクロス スペクトログラム

3 kHz でそれぞれ 1 秒間サンプリングされた 2 つの信号を生成します。最初の信号は 2 次チャープで、測定中に周波数が 300 Hz から 1300 Hz に増加します。チャープはホワイト ガウス ノイズに組み込まれます。2 番目の信号もホワイト ノイズに組み込まれますが、これは正弦関数的に変化する周波数成分をもつチャープです。

fs = 3000;
t = 0:1/fs:1-1/fs;

x1 = chirp(t,300,t(end),1300,'quadratic')+randn(size(t))/100;

x2 = exp(2j*pi*100*cos(2*pi*2*t))+randn(size(t))/100;

2 つの信号のクロス スペクトログラムを計算し、プロットします。隣接するセグメント間で 255 サンプルがオーバーラップする 256 サンプルのセグメントに信号を分割します。形状係数が β = 30 のカイザー ウィンドウを使用して、セグメントにウィンドウを適用します。既定の DFT の点の数を使用します。クロス スペクトログラムの中央を周波数ゼロに揃えます。

nwin = 256;

xspectrogram(x1,x2,kaiser(nwin,30),nwin-1,[],fs,'centered','yaxis')

パワー スペクトル密度の代わりにパワー スペクトルを計算します。–40 dB より小さい値をゼロに設定します。プロットの中央をナイキスト周波数に揃えます。

xspectrogram(x1,x2,kaiser(nwin,30),nwin-1,[],fs, ...
    'power','MinThreshold',-40,'yaxis')
title('Cross-Spectrogram of Quadratic Chirp and Complex Chirp')

しきい値処理により、共通周波数の領域がさらに強調表示されます。

2 つのシーケンスのクロス スペクトログラム

2 つのシーケンスのクロス スペクトログラムを計算し、プロットします。

各シーケンスを 4096 サンプルの長さに指定します。

N = 4096;

最初のシーケンスを作成するために、ホワイト ガウス ノイズに組み込まれる凸 2 次チャープを生成し、バンドパス フィルターを適用します。

• チャープの初期正規化周波数は 0.1π で、測定の最後には 0.8π に増加します。

• 16 次フィルターは 0.2π ～ 0.4π ラジアン/サンプルの正規化周波数を通過させ、60 dB の阻止帯域の減衰量をもちます。

rx = chirp(0:N-1,0.1/2,N,0.8/2,'quadratic',[],'convex')'+randn(N,1)/100;
dx = designfilt('bandpassiir','FilterOrder',16, ...
    'StopbandFrequency1',0.2,'StopbandFrequency2',0.4, ...
    'StopbandAttenuation',60);
x = filter(dx,rx);

2 番目のシーケンスを作成するために、ホワイト ガウス ノイズに組み込まれる線形チャープを生成し、バンドストップ フィルターを適用します。

• チャープの初期正規化周波数は 0.9π で、測定の最後には 0.1π に減少します。

• 16 次フィルターは 0.6π ～ 0.8π ラジアン/サンプルの正規化周波数を阻止し、1 dB の通過帯域リップルをもちます。

ry = chirp(0:N-1,0.9/2,N,0.1/2)'+randn(N,1)/100;
dy = designfilt('bandstopiir','FilterOrder',16, ...
    'PassbandFrequency1',0.6,'PassbandFrequency2',0.8, ...
    'PassbandRipple',1);
y = filter(dy,ry);

2 つのシーケンスをプロットします。2 番目のシーケンスを垂直方向にオフセットして、両方が表示されるようにします。

plot([x y+2])

x と y のクロス スペクトログラムを計算し、プロットします。512 サンプルのハミング ウィンドウを使用します。隣り合ったセグメント間のオーバーラップを 500 サンプル、DFT 点を 2048 に指定します。

xspectrogram(x,y,hamming(512),500,2048,'yaxis')

–50 dB より小さいクロス スペクトログラム値をゼロに設定します。

xspectrogram(x,y,hamming(512),500,2048,'MinThreshold',-50,'yaxis')

スペクトログラムに、フィルターによって強調または抑制された周波数領域が示されます。

入力引数

x, y — 入力信号

ベクトル

入力信号。ベクトルで指定します。

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

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

window — ウィンドウ

整数 | ベクトル | []

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

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

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

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

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

データ型: double | single

nfft — DFT 点の数

正の整数 | []

DFT 点の数。正の整数スカラーとして指定します。nfft を空として指定した場合、xspectrogram により DFT 長が max(256,2^p) に設定されます。ここで p = ⌈log₂ N_w⌉ で、

• window がスカラーの場合は、N_w = window です。

• window がベクトルの場合は、N_w = length(window) です。

データ型: single | double

w — 正規化周波数

ベクトル

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

例: pi./[2 4]

データ型: double | single

f — 周波数

ベクトル

周波数。ベクトルで指定します。f は少なくとも 2 つの要素をもたなければなりません。f の単位はサンプルレート fs により指定されます。

データ型: double | single

fs — サンプルレート

1 Hz (既定値) | 正のスカラー

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

データ型: double | single

freqrange — クロス スペクトル推定の周波数範囲

'onesided' | 'twosided' | 'centered'

クロス スペクトル推定の周波数範囲。'onesided'、'twosided' または 'centered' として指定します。実数値信号の場合、既定の設定は 'onesided' です。複素数値信号の場合、既定の設定は 'twosided' で、'onesided' を指定するとエラーになります。

• 'onesided' — 実数の入力信号の片側クロス スペクトログラムを返します。nfft が偶数の場合、s は nfft/2 + 1 行をもち、計算区間は [0, π] ラジアン/サンプルです。nfft が奇数の場合、s は (nfft + 1)/2 行をもち、計算区間は [0, π) ラジアン/サンプルです。fs を指定すると、それぞれの場合の計算区間は [0, fs/2] サイクル/単位時間、[0, fs/2) サイクル/単位時間となります。

• 'twosided' — 実信号または複素信号の両側クロス スペクトログラムを返します。s は nfft 行をもち、計算区間は [0, 2π) ラジアン/サンプルです。fs を指定した場合、計算区間は [0, fs) サイクル/単位時間となります。

• 'centered' — 中央に揃えた、実信号または複素信号の両側クロス スペクトログラムを返します。s は nfft 行をもちます。nfft が偶数の場合、s の計算区間は (–π, π] ラジアン/サンプルです。nfft が奇数の場合、s の計算区間は (–π, π) ラジアン/サンプルです。fs を指定すると、それぞれの場合の計算区間は (–fs/2, fs/2] サイクル/単位時間、(–fs/2, fs/2) サイクル/単位時間となります。

spectrumtype — クロス パワー スペクトルのスケーリング

'psd' (既定値) | 'power'

クロス パワー スペクトルのスケーリング。'psd' または 'power' として指定します。

• spectrumtype を省略するか、'psd' を指定すると、クロス パワー スペクトル密度が返されます。

• 'power' を指定すると、クロス パワー スペクトル密度の各推定を分解能帯域幅でスケーリングします。これはウィンドウの等価ノイズ帯域幅とセグメント持続時間に依存します。結果は、各周波数のパワーの推定です。

thresh — しきい値

-Inf (既定値) | 実数スカラー

しきい値。デシベル単位で表される実数スカラーで指定します。xspectrogram は、10 log₁₀(s) ≤ thresh となる s の要素をゼロに設定します。

freqloc — 周波数の表示軸

'xaxis' (既定値) | 'yaxis'

周波数の表示軸。'xaxis' または 'yaxis' で指定します。

• 'xaxis' — 周波数が x 軸に、時間が y 軸に表示されます。

• 'yaxis' — 周波数が y 軸に、時間が x 軸に表示されます。

この引数は出力引数で xspectrogram を呼び出している場合に無視されます。

出力引数

s — クロス スペクトログラム

行列

クロス スペクトログラム。行列として返されます。時間は s の列方向に、周波数は行方向に下がって 0 から増加します。

• 入力信号 x と y が長さ N の場合、s は k 列をもちます。ここで、以下のようになります。

  • window がスカラーの場合は、k = ⌊(N – noverlap)/(window – noverlap)⌋。

  • window がベクトルの場合は、k = ⌊(N – noverlap)/(length(window) – noverlap)⌋。

• 入力信号が実数で nfft が偶数の場合は、s は (nfft/2 + 1) 行になります。

• 入力信号が実数で nfft が奇数の場合は、s は (nfft + 1)/2 行になります。

• 入力信号が複素数の場合は、s は nfft 行になります。

データ型: double | single

w — 正規化周波数

ベクトル

正規化周波数。ベクトルとして返されます。w の長さは s の行数と等しくなります。

データ型: double | single

t — 時点

ベクトル

時点。ベクトルとして返されます。t の時間値は、window を使用して指定された各セグメントの中間点に対応します。

データ型: double | single

f — 巡回周波数

ベクトル

巡回周波数。ベクトルとして返されます。f の長さは s の行数と等しくなります。

データ型: double | single

c — 時変複素クロス スペクトル

行列

時変複素クロス スペクトル。行列として返されます。クロス スペクトログラム s は c の振幅です。

データ型: double | single

参照

[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] Mitra, Sanjit K. Digital Signal Processing: A Computer-Based Approach. 2nd Ed. New York: McGraw-Hill, 2001.

参考

cpsd | mscohere | spectrogram

R2017a で導入