xspectrogram
短時間フーリエ変換を使用したクロス スペクトログラム
構文
説明
[___] = xspectrogram(___, は、freqrange)freqrange で指定される周波数範囲にわたるクロス スペクトログラムを返します。freqrange の有効なオプションは、"onesided"、"twosided" および "centered" です。
[___] = xspectrogram(___, は、名前と値の引数を使用して追加オプションを指定します。オプションには、最小しきい値と出力時間次元が含まれます。Name,Value)
[___] = xspectrogram(___, は、spectrumtype)spectrumtype を "psd" として指定した場合は、短時間のクロス パワー スペクトル密度推定を返し、spectrumtype を "power" として指定した場合は、短時間のクロス パワー スペクトル推定を返します。
出力引数なしで 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('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" という音節が特に顕著になります。
それぞれ 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 つのシーケンスのクロス スペクトログラムを計算し、プロットします。
各シーケンスを 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')
スペクトログラムに、フィルターによって強調または抑制された周波数領域が示されます。
入力引数
名前と値の引数
出力引数
参照
[1] Mitra, Sanjit K. Digital Signal Processing: A Computer-Based Approach. 2nd Ed. New York: McGraw-Hill, 2001.
[2] Oppenheim, Alan V., and Ronald W. Schafer, with John R. Buck. Discrete-Time Signal Processing. 2nd Ed. Upper Saddle River, NJ: Prentice Hall, 1999.
