このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。
Signal Processing Toolbox を使用したスペクトログラムの計算
Signal Processing Toolbox™ には、非定常信号のスペクトログラムを計算する 3 つの関数が用意されています。各関数には、異なる入力引数、既定値、および出力があります。最適な選択は、特定の用途によって異なります。
"非定常" 信号は、周波数成分が時間の経過とともに変化する信号です。"短時間フーリエ変換" (STFT) は、この周波数成分が信号の変化に伴って変化する様子を解析するために使用します。STFT の振幅の 2 乗は、信号の "スペクトログラム" 時間-周波数表現と呼ばれます。
スペクトログラムは、時間-周波数表現の候補の 1 つにすぎません。Signal Processing Toolbox および Wavelet Toolbox™ で利用可能な他の時間-周波数表現の概要については、時間-周波数ギャラリーを参照してください。関数 periodogram
を使用した定常信号の処理については、FFT を使用したパワー スペクトル密度推定を参照してください。
スペクトログラム計算用の関数
Signal Processing Toolbox には、スペクトログラムの計算に使用できる次の関数があります。
spectrogram
— 最大限の柔軟性を目指して設計されています。STFT と、セグメントごとのパワー スペクトル密度またはパワー スペクトルを計算します。再割り当てをサポートします。stft
— 可逆性と最大限の制御を目指して設計されています。STFT を計算します。dlstft
およびstftLayer
によって使用されます。マルチチャネル入力をサポートします。これに対応する関数istft
は、逆 STFT を計算します。pspectrum
— 使いやすさを目指して設計されています。パワー スペクトルを計算します。信号アナライザーによって生成された解析スクリプトで使用されます。定常信号のスペクトルとパーシステンス スペクトルを計算できます。再割り当てをサポートします。
関数 spectrogram
は、以下の説明で参考として使用されます。
カテゴリ | パラメーター | 関数 | ||
---|---|---|---|---|
spectrogram | stft | pspectrum | ||
入力 | Signal | Nx 個の要素をもつベクトル |
|
|
ウィンドウ (g(n)) | 2 番目の位置引数 (既定: ハミング ウィンドウ) | 名前と値の引数 (既定: 周期的ハン ウィンドウ) | カイザー ウィンドウのみ | |
ウィンドウの長さ (M) | サンプル数として指定 (既定: ⌊Nx/4.5⌋) | サンプル数として指定 (既定: 128) | 名前と値の引数 TimeResolution | |
漏れ (ℓ) |
|
| カイザー ウィンドウの形状係数 β に関連付けられた名前と値の引数 Leakage : ℓ = 1 – β/40 | |
オーバーラップ (L) | 3 番目の位置引数として指定されるサンプル数 (既定: ウィンドウの長さの 50%) | 名前と値の引数 (既定: ウィンドウの長さの 75%) | 名前と値の引数 (既定: ここで、ENBW はウィンドウの等価ノイズ帯域幅) | |
DFT 点の数 (NDFT) | 4 番目の位置引数 (既定: max(256,2⌈log2M⌉)) | 名前と値の引数 (既定: 128) | 常に 1024 | |
時間情報 | 5 番目の位置引数として指定されるサンプル レート | 2 番目の位置引数として指定されるサンプル レートまたは時間ベクトル | 2 番目の位置引数として指定されるサンプル レートまたは時間ベクトル | |
関数呼び出し | fs = 100; x = exp(2j*pi*20* ... (0:1/fs:2-1/fs)); M = 200; lk = 0.5; g = kaiser(M,40*(1-lk)); L = 100; Ndft = 1024; | sps = abs( ... spectrogram(x,g,L, ... Ndft,fs,"centered") ... ).^2; | sts = abs( ... stft(x,fs,Window=g, ... OverlapLength=L, ... FFTLength=Ndft) ... ).^2; | pss = pspectrum(x,fs, ... "spectrogram", ... TimeResolution=M/fs, ... Leakage=lk, ... OverlapPercent=L/M*100 ... )*sum(g)^2; |
簡易プロット |
fs = 6e2; ts = 0:1/fs:2.05; x = vco(sin(2*pi*ts).* ... exp(-ts),[0.1 0.4]*fs,fs); M = 32; lk = 0.9; g = kaiser(M,40*(1-lk)); L = 22; Ndft = 1024;
| | | |
出力 | 周波数範囲 |
| 名前と値の引数
| logical の名前と値の引数
|
時間間隔 |
|
|
| |
正規化 |
| 最初の出力引数は STFT。その振幅の 2 乗がスペクトログラム。 |
| |
PSD とパワー スペクトル |
| 出力引数は STFT。 |
| |
例 |
STFT とスペクトログラムの定義
信号の STFT は、信号上の長さ M の "解析ウィンドウ" g(n) をスライドして、ウィンドウが適用されたデータの各セグメントの離散フーリエ変換 (DFT) を計算することによって算出されます。ウィンドウは、R サンプルの間隔で元の信号を飛び越えます。これは、隣り合ったセグメント間の L = M – R 個のサンプルのオーバーラップに相当します。ほとんどのウィンドウ関数は、スペクトル リンギングを回避するためにエッジで小さくなります。ウィンドウが適用された各セグメントの DFT は、時間と周波数の各点の振幅と位相を含む複素数値行列に対して追加されます。STFT 行列は次の列数をもちます。
ここで、Nx は信号 x(n) の長さです。⌊⌋ 記号は床関数を表します。行列内の行数は、中央変換および両側変換の場合は DFT 点の数である NDFT と同じで、実数値信号の片側変換の場合は NDFT/2 に近い奇数に等しくなります。
STFT 行列 の m 番目の列には、時間 mR 付近を中心としたウィンドウが適用されたデータの DFT が含まれます。
関数 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 を計算します。
信号を、それぞれ サンプルの長さのセグメントに分割します。
隣接するセグメント間に 個のサンプルのオーバーラップを指定します。
最後の短いセグメントを破棄します。
各セグメントにバートレット ウィンドウを適用します。
点で各セグメントの離散フーリエ変換を評価します。既定では、複素数値信号の場合、
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 の定義
定義を使用して サンプルの信号の STFT を計算します。信号を 個のオーバーラップ セグメントに分割します。各セグメントにウィンドウを適用し、それぞれについて、 点での離散フーリエ変換を評価します。
[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
の既定の値を使用して計算を繰り返します。
信号を の長さのセグメントに分割します。ここで、 は信号の長さです。各セグメントにハミング ウィンドウを適用します。
セグメント間で 50% のオーバーラップを指定します。
FFT を計算するには、 点を使用します。正の正規化周波数のみ、スペクトログラムを計算します。
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 つの出力をプロットします。どちらの場合も、周波数軸を で除算します。stft
出力の場合、サンプル数を有効なサンプル レート で除算します。
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
と pspectrum
の比較
電圧制御発振器と 3 つのガウス原子で構成される信号を生成します。信号は 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 を計算します。
サンプル信号を、 ミリ秒の時間分解能に対応する長さ のサンプルのセグメントに分割します。
個のサンプル、または隣接するセグメント間の 20% のオーバーラップを指定します。
カイザー ウィンドウを使用して各セグメントにウィンドウを適用し、漏れには を指定します。
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
は常にカイザー ウィンドウを として使用します。漏れ とウィンドウの形状係数 は、 の関係にあります。pspectrum
は、離散フーリエ変換の計算に常に 個の点を使用します。両側または中央揃えの周波数範囲の変換を計算したい場合にこの数値を指定することができます。一方、実信号の既定値である片側変換の場合、spectrogram
は 個の点を使用します。あるいは、この例にあるように、変換を計算したい周波数のベクトルを指定することができます。信号を 個のセグメントに厳密に分割できない場合、
spectrogram
は信号を切り捨て、pspectrum
は信号をゼロでパディングして追加のセグメントを作成します。出力を等価にするには、最後のセグメントと時間ベクトルの最後の要素を削除します。spectrogram
は、振幅の 2 乗をスペクトログラムとする STFT を返します。pspectrum
は、あらかじめ係数 で除算してから 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 乗されており、正規化係数 を含みます。実信号の片側スペクトログラムの場合も、追加の係数 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 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 を計算してプロットします。
信号を、それぞれ サンプルの長さのセグメントに分割します。
隣接するセグメント間に 個のサンプルのオーバーラップを指定します。
最後の短いセグメントを破棄します。
各セグメントにフラットトップ ウィンドウを適用します。
点で各セグメントの離散フーリエ変換を評価します。このとき、点の数が奇数であることに注意します。
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")
定義を使用して両側スペクトログラムを計算します。隣接するセグメント間で サンプルがオーバーラップする サンプルのセグメントに信号を分割します。各セグメントにウィンドウを適用し、それぞれについて、 点での離散フーリエ変換を計算します。
[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 のゼロ周波数成分をスペクトルの中央にシフトします。が奇数の場合、周波数範囲は両側でオープンとなります。 が偶数の場合、周波数範囲は下端でオープン、上端でクローズとなります。
出力を比較してスペクトログラムを表示します。
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 で使用したのと同じ時間値を使用します。
が奇数の場合、片側 STFT は両側 STFT の最初の 行により構成されます。 が偶数の場合、片側 STFT は両側 STFT の最初の 行により構成されます。
が奇数の場合、周波数範囲はゼロ周波数でクローズ、ナイキスト周波数でオープンとなります。 が偶数の場合、周波数範囲は両側でクローズとなります。
出力を比較してスペクトログラムを表示します。実数値信号の場合、"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 とパワー スペクトルのセグメント化
計算を繰り返しますが、今度は正規化周波数で計算します。サンプル レートを として指定した場合、結果は同じになります。
[~,~,~,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
参考
アプリ
関数
pspectrum
|spectrogram
|stft
|istft
|xspectrogram