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(___,spectrumtype)

[___] = spectrogram(___,"reassigned")

[___,ps,fc,tc] = spectrogram(___)

[___] = spectrogram(___,freqrange)

[___] = spectrogram(___,Name=Value)

spectrogram(___)

spectrogram(___,freqloc)

説明

s = spectrogram(x) は、入力信号 x の短時間フーリエ変換 (STFT) を返します。s の各列は、x の短時間の、時間が局所化された周波数成分の推定を含みます。s の振幅の 2 乗は、x の "スペクトログラム" 時間-周波数表現と呼ばれます[1]。

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

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

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

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

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

[s,w,t] = spectrogram(x,window,noverlap,w) は、w で指定した正規化周波数における STFT を返します。w は少なくとも 2 つの要素をもたなければなりません。そうでない場合は、関数が nfft として解釈するためです。

[s,f,t] = spectrogram(x,window,noverlap,f,fs) は、f で指定した巡回周波数における STFT を返します。f は少なくとも 2 つの要素をもたなければなりません。そうでない場合は、関数が nfft として解釈するためです。

また、[___,ps] = spectrogram(___,spectrumtype) は x のスペクトログラムに比例する行列 ps も返します。

• spectrumtype を "psd" として指定すると、ps の各列には、ウィンドウが適用されたセグメントのパワー スペクトル密度 (PSD) の推定が格納されます。

• spectrumtype を "power" として指定すると、ps の各列には、ウィンドウが適用されたセグメントのパワー スペクトルの推定が格納されます。

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

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

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

[___] = spectrogram(___,Name=Value) は、名前と値の引数を使用して追加オプションを指定します。オプションには、最小しきい値と出力時間次元が含まれます。

出力引数なしで spectrogram(___) を使用すると、現在の Figure ウィンドウに ps がデシベル単位でプロットされます。

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

例

スペクトログラムの既定値

正弦波の和から構成される信号のサンプルを $$N_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')

計算を繰り返します。

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

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

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

• FFT を計算するには、$$\max (256,2^p)$$ 点を使用します。ここで、$$p=\lceil {\log }_{2}nsc\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

関数 spectrogram と STFT 定義の比較

600 Hz で 2 秒間サンプリングされた複素数値凸 2 次チャープで構成される信号を生成します。チャープの初期周波数は 250 Hz、最終周波数は 50 Hz です。

fs = 6e2;
ts = 0:1/fs:2;
x = chirp(ts,250,ts(end),50,"quadratic",0,"convex","complex");

関数 spectrogram

関数 spectrogram を使用して信号の STFT を計算します。

• 信号を、それぞれ $\mathit{M}=49$ サンプルの長さのセグメントに分割します。

• 隣接するセグメント間に $\mathit{L}=11$ 個のサンプルのオーバーラップを指定します。

• 最後の短いセグメントを破棄します。

• 各セグメントにバートレット ウィンドウを適用します。

• ${\mathit{N}}_{\mathrm{DFT}}=1024$ 点で各セグメントの離散フーリエ変換を評価します。既定では、複素数値信号の場合、spectrogram は両側変換を計算します。

M = 49;
L = 11;
g = bartlett(M);
Ndft = 1024;

[s,f,t] = spectrogram(x,g,L,Ndft,fs);

関数 waterplot を使用して、STFT の振幅の 2 乗として定義されるスペクトログラムを計算し、表示します。

waterplot(s,f,t)

STFT の定義

定義を使用して ${\mathit{N}}_{\mathit{x}}$ サンプルの信号の STFT を計算します。信号を $$\lfloor \frac{N_x-L}{M-L}\rfloor$$ 個のオーバーラップ セグメントに分割します。各セグメントにウィンドウを適用し、それぞれについて、${\mathit{N}}_{\mathrm{DFT}}$ 点での離散フーリエ変換を評価します。

[segs,~] = buffer(1:length(x),M,L,"nodelay");

X = fft(x(segs).*g,Ndft);

STFT の時間範囲と周波数範囲を計算します。

• 時間値を求めるには、時間ベクトルをオーバーラップ セグメントに分割します。時間値は、下端でオープンの区間として扱われる各セグメントの中間点です。

• 周波数値を求めるには、ゼロ周波数でクローズ、下端でオープンとなるナイキスト区間を指定します。

tbuf = ts(segs);
tint = mean(tbuf(2:end,:));

fint = 0:fs/Ndft:fs-fs/Ndft;

