ドキュメンテーション

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

spectrogram

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

構文

  • s = spectrogram(x)
  • s = spectrogram(x,window)
  • s = spectrogram(x,window,noverlap)
  • s = spectrogram(x,window,noverlap,nfft)
  • [s,w,t] = spectrogram(___)
  • [s,f,t] = spectrogram(___,fs)
  • [s,w,t] = spectrogram(x,window,noverlap,w)
  • [s,f,t] = spectrogram(x,window,noverlap,f,fs)
  • [___,ps] = spectrogram(___)
  • [___] = spectrogram(___,'reassigned')
  • [___,ps,fc,tc] = spectrogram(___)
  • [___] = spectrogram(___,freqrange)
  • [___] = spectrogram(___,spectrumtype)
  • [___] = spectrogram(___,'MinThreshold',thresh)
  • spectrogram(___)
  • spectrogram(___,freqloc)

説明

s = spectrogram(x) は、入力信号 x の短時間フーリエ変換を返します。s の各列は、x の短時間の、時間が局在化された周波数成分の推定を含みます。

s = spectrogram(x,window) は、window を使用して信号をセクションに分割し、ウィンドウ処理を実行します。

s = spectrogram(x,window,noverlap) は、隣接するセクション間で noverlap サンプルのオーバーラップを使用します。

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

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

[s,f,t] = spectrogram(___,fs) は、巡回周波数 f のベクトルをサンプルレート fs で表して返します。

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

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

[___,ps] = spectrogram(___) は、各セクションのパワー スペクトル密度 (PSD) 推定またはパワー スペクトル推定を含む行列 ps も返します。

[___] = spectrogram(___,'reassigned') は、各 PSD またはパワー スペクトル推定を、そのエネルギーの中心の位置に再代入します。十分に局在化された時相成分またはスペクトル成分が信号に含まれる場合、このオプションは、よりシャープなスペクトログラムを生成します。

[___,ps,fc,tc] = spectrogram(___) 各 PSD またはパワー スペクトル推定のエネルギー中心における周波数および時間を含む 2 つの行列 fc および tc も返します。

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

[___] = spectrogram(___,spectrumtype) では、spectrumtype'psd' に指定されている場合は PSD 推定が、spectrumtype'power' に指定されている場合はパワー スペクトル推定が返されます。

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

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

spectrogram(___,freqloc) では、周波数をプロットする軸を指定します。

すべて折りたたむ

正弦波の和から構成される信号のサンプルを $N_{\tt x}=1024$ 個生成します。正弦波の正規化周波数は、 $2\pi/5$ ラジアン/サンプルおよび $4\pi/5$ ラジアン/サンプルです。高周波数の正弦波の振幅は、他の正弦波の振幅の 10 倍です。

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

w0 = 2*pi/5;
x = sin(w0*n)+10*sin(2*w0*n);

関数の既定の設定を使用して、短時間フーリエ変換を計算します。スペクトログラムをプロットします。

s = spectrogram(x);

spectrogram(x,'yaxis')

計算を繰り返します。

  • 長さ ${\tt nsc}=\lfloor N_{\tt x}/4.5\rfloor$ のセクションに信号を分割します。

  • ハミング ウィンドウを使用して、セクションにウィンドウを適用します。

  • 隣接するセクション間で 50% のオーバーラップを指定します。

  • FFT を計算するには、 $\max(256,2^p)$ 点を使用します。ここで、 $p=\lceil\log_2N_{\tt x}\rceil$ です。

2 つの方法の結果が同じになることを確認します。

Nx = length(x);
nsc = floor(Nx/4.5);
nov = floor(nsc/2);
nff = max(256,2^nextpow2(nsc));

t = spectrogram(x,hamming(nsc),nov,nff);

maxerr = max(abs(abs(t(:))-abs(s(:))))
maxerr =

     0

長さが等しく、セクション間のオーバーラップが 50% である 8 つのセクションに信号を分割します。前のステップと同じ FFT 長を指定します。短時間フーリエ変換を計算し、前の 2 つの手順と結果が同じであることを確認します。

ns = 8;
ov = 0.5;
lsc = floor(Nx/(ns-(ns-1)*ov));

t = spectrogram(x,lsc,floor(ov*lsc),nff);

maxerr = max(abs(abs(t(:))-abs(s(:))))
maxerr =

     0

1 kHz で 2 秒間サンプリングされた 2 次チャープ x を生成します。チャープの初期周波数は 100 Hz で、t = 1 秒で 200 Hz になります。

t = 0:0.001:2;
x = chirp(t,100,1,200,'quadratic');

x のスペクトログラムを計算して表示します。

  • 信号にハミング ウィンドウを適用し、長さ 128 のセクションに分割します。

  • 隣り合ったセクション間のオーバーラップのサンプルを 120 個に指定します。

  • 周波数 $\lfloor128/2+1\rfloor=65$ および時間ビン $\lfloor({\tt length(x)}-120)/(128-120)\rfloor=235$ でスペクトルを評価します。

spectrogram(x,128,120,128,1e3)

ハミング ウィンドウをブラックマン ウィンドウに置き換えます。オーバーラップを 60 サンプルに減らします。値が上から下へ向かって増加するように時間軸をプロットします。

spectrogram(x,blackman(128),60,128,1e3)
ax = gca;
ax.YDir = 'reverse';

100 Hz から開始し、t = 1 秒で 200 Hz になる 2 次チャープの各セグメントの PSD を計算して表示します。サンプルレートを 1 kHz、セグメントの長さを 128 サンプル、オーバーラップのサンプルを 120 個に指定します。128 個の DFT 点と既定のハミング ウィンドウを使用します。

t = 0:0.001:2;
x = chirp(t,100,1,200,'quadratic');

spectrogram(x,128,120,128,1e3,'yaxis')
title('Quadratic Chirp')

DC から開始し、t = 1 秒で 150 Hz になる線形チャープの各セグメントの PSD を計算して表示します。サンプルレートを 1 kHz、セグメントの長さを 256 サンプル、オーバーラップのサンプルを 250 個に指定します。既定のハミング ウィンドウと 256 個の DFT 点を使用します。

t = 0:0.001:2;
x = chirp(t,0,1,150);

spectrogram(x,256,250,256,1e3,'yaxis')
title('Linear Chirp')

20 Hz で開始し、t = 1 秒で 60 Hz になる対数チャープを 1 kHz でサンプリングし、その対数チャープの各セグメントの PSD を計算して表示します。セグメントの長さを 256 サンプル、オーバーラップのサンプルを 250 個に指定します。既定のハミング ウィンドウと 256 個の DFT 点を使用します。

t = 0:0.001:2;
x = chirp(t,20,1,60,'logarithmic');

spectrogram(x,256,250,[],1e3,'yaxis')
title('Logarithmic Chirp')

周波数軸に対数スケールを使用します。スペクトログラムは直線になります。

ax = gca;
ax.YScale = 'log';

関数 spectrogram を使用して信号の瞬間周波数を測定し、追跡します。

1 kHz で 2 秒間サンプリングされた 2 次チャープを生成します。初期周波数が 100 Hz で、1 秒後には 200 Hz に増加するようにチャープを指定します。

Fs = 1000;
t = 0:1/Fs:2-1/Fs;
y = chirp(t,100,1,200,'quadratic');

関数 spectrogram に実装された短時間フーリエ変換を使用して、チャープのスペクトルを推定します。信号にハミング ウィンドウを適用し、長さ 100 のセクションに分割します。隣り合ったセクション間のオーバーラップのサンプル 80 個を指定し、周波数 $\lfloor100/2+1\rfloor=51$ におけるスペクトルを評価します。既定のカラー バーを非表示にします。

spectrogram(y,100,80,100,Fs,'yaxis')
view(-77,72)
shading interp
colorbar off

$\lfloor(2000-80)/(100-80)\rfloor = 96$ の各時点でのパワー スペクトル密度の最大値を検出することにより、チャープの周波数を追跡します。スペクトログラムを 2 次元グラフィックスで表示します。カラー バーを元に戻します。

[s,f,t,p] = spectrogram(y,100,80,100,Fs);

[q,nd] = max(10*log10(p));

hold on
plot3(t,f(nd),q,'r','linewidth',4)
hold off
colorbar
view(2)

1 kHz で 2 秒間サンプリングされたチャープ信号を生成します。初期周波数が 100 Hz で、1 秒後には 200 Hz に増加するようにチャープを指定します。

Fs = 1000;
t = 0:1/Fs:2;
y = chirp(t,100,1,200,'quadratic');

再割り当てされた信号のスペクトログラムを推定します。

  • 信号に形状パラメーター $\beta=18$ のカイザー ウィンドウを適用し、長さ 128 のセクションに分割します。

  • 隣り合ったセクション間のオーバーラップのサンプルを 120 個に指定します。

  • 周波数 $\lfloor128/2\rfloor=65$ および時間ビン $\lfloor({\tt length(x)}-120)/(128-120)\rfloor=235$ でスペクトルを評価します。

spectrogram(y,kaiser(128,18),120,128,Fs,'reassigned','yaxis')