spectrogram の出力を定義と比較します。スペクトログラムを表示します。

maxdiff = max(max(abs(s-X)))

maxdiff = 0

waterplot(X,fint,tint)

function waterplot(s,f,t)
% Waterfall plot of spectrogram
    waterfall(f,t,abs(s)'.^2)
    set(gca,XDir="reverse",View=[30 50])
    xlabel("Frequency (Hz)")
    ylabel("Time (s)")
end

関数 spectrogram と stft の比較

1.4 kHz で 2 秒間サンプリングされたチャープで構成される信号を生成します。チャープの周波数は、測定時間中に 600 Hz から 100 Hz に線形的に減少します。

fs = 1400;
x = chirp(0:1/fs:2,600,2,100);

stft 既定

関数 spectrogram と stft を使用して信号の STFT を計算します。関数 stft の既定の値を使用します。

• 信号を 128 サンプルのセグメントに分割し、各セグメントに周期的ハン ウィンドウを適用します。

• 隣接するセグメント間に 96 個のサンプルのオーバーラップを指定します。この長さはウィンドウの長さの 75% と等価です。

• 128 の DFT 点を指定し、STFT の中央を周波数ゼロに揃え、周波数は Hz 単位で表します。

2 つの結果が等価であることを確認します。

M = 128;
g = hann(M,"periodic");
L = 0.75*M;
Ndft = 128;

[sp,fp,tp] = spectrogram(x,g,L,Ndft,fs,"centered");

[s,f,t] = stft(x,fs);

dff = max(max(abs(sp-s)))

dff = 0

関数 mesh を使用して、2 つの出力をプロットします。

nexttile
mesh(tp,fp,abs(sp).^2)
title("spectrogram")
view(2), axis tight

nexttile
mesh(t,f,abs(s).^2)
title("stft")
view(2), axis tight

spectrogram 既定

関数 spectrogram の既定の値を使用して計算を繰り返します。

• 信号を $\mathit{M}=\lfloor {\mathit{N}}_{\mathit{x}}/4.5\rfloor$ の長さのセグメントに分割します。ここで、${\mathit{N}}_{\mathit{x}}$ は信号の長さです。各セグメントにハミング ウィンドウを適用します。

• セグメント間で 50% のオーバーラップを指定します。

• FFT を計算するには、$\max \left(256,2^{\lceil {\log }_2\mathit{M}\rceil }\right)$ 点を使用します。正の正規化周波数のみ、スペクトログラムを計算します。

M = floor(length(x)/4.5);
g = hamming(M);
L = floor(M/2);
Ndft = max(256,2^nextpow2(M));

[sx,fx,tx] = spectrogram(x);

[st,ft,tt] = stft(x,Window=g,OverlapLength=L, ...
    FFTLength=Ndft,FrequencyRange="onesided");

dff = max(max(sx-st))

dff = 0

関数 waterplot を使用して、2 つの出力をプロットします。どちらの場合も、周波数軸を $\pi$ で除算します。stft 出力の場合、サンプル数を有効なサンプル レート $2\pi$ で除算します。

figure
nexttile
waterplot(sx,fx/pi,tx)
title("spectrogram")

nexttile
waterplot(st,ft/pi,tt/(2*pi))
title("stft")

function waterplot(s,f,t)
% Waterfall plot of spectrogram
    waterfall(f,t,abs(s)'.^2)
    set(gca,XDir="reverse",View=[30 50])
    xlabel("Frequency/\pi")
    ylabel("Samples")
end

スペクトログラムと瞬時周波数

関数 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 個を指定し、周波数 $$\lfloor 100/2+1\rfloor =51$$ におけるスペクトルを評価します。

spectrogram(y,100,80,100,fs,'yaxis')

$$\lfloor (2000-80)/(100-80)\rfloor =96$$ の時間点にわたる最大エネルギーをもつ時間-周波数リッジを検出することにより、チャープの周波数を追跡します。スペクトログラム プロットに瞬時周波数を重ね合わせます。

[~,f,t,p] = spectrogram(y,100,80,100,fs);

[fridge,~,lr] = tfridge(p,f);

hold on
plot3(t,fridge,abs(p(lr)),'LineWidth',4)
hold off

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

正弦関数的に変化する周波数成分をもつチャープのサンプルを 512 個生成します。

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

x = exp(1j*pi*sin(8*n/N)*32);

チャープの中央揃えの短時間フーリエ変換を計算します。信号を 16 サンプルがオーバーラップする 32 サンプルのセグメントに分割します。64 個の DFT ポイントを指定します。スペクトログラムをプロットします。

[scalar,fs,ts] = spectrogram(x,32,16,64,'centered');

spectrogram(x,32,16,64,'centered','yaxis')

区間 $(-\pi ,\pi ]$ に対して 64 等間隔周波数のスペクトログラムを計算すると、同じ結果が得られます。'centered' オプションは必要ありません。

fintv = -pi+pi/32:pi/32:pi;

[vector,fv,tv] = spectrogram(x,32,16,fintv);

spectrogram(x,32,16,fintv,'yaxis')

関数 spectrogram と pspectrum の比較

電圧制御発振器と 3 つのガウス原子で構成される信号を生成します。信号は ${\mathit{f}}_{\mathit{s}}=2$ kHz で 2 秒間サンプリングされます。

fs = 2000;
tx = 0:1/fs:2;
gaussFun = @(A,x,mu,f) exp(-(x-mu).^2/(2*0.03^2)).*sin(2*pi*f.*x)*A';
s = gaussFun([1 1 1],tx',[0.1 0.65 1],[2 6 2]*100)*1.5;
x = vco(chirp(tx+.1,0,tx(end),3).*exp(-2*(tx-1).^2),[0.1 0.4]*fs,fs);
x = s+x';

短時間フーリエ変換

関数 pspectrum を使用して STFT を計算します。

• ${\mathit{N}}_{\mathit{x}}$ サンプル信号を、$80/2000=40$ ミリ秒の時間分解能に対応する長さ $\mathit{M}=80$ のサンプルのセグメントに分割します。

• $\mathit{L}=16$ 個のサンプル、または隣接するセグメント間の 20% のオーバーラップを指定します。

• カイザー ウィンドウを使用して各セグメントにウィンドウを適用し、漏れには $\ell =0.7$ を指定します。

M = 80;
L = 16;
lk = 0.7;

[S,F,T] = pspectrum(x,fs,"spectrogram", ...
    TimeResolution=M/fs,OverlapPercent=L/M*100, ...
    Leakage=lk);

関数 spectrogram で得られた結果と比較します。

• ウィンドウの長さとオーバーラップをサンプル単位で直接指定します。

• pspectrum は常にカイザー ウィンドウを $\mathit{g}\left(\mathit{n}\right)$ として使用します。漏れ $\ell$ とウィンドウの形状係数 $\beta$ は、$\beta =40\times \left(1-\ell \right)$ の関係にあります。

• pspectrum は、離散フーリエ変換の計算に常に ${\mathit{N}}_{\mathrm{DFT}}=1024$ 個の点を使用します。両側または中央揃えの周波数範囲の変換を計算したい場合にこの数値を指定することができます。一方、実信号の既定値である片側変換の場合、spectrogram は $1024/2+1=513$ 個の点を使用します。あるいは、この例にあるように、変換を計算したい周波数のベクトルを指定することができます。

• 信号を $\mathit{k}=\lfloor \frac{{\mathit{N}}_{\mathit{x}}-\mathit{L}}{\mathit{M}-\mathit{L}}\rfloor$ 個のセグメントに厳密に分割できない場合、spectrogram は信号を切り捨て、pspectrum は信号をゼロでパディングして追加のセグメントを作成します。出力を等価にするには、最後のセグメントと時間ベクトルの最後の要素を削除します。

• spectrogram は、振幅の 2 乗をスペクトログラムとする STFT を返します。pspectrum は、あらかじめ係数 ${\sum }_{\mathit{n}}\mathit{g}\left(\mathit{n}\right)$ で除算してから 2 乗した各セグメントのパワー スペクトルを返します。

• 片側変換の場合、pspectrum はスペクトログラムに追加の係数 2 を追加します。

g = kaiser(M,40*(1-lk));

k = (length(x)-L)/(M-L);
if k~=floor(k)
    S = S(:,1:floor(k));
    T = T(1:floor(k));
end

[s,f,t] = spectrogram(x/sum(g)*sqrt(2),g,L,F,fs);

関数 waterplot を使用して、これら 2 つの関数によって計算されたスペクトログラムを表示します。

subplot(2,1,1)
waterplot(sqrt(S),F,T)
title("pspectrum")

subplot(2,1,2)
waterplot(s,f,t)
title("spectrogram")

maxd = max(max(abs(abs(s).^2-S)))

maxd = 2.4419e-08

パワー スペクトルと簡易プロット

関数 spectrogram は、各セグメントのパワー スペクトルまたはパワー スペクトル密度に対応する 4 番目の引数をもちます。pspectrum の出力と同様、ps 引数はあらかじめ 2 乗されており、正規化係数 ${\sum }_{\mathit{n}}\mathit{g}\left(\mathit{n}\right)$ を含みます。実信号の片側スペクトログラムの場合も、追加の係数 2 を含める必要があります。関数のスケーリング引数を "power" に設定します。

[~,~,~,ps] = spectrogram(x*sqrt(2),g,L,F,fs,"power");

max(abs(S(:)-ps(:)))

ans = 2.4419e-08

pspectrum と spectrogram は、出力引数なしで呼び出された場合に、信号のスペクトログラムをデシベル単位でプロットします。片側スペクトログラムの場合は係数 2 を含めます。両方のプロットで同じになるようにカラーマップを設定します。x の範囲を同じ値に設定して、pspectrum プロットの最後にある追加のセグメントが表示されるようにします。spectrogram プロットの y 軸に周波数を表示します。

subplot(2,1,1)
pspectrum(x,fs,"spectrogram", ...
    TimeResolution=M/fs,OverlapPercent=L/M*100, ...
    Leakage=lk)
title("pspectrum")
cc = clim;
xl = xlim;

subplot(2,1,2)
spectrogram(x*sqrt(2),g,L,F,fs,"power","yaxis")
title("spectrogram")
clim(cc)
xlim(xl)

function waterplot(s,f,t)
% Waterfall plot of spectrogram
    waterfall(f,t,abs(s)'.^2)
    set(gca,XDir="reverse",View=[30 50])
    xlabel("Frequency (Hz)")
    ylabel("Time (s)")
end

再割り当てされた 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 個に指定します。

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

出力引数なしで関数 spectrogram を使用して、再割り当てされたスペクトログラムをプロットします。y 軸に周波数、x 軸に時間を表示します。

spectrogram(y,kaiser(128,18),120,128,fs, ...
    "reassigned","yaxis")

関数 imagesc を使用してプロットをやり直します。周波数値が下から上に増加するように、y 軸の方向を指定します。デシベル変換時に負の無限大となることを避けるため、再割り当てされたスペクトログラムに eps を追加します。

[~,fr,tr,pxx] = spectrogram(y,kaiser(128,18),120,128,fs, ...
    "reassigned");

imagesc(tr,fr,pow2db(pxx+eps))
axis xy
xlabel("Time (s)")
ylabel("Frequency (Hz)")
colorbar

しきい値を使用したスペクトログラム

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 個に指定します。

• $$\lfloor 128/2\rfloor =65$$ 周波数および $$\lfloor (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),'.')

中央揃えのスペクトログラムと片側スペクトログラムの計算

2 kHz で 2 秒間サンプリングされた実数値のチャープで構成される信号を生成します。

fs = 2000;
tx = 0:1/fs:2;
x = vco(-chirp(tx,0,tx(end),2).*exp(-3*(tx-1).^2), ...
    [0.1 0.4]*fs,fs).*hann(length(tx))';

両側スペクトログラム

信号の両側 STFT を計算してプロットします。

• 信号を、それぞれ $\mathit{M}=73$ サンプルの長さのセグメントに分割します。

• 隣接するセグメント間に $\mathit{L}=24$ 個のサンプルのオーバーラップを指定します。

• 最後の短いセグメントを破棄します。

• 各セグメントにフラットトップ ウィンドウを適用します。

• ${\mathit{N}}_{\mathrm{DFT}}=895$ 点で各セグメントの離散フーリエ変換を評価します。このとき、点の数が奇数であることに注意します。

M = 73;
L = 24;
g = flattopwin(M);
Ndft = 895;
neven = ~mod(Ndft,2);

[stwo,f,t] = spectrogram(x,g,L,Ndft,fs,"twosided");

出力引数なしで関数 spectrogram を使用し、両側スペクトログラムをプロットします。

spectrogram(x,g,L,Ndft,fs,"twosided","power","yaxis")

定義を使用して両側スペクトログラムを計算します。隣接するセグメント間で $\mathit{L}$ サンプルがオーバーラップする $\mathit{M}$ サンプルのセグメントに信号を分割します。各セグメントにウィンドウを適用し、それぞれについて、${\mathit{N}}_{\mathrm{DFT}}$ 点での離散フーリエ変換を計算します。

[segs,~] = buffer(1:length(x),M,L,"nodelay");

Xtwo = fft(x(segs).*g,Ndft);

時間範囲と周波数範囲を計算します。

• 時間値を求めるには、時間ベクトルをオーバーラップ セグメントに分割します。時間値は、下端でオープンの区間として扱われる各セグメントの中間点です。

• 周波数値を求めるには、ゼロ周波数でクローズ、上端でオープンとなるナイキスト区間を指定します。

tbuf = tx(segs);
ttwo = mean(tbuf(2:end,:));

ftwo = 0:fs/Ndft:fs*(1-1/Ndft);

spectrogram の出力を定義と比較します。関数 waterplot を使用して、スペクトログラムを表示します。

diffs = [max(max(abs(stwo-Xtwo))) max(abs(f-ftwo')) max(abs(t-ttwo))]

diffs = 1×3
10-12 ×

         0    0.2274    0.0002

figure
nexttile
waterplot(Xtwo,ftwo,ttwo)
title("Two-Sided, Definition")

nexttile
waterplot(stwo,f,t)
title("Two-Sided, spectrogram Function")

中央揃えのスペクトログラム

信号の中央揃えのスペクトログラムを計算します。

• 両側 STFT で使用したのと同じ時間値を使用します。

• 関数 fftshift を使用して、STFT のゼロ周波数成分をスペクトルの中央にシフトします。

• ${\mathit{N}}_{\mathrm{DFT}}$ が奇数の場合、周波数範囲は両側でオープンとなります。${\mathit{N}}_{\mathrm{DFT}}$ が偶数の場合、周波数範囲は下端でオープン、上端でクローズとなります。

出力を比較してスペクトログラムを表示します。

tcen = ttwo;

if ~neven
    Xcen = fftshift(Xtwo,1);
    fcen = -fs/2*(1-1/Ndft):fs/Ndft:fs/2;
else
    Xcen = fftshift(circshift(Xtwo,-1),1);
    fcen = (-fs/2*(1-1/Ndft):fs/Ndft:fs/2)+fs/Ndft/2;
end

[scen,f,t] = spectrogram(x,g,L,Ndft,fs,"centered");

diffs = [max(max(abs(scen-Xcen))) max(abs(f-fcen')) max(abs(t-tcen))]

diffs = 1×3
10-12 ×

         0    0.2274    0.0002

figure
nexttile
waterplot(Xcen,fcen,tcen)
title("Centered, Definition")

nexttile
waterplot(scen,f,t)
title("Centered, spectrogram Function")

片側スペクトログラム

信号の片側スペクトログラムを計算します。

• 両側 STFT で使用したのと同じ時間値を使用します。

• ${\mathit{N}}_{\mathrm{DFT}}$ が奇数の場合、片側 STFT は両側 STFT の最初の $\left({\mathit{N}}_{\mathrm{DFT}}+1\right)/2$ 行により構成されます。${\mathit{N}}_{\mathrm{DFT}}$ が偶数の場合、片側 STFT は両側 STFT の最初の ${\mathit{N}}_{\mathrm{DFT}}/2+1$ 行により構成されます。

• ${\mathit{N}}_{\mathrm{DFT}}$ が奇数の場合、周波数範囲はゼロ周波数でクローズ、ナイキスト周波数でオープンとなります。${\mathit{N}}_{\mathrm{DFT}}$ が偶数の場合、周波数範囲は両側でクローズとなります。

出力を比較してスペクトログラムを表示します。実数値信号の場合、"onesided" 引数はオプションです。

tone = ttwo;

if ~neven
    Xone = Xtwo(1:(Ndft+1)/2,:);
else
    Xone = Xtwo(1:Ndft/2+1,:);
end

fone = 0:fs/Ndft:fs/2;

[sone,f,t] = spectrogram(x,g,L,Ndft,fs);

diffs = [max(max(abs(sone-Xone))) max(abs(f-fone')) max(abs(t-tone))]

diffs = 1×3
10-12 ×

         0    0.1137    0.0002

figure
nexttile
waterplot(Xone,fone,tone)
title("One-Sided, Definition")

nexttile
waterplot(sone,f,t)
title("One-Sided, spectrogram Function")

function waterplot(s,f,t)
% Waterfall plot of spectrogram
waterfall(f,t,abs(s)'.^2)
set(gca,XDir="reverse",View=[30 50])
xlabel("Frequency (Hz)")
ylabel("Time (s)")
end

セグメントの PSD とパワー スペクトルの計算

関数 spectrogram は、4 番目の出力引数として、各セグメントのパワー スペクトル密度 (PSD) またはパワー スペクトルのいずれかを含む行列をもちます。パワー スペクトルは、PSD にウィンドウの等価ノイズ帯域幅 (ENBW) を乗算した値と等価です。

1 kHz で 1 秒間サンプリングされた対数チャープで構成される信号を生成します。チャープの初期周波数は 400 Hz で、測定の最後には 10 Hz に減少します。

fs = 1000;
tt = 0:1/fs:1-1/fs;
y = chirp(tt,400,tt(end),10,"logarithmic");

サンプル レートによる PSD とパワー スペクトルのセグメント化

信号を 102 サンプルのセグメントに分割し、各セグメントにハン ウィンドウを適用します。隣り合ったセグメント間のオーバーラップを 12 サンプル、DFT 点を 1024 に指定します。

M = 102;
g = hann(M);
L = 12;
Ndft = 1024;

既定の PSD スペクトル タイプで信号のスペクトログラムを計算します。STFT、およびセグメントのパワー スペクトル密度の配列を出力します。

[s,f,t,p] = spectrogram(y,g,L,Ndft,fs);

スペクトル タイプを "power" に指定して計算を繰り返します。STFT、およびセグメントのパワー スペクトルの配列を出力します。

[r,~,~,q] = spectrogram(y,g,L,Ndft,fs,"power");

どちらの場合もスペクトログラムが同じになることを確認します。周波数に対数スケールを使用して、スペクトログラムをプロットします。

max(max(abs(s).^2-abs(r).^2))

ans = 0

waterfall(f,t,abs(s)'.^2)
set(gca,XScale="log",...
    XDir="reverse",View=[30 50])

パワー スペクトルが、パワー スペクトル密度にウィンドウの ENBW を乗算した値と等価であることを確認します。

max(max(abs(q-p*enbw(g,fs))))

ans = 1.1102e-16

セグメントのパワー スペクトルの行列が、スペクトログラムに比例していることを確認します。比例係数は、ウィンドウ要素の和の 2 乗です。

max(max(abs(s).^2-q*sum(g)^2))

ans = 3.4694e-18

正規化周波数による PSD とパワー スペクトルのセグメント化

計算を繰り返しますが、今度は正規化周波数で計算します。サンプル レートを $2\pi$ として指定した場合、結果は同じになります。

[~,~,~,pn] = spectrogram(y,g,L,Ndft);
[~,~,~,qn] = spectrogram(y,g,L,Ndft,"power");

max(max(abs(qn-pn*enbw(g,2*pi))))

ans = 1.1102e-16

オーディオ信号のチャープの追跡

減少する 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

3D スペクトログラムの可視化

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 に設定します。

view(-45,65)
colormap bone

入力引数

x — 入力信号

ベクトル

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

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

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

window — ウィンドウ

正の整数 | ベクトル | []

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

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

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

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

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

利用可能なウィンドウのリストについては、ウィンドウを参照してください。

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

noverlap — オーバーラップするサンプルの数

非負の整数 | []

オーバーラップするサンプルの数。非負の整数として指定します。

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

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

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

nfft — DFT 点の数

正の整数 | []

DFT 点の数。正の整数スカラーとして指定します。nfft を空として指定した場合、spectrogram によりパラメーターが max(256,2^p) に設定されます。ここで p = ⌈log₂ N_w⌉ で、⌈⌉ 記号は天井関数を表し、

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

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

w — 正規化周波数

ベクトル

正規化周波数。ベクトルとして指定します。w は少なくとも 2 つの要素をもたなければなりません。そうでない場合は、関数が nfft として解釈するためです。正規化周波数の単位はラジアン/サンプルです。

例: pi./[2 4]

データ型: single | double

f — 巡回周波数

ベクトル

巡回周波数。ベクトルとして指定します。f は少なくとも 2 つの要素をもたなければなりません。そうでない場合は、関数が nfft として解釈するためです。f の単位はサンプル レート fs により指定されます。

データ型: single | double

fs — サンプル レート

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

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

freqrange — PSD 推定の周波数範囲

"onesided" | "twosided" | "centered"

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 | string

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

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

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

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

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

データ型: char | string

freqloc — 周波数の表示軸

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

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

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

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

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

データ型: char | string

名前と値の引数

オプションの引数のペアを Name1=Value1,...,NameN=ValueN として指定します。ここで、Name は引数名で、Value は対応する値です。名前と値の引数は他の引数の後に指定しなければなりませんが、ペアの順序は重要ではありません。

例: spectrogram(x,100,OutputTimeDimension="downrows") は x を長さ 100 のセグメントに分割し、各セグメントにその長さのハミング ウィンドウを適用します。スペクトログラムの出力には、行方向に沿った時間次元があります。

R2021a より前では、コンマを使用して名前と値をそれぞれ区切り、Name を引用符で囲みます。

例: spectrogram(x,100,'OutputTimeDimension','downrows') は x を長さ 100 のセグメントに分割し、各セグメントにその長さのハミング ウィンドウを適用します。スペクトログラムの出力には、行方向に沿った時間次元があります。

MinThreshold — しきい値

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

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

OutputTimeDimension — 出力時間次元

"acrosscolumns" (既定値) | "downrows"

出力時間次元。"acrosscolumns" または "downrows" として指定します。行に沿った s、ps、fc、および tc の時間次元と、列に沿った周波数次元が必要な場合は、この値を "downrows" に設定します。列に沿った s、ps、fc、および tc の時間次元と、行に沿った周波数次元が必要な場合は、この値を "acrosscolumns" に設定します。関数が出力引数なしで呼び出された場合、この入力は無視されます。

データ型: char | string

出力引数

s — 短時間フーリエ変換

行列

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

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

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

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

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

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

• x が複素数値の場合、s は nfft 行になります。

メモ

freqrange が "onesided" に設定されている場合、spectrogram は、正のナイキスト範囲の s 値を合計パワーを保存せずに出力します。

s は "reassigned" オプションの影響は受けません。

w — 正規化周波数

ベクトル

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

t — 時点

ベクトル

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

f — 巡回周波数

ベクトル

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

ps — パワー スペクトル密度またはパワー スペクトル

行列

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

• x が実数で freqrange が指定されない場合、または "onesided" に設定されている場合、ps は各セグメントの PSD またはパワー スペクトルの片側修正ピリオドグラム推定を含みます。この関数は 0 とナイキスト周波数以外のすべての周波数でパワーを 2 倍にして、合計パワーを保存します。

• x が複素数値の場合、もしくは freqrange が "twosided" または "centered" に設定されている場合、ps は各セグメントの PSD またはパワー スペクトルの両側修正ピリオドグラム推定を含みます。

• w に正規化周波数のベクトルを指定した場合、または f に巡回周波数のベクトルを指定した場合、ps は入力周波数で評価された各セグメントの PSD またはパワー スペクトルの修正ピリオドグラム推定を含みます。

fc, tc — エネルギー中心の周波数および時間

行列

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

詳細

短時間フーリエ変換

短時間フーリエ変換 (STFT) を使用して、非定常信号の周波数成分が時間の経過と共に変化する様子を解析します。STFT の振幅の 2 乗は、信号の "スペクトログラム" 時間-周波数表現と呼ばれます。スペクトログラムの詳細や、関数 Signal Processing Toolbox™ を使用したスペクトログラムの計算方法については、Signal Processing Toolbox を使用したスペクトログラムの計算を参照してください。

信号の STFT は、信号上の長さ M の "解析ウィンドウ" g(n) をスライドして、ウィンドウが適用されたデータの各セグメントの離散フーリエ変換 (DFT) を計算することによって算出されます。ウィンドウは、R サンプルの間隔で元の信号を飛び越えます。これは、隣り合ったセグメント間の L = M – R 個のサンプルのオーバーラップに相当します。ほとんどのウィンドウ関数は、スペクトル リンギングを回避するためにエッジで小さくなります。ウィンドウが適用された各セグメントの DFT は、時間と周波数の各点の振幅と位相を含む複素数値行列に対して追加されます。STFT 行列は次の列数をもちます。

$$k=\lfloor \frac{N_x-L}{M-L}\rfloor$$

ここで、N_x は信号 x(n) の長さです。⌊⌋ 記号は床関数を表します。行列内の行数は、中央変換および両側変換の場合は DFT 点の数である N_DFT と同じで、実数値信号の片側変換の場合は N_DFT/2 に近い奇数に等しくなります。

STFT 行列 $$X(f)=\left[\begin{array}{ccccc}{X}_1(f)& X_2(f)& X_3(f)& \cdots & X_k(f)\end{array}\right]$$ の m 番目の列には、時間 mR 付近を中心としたウィンドウが適用されたデータの DFT が含まれます。

$$X_m(f)={\displaystyle \sum _{n=-\infty }^{\infty }x(n)g(n-mR)e^{-j2\pi fn}}.$$

• 短時間フーリエ変換は可逆変換です。逆変換プロセスは、ウィンドウが適用されたセグメントのオーバーラップ加算により、ウィンドウ エッジでの信号の減衰を補正します。詳細については、逆短時間フーリエ変換を参照してください。

• 関数 istft は、信号の STFT の逆変換を行います。

• 特定の条件下では、信号の "完全再構成" を実現することができます。詳細については、完全再構成を参照してください。

• 関数 stftmag2sig は、STFT の振幅から再構成された信号の推定を返します。

ヒント

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

参照

[1] Boashash, Boualem, ed. Time Frequency Signal Analysis and Processing: A Comprehensive Reference. Second edition. EURASIP and Academic Press Series in Signal and Image Processing. Amsterdam and Boston: Academic Press, 2016.

[2] 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.

[3] 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.

[4] Oppenheim, Alan V., and Ronald W. Schafer, with John R. Buck. Discrete-Time Signal Processing. Second edition. Upper Saddle River, NJ: Prentice Hall, 1999.

[5] Rabiner, Lawrence R., and Ronald W. Schafer. Digital Processing of Speech Signals. Englewood Cliffs, NJ: Prentice-Hall, 1978.

拡張機能

tall 配列

メモリの許容量を超えるような多数の行を含む配列を計算します。

使用上の注意および制限:

• 入力は tall 列ベクトルでなければなりません。

• 引数 window は常に指定しなければなりません。

• OutputTimeDimension は常に指定し、"downrows" に設定しなければなりません。

• reassigned オプションはサポートされません。

• 出力引数のない構文はサポートされません。

詳細については、tall 配列を参照してください。

C/C++ コード生成

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

使用上の注意および制限:

• 名前と値の引数を使用して指定する引数はコンパイル時の定数でなければなりません。

GPU コード生成

GPU Coder™ を使用して NVIDIA® GPU のための CUDA® コードを生成します。

使用上の注意および制限:

• 名前と値の引数を使用して指定する引数はコンパイル時の定数でなければなりません。

• 可変サイズの window は倍精度でなければなりません。

スレッドベースの環境

MATLAB® の backgroundPool を使用してバックグラウンドでコードを実行するか、Parallel Computing Toolbox™ の ThreadPool を使用してコードを高速化します。

使用上の注意および制限:

• 出力引数のない構文はサポートされません。

• 周波数ベクトルは等間隔でなければなりません。

詳細については、スレッドベースの環境での MATLAB 関数の実行を参照してください。

GPU 配列

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

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

バージョン履歴

R2006a より前に導入

R2024a: 単精度可変サイズ ウィンドウ入力のためのコード生成サポート

関数 spectrogram は、コード生成用の単精度可変サイズ ウィンドウ入力をサポートします。

R2023b: spectrogram は単精度データと GPU コード生成をサポートする

関数 spectrogram は、単精度入力およびグラフィックス処理装置 (GPU) 用のコード生成をサポートします。CUDA^® コードを生成するには、MATLAB^® Coder™ および GPU Coder™ が必要です。

R2023a: [プロットの作成] ライブ エディター タスクを使用した関数出力の可視化

[プロットの作成] ライブ エディター タスクを使用して、spectrogram の出力を対話的に可視化することができるようになりました。さまざまなグラフ タイプを選択し、オプションのパラメーターを設定できます。タスクは、ライブ スクリプトの一部となるコードも自動的に生成します。

参考

アプリ

• 信号アナライザー

関数

• goertzel | istft | periodogram | pspectrum | pwelch | stft | xspectrogram

トピック

• Signal Processing Toolbox を使用したスペクトログラムの計算
• 時間-周波数ギャラリー
• 時間-周波数解析の実践的基礎
• LPC 係数によるフォルマント推定

外部の Web サイト

• Fourier Analysis (MathWorks Teaching Resources)