1 kHz で 2 秒間サンプリングされたチャープ信号を生成します。初期周波数が 100 Hz で、1 秒後には 200 Hz に増加するようにチャープを指定します。

Fs = 1000;
t = 0:1/Fs:2;
y = chirp(t,100,1,200,'quadratic');

信号の時間依存のパワー スペクトル密度 (PSD) を推定します。

  • 信号に形状パラメーター $\beta=18$ のカイザー ウィンドウを適用し、長さ 128 のセクションに分割します。

  • 隣り合ったセクション間のオーバーラップのサンプルを 120 個に指定します。

  • 周波数 $\lfloor128/2\rfloor=65$ および時間ビン $\lfloor({\tt length(x)}-120)/(128-120)\rfloor=235$ でスペクトルを評価します。

各 PSD 推定の重心の周波数および時間を出力します。 $-30$ dB より小さい PSD の要素をゼロに設定します。

[~,~,~,pxx,fc,tc] = spectrogram(y,kaiser(128,18),120,128,Fs, ...
    'MinThreshold',-30);

重心の周波数および時間の関数として非ゼロの要素をプロットします。

plot(tc(pxx>0),fc(pxx>0),'.')

1024 Hz で 2 秒間サンプリングされた信号を生成します。

nSamp = 2048;
Fs = 1024;
t = (0:nSamp-1)'/Fs;

最初の 1 秒間、信号は 400 Hz の正弦波と凹の 2 次チャープで構成されています。開始および終了時の周波数が 250 Hz で、最小でも 150 Hz に到達するチャープを区間の中間点で対称となるように指定します。

t1 = t(1:nSamp/2);

x11 = sin(2*pi*400*t1);
x12 = chirp(t1-t1(nSamp/4),150,nSamp/Fs,1750,'quadratic');
x1 = x11+x12;

信号の他の部分は、減少する周波数の 2 つの線形チャープで構成されています。一方のチャープの初期周波数は 250 Hz で、100 Hz まで減少します。もう一方のチャープの初期周波数は 400 Hz で、250 Hz まで減少します。

t2 = t(nSamp/2+1:nSamp);

x21 = chirp(t2,400,nSamp/Fs,100);
x22 = chirp(t2,550,nSamp/Fs,250);
x2 = x21+x22;

ホワイト ガウス ノイズを信号に付加します。20 dB の S/N 比を指定します。再現可能な結果が必要な場合は、乱数発生器をリセットします。

SNR = 20;
rng('default')

sig = [x1;x2];
sig = sig + randn(size(sig))*std(sig)/db2mag(SNR);

信号のスペクトログラムを計算してプロットします。カイザー ウィンドウを長さ 63、形状パラメーター $\beta=17$ に、隣接するセグメント間のオーバーラップのサンプルを 10 個少ない数に、FFT 長を 256 にそれぞれ指定します。

nwin = 63;
wind = kaiser(nwin,17);
nlap = nwin-10;
nfft = 256;

spectrogram(sig,wind,nlap,nfft,Fs,'yaxis')

SNR より小さい値をもつ要素がゼロに設定されるように、スペクトログラムをしきい値処理します。

spectrogram(sig,wind,nlap,nfft,Fs,'MinThreshold',-SNR,'yaxis')

各 PSD 推定を、そのエネルギー中心の位置に再代入します。

spectrogram(sig,wind,nlap,nfft,Fs,'reassign','yaxis')

SNR より小さい値をもつ要素がゼロに設定されるように、再割り当てされたスペクトログラムをしきい値処理します。

spectrogram(sig,wind,nlap,nfft,Fs,'reassign','MinThreshold',-SNR,'yaxis')

減少する 2 つのチャープと広帯域のスプラッター音を含むオーディオ信号を読み込みます。短時間フーリエ変換を計算します。波形を 300 サンプルがオーバーラップする 400 サンプルのセグメントに分割します。スペクトログラムをプロットします。

load splat

% To hear, type soundsc(y,Fs)

sg = 400;
ov = 300;

spectrogram(y,sg,ov,[],Fs,'yaxis')
colormap bone

関数 spectrogram を使用し、信号のパワー スペクトル密度 (PSD) を出力します。

[s,f,t,p] = spectrogram(y,sg,ov,[],Fs);

関数 medfreq を使用して 2 つのチャープを追跡します。より強い低周波のチャープを検出するには、探索を 100 Hz より高い周波数と広帯域音が始まる前の時間に制限します。

f1 = f > 100;
t1 = t < 0.75;

m1 = medfreq(p(f1,t1),f(f1));

微弱な高周波のチャープを検出するには、探索を 2500 Hz より高い周波数と 0.3 秒 ~ 0.65 秒の時間に制限します。

f2 = f > 2500;
t2 = t > 0.3 & t < 0.65;

m2 = medfreq(p(f2,t2),f(f2));

結果をスペクトログラムに重ね合わせます。kHz の単位で表すために周波数の値を 1000 で除算します。

hold on
plot(t(t1),m1/1000,'linewidth',4)
plot(t(t2),m2/1000,'linewidth',4)
hold off

10 kHz でサンプリングされた 2 秒間の信号を生成します。信号の瞬間周波数を時間の三角形関数として指定します。

fs = 10e3;
t = 0:1/fs:2;
x1 = vco(sawtooth(2*pi*t,0.5),[0.1 0.4]*fs,fs);

信号のスペクトログラムを計算してプロットします。長さ 256 で形状パラメーター $\beta=5$ のカイザー ウィンドウを使用します。セクション間のオーバーラップのサンプル 220 個と DFT 点 512 個を指定します。周波数を y 軸にプロットします。既定のカラーマップとビューを使用します。

spectrogram(x1,kaiser(256,5),220,512,fs,'yaxis')

スペクトログラムをウォーターフォール プロットとして表示するようにビューを変更します。カラーマップを bone に設定します。

colormap bone
view(-45,65)

関連する例

入力引数

すべて折りたたむ

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

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

データ型: single | double
複素数のサポート: はい

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

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

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

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

window を空として指定した場合、spectrogram は、xnoverlap 個のオーバーラップ サンプルをもつ 8 個のセクションに分割されているハミング ウィンドウを使用します。

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

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

データ型: single | double

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

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

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

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

データ型: double | single

DFT 点の数。正の整数スカラーとして指定します。nfft を空として指定した場合、spectrogram によりパラメーターが max(256,2p) に設定されます。ここで p = ⌈log2 Nx⌉ は入力信号の長さ Nx についての値です。

データ型: single | double

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

例: pi./[2 4]

データ型: double | single

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

データ型: double | single

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

データ型: double | single

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

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

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

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

データ型: char

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

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

  • 'power' を指定すると、ウィンドウの等価ノイズ帯域幅ごとに PSD 推定をスケーリングします。結果は、各周波数のパワーの推定です。'reassigned' オプションがオンの場合、関数は、再代入の前に PSD を各周波数ビンの幅で積分します。

データ型: char

しきい値。dB 単位で表される実数スカラーで指定します。spectrogram は、10 log10(ps) ≤ thresh となる ps の要素をゼロに設定します。

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

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

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

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

データ型: char

出力引数

すべて折りたたむ

短時間フーリエ変換。行列として返されます。時間は s の列方向に、周波数は行方向に下がって 0 から増加します。

  • x が長さ Nx の信号の場合、s は k 列になります。ここで

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

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

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

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

  • x が複素数の場合、snfft 行になります。

データ型: double | single

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

データ型: double | single

時点。ベクトルとして返されます。t の時間値は、各セクションの中間点に対応します。

データ型: double | single

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

データ型: double | single

パワー スペクトル密度 (PSD) またはパワー スペクトル。行列として返されます。

  • x が実数の場合、ps は各セクションの PSD またはパワー スペクトルの片側修正ピリオドグラム推定を含みます。

  • x が複素数の場合または周波数のベクトルを指定した場合には、ps は各セクションの PSD またはパワー スペクトルの両側修正ピリオドグラム推定を含みます。

データ型: double | single

エネルギー中心の周波数および時間。短時間フーリエ変換と同じサイズの行列として返されます。サンプルレートを指定しない場合、fc の要素は正規化周波数として返されます。

詳細

すべて折りたたむ

ヒント

短時間フーリエ変換にゼロがある場合、デシベルに変換すると、プロットできない負の無限大になります。この問題の発生を避けるため、spectrogram を出力引数なしで呼び出した場合には短時間フーリエ変換に eps が加えられます。

参照

[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 Ronald W. Schafer. Digital Processing of Speech Signals. Englewood Cliffs, NJ: Prentice-Hall, 1978.

[3] Chassande-Motin, Éric, François Auger, and Patrick Flandrin. “Reassignment.” In Time-Frequency Analysis: Concepts and Methods. Edited by Franz Hlawatsch and François Auger. London: ISTE/John Wiley and Sons, 2008.

[4] Fulop, Sean A., and Kelly Fitz. “Algorithms for computing the time-corrected instantaneous frequency (reassigned) spectrogram, with applications.” Journal of the Acoustical Society of America. Vol. 119, January 2006, pp. 360–371.

R2006a より前に導入

この情報は役に立ちましたか?