stft

短時間フーリエ変換

構文

s = stft(x)

s = stft(x,fs)

s = stft(x,ts)

s = stft(___,Name=Value)

[s,f] = stft(___)

[s,f,t] = stft(___)

stft(___)

説明

s = stft(x) は x の短時間フーリエ変換 (STFT) を返します。

s = stft(x,fs) は、サンプル レート fs を使用して x の STFT を返します。

s = stft(x,ts) は、サンプル時間 ts を使用して x の STFT を返します。

s = stft(___,Name=Value) は、名前と値の引数を使用して追加オプションを指定します。オプションには FFT のウィンドウと長さが含まれます。これらの引数を前の入力構文のいずれかに追加できます。

[s,f] = stft(___) は、STFT が評価される周波数 f を返します。

[s,f,t] = stft(___) は、STFT が評価される時間を返します。

出力引数を設定せずに stft(___) を使用すると、現在の Figure ウィンドウに STFT の振幅の 2 乗がデシベル単位でプロットされます。

例

短時間フーリエ変換

正弦関数的に変化する周波数をもつチャープを生成します。この信号は 10 kHz で 2 秒間サンプリングされます。

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

チャープの短時間フーリエ変換を計算します。信号を 256 サンプルのセグメントに分割し、形状パラメーター β = 5 のカイザー ウィンドウを使用して各セグメントにウィンドウを適用します。隣接するセグメント間の 220 サンプルのオーバーラップおよび 512 の DFT 長を指定します。STFT が計算される周波数と時間の値を出力します。

[s,f,t] = stft(x,fs,Window=kaiser(256,5),OverlapLength=220,FFTLength=512);

STFT の振幅の 2 乗は、"スペクトログラム" とも呼ばれます。スペクトログラムをデシベル単位でプロットします。60 dB の範囲をもち、最後の色がスペクトログラムの最大値に対応するカラーマップを指定します。

sdb = mag2db(abs(s));
mesh(t,f/1000,sdb);

cc = max(sdb(:))+[-60 0];
ax = gca;
ax.CLim = cc;
view(2)
colorbar

出力引数なしで関数 stft を呼び出して、同じプロットを取得します。

stft(x,fs,Window=kaiser(256,5),OverlapLength=220,FFTLength=512)

2 次チャープの STFT

1 kHz で 2 秒間サンプリングされた 2 次チャープを生成します。瞬時周波数は、$\mathit{t}=0$ では 100 Hz であり、$\mathit{t}=1$ 秒で 200 Hz になります。

ts = 0:1/1e3:2;

f0 = 100;
f1 = 200;

x = chirp(ts,f0,1,f1,"quadratic",[],"concave");

持続時間が 1 ms の 2 次チャープの STFT を計算して表示します。

d = seconds(1e-3);
win = hamming(100,"periodic");

stft(x,d,Window=win,OverlapLength=98,FFTLength=128);

関数 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

STFT 周波数範囲

5 kHz で 4 秒間サンプリングされた信号を生成します。信号は、振動する振幅の領域および増加トレンドで変動する周波数の領域によって分離された、持続時間が減少していく一連のパルスで構成されています。信号をプロットします。

fs = 5000;
t = 0:1/fs:4-1/fs;

x = besselj(0,600*(sin(2*pi*(t+1).^3/30).^5));

plot(t,x)

信号の片側、両側、および中心の短時間フーリエ変換を計算します。いずれの場合も、形状係数 $\beta =10$ で 202 サンプルのカイザー ウィンドウを使用して、信号セグメントにウィンドウを適用します。各変換の計算に使用される周波数範囲を表示します。

rngs = ["onesided" "twosided" "centered"];

for kj = 1:length(rngs)
    
    opts = {'Window',kaiser(202,10),'FrequencyRange',rngs(kj)};

    [~,f] = stft(x,fs,opts{:});
    subplot(length(rngs),1,kj)
    stft(x,fs,opts{:})
    title(sprintf('''%s'': [%5.3f, %5.3f] kHz',rngs(kj),[f(1) f(end)]/1000))

end

計算を繰り返しますが、今回はカイザー ウィンドウの長さを奇数の 203 に変更します。'twosided' の周波数範囲は変わりません。他の 2 つの周波数範囲は、上限でオープンになります。

for kj = 1:length(rngs)
    
    opts = {'Window',kaiser(203,10),'FrequencyRange',rngs(kj)};

    [~,f] = stft(x,fs,opts{:});
    subplot(length(rngs),1,kj)
    stft(x,fs,opts{:})
    title(sprintf('''%s'': [%5.3f, %5.3f] kHz',rngs(kj),[f(1) f(end)]/1000))

end

マルチチャネル信号の STFT

1 kHz で 1 秒間サンプリングされた 3 つの異なるチャープで構成される 3 チャネルの信号を生成します。

1. 最初のチャネルは凹二次チャープで構成され、瞬時周波数は t = 0 では 100 Hz であり、t = 1 秒では 300 Hz になります。初期位相は 45 度です。

2. 2 番目のチャネルは凸二次チャープで構成され、瞬時周波数は t = 0 では 100 Hz であり、t = 1 秒では 500 Hz になります。

3. 3 番目のチャネルは対数チャープで構成され、瞬時周波数は t = 0 では 300 Hz であり、t = 1 秒では 500 Hz になります。

長さが 128 の周期的ハミング ウィンドウと 50 サンプルのオーバーラップ長を使用して、マルチチャネル信号の STFT を計算します。

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

x = [chirp(t,100,1,300,'quadratic',45,'concave');
      chirp(t,100,1,500,'quadratic',[],'convex');
      chirp(t,300,1,500,'logarithmic')]'; 
  
[S,F,T] = stft(x,fs,'Window',hamming(128,'periodic'),'OverlapLength',50);

各チャネルの STFT をウォーターフォール プロットとして可視化します。補助関数 helperGraphicsOpt を使用して、座標軸の動作を制御します。

waterfall(F,T,abs(S(:,:,1))')
helperGraphicsOpt(1)

waterfall(F,T,abs(S(:,:,2))')
helperGraphicsOpt(2)

waterfall(F,T,abs(S(:,:,3))')
helperGraphicsOpt(3)

この補助関数は、現在の座標軸の外観や動作を設定します。

function helperGraphicsOpt(ChannelId)
ax = gca;
ax.XDir = 'reverse';
ax.ZLim = [0 30];
ax.Title.String = ['Input Channel: ' num2str(ChannelId)];
ax.XLabel.String = 'Frequency (Hz)';
ax.YLabel.String = 'Time (seconds)';
ax.View = [30 45];
end

入力引数

x — 入力信号

ベクトル | 行列 | timetable

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

メモ

x と s を同じ長さになるようにする場合、(length(x)-noverlap)/(length(window)-noverlap) の値は整数でなければなりません。Window を使用して window の長さを指定し、OverlapLength を使用して noverlap を指定します。

• 入力に複数のチャネルがある場合は、各列がチャネルに対応する行列として x を指定します。

• timetable 入力の場合、x に等間隔に増加する有限の行時間を含めなければなりません。timetable が欠損している場合や時間点が重複している場合、欠損または重複する時間および非等間隔の時間をもつ timetable の整理のヒントを使用して修正できます。

• マルチチャネル timetable 入力の場合は、行列を含む単一の変数を使用した timetable、またはそれぞれが列ベクトルを含む複数の変数を使用した timetable として x を指定します。すべての変数は同じ精度でなければなりません。

x の各チャネルの長さは、ウィンドウの長さ以上でなければなりません。

例: chirp(0:1/4e3:2,250,1,500,"quadratic") は、単一チャネルのチャープを指定します。

例: timetable(rand(5,2),SampleRate=1) は 1 Hz で 4 秒間サンプリングされた 2 チャネルの確率変数を指定します。

例: timetable(rand(5,1),rand(5,1),SampleRate=1) は 1 Hz で 4 秒間サンプリングされた 2 チャネルの確率変数を指定します。

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

fs — サンプル レート

2π (既定値) | 正のスカラー

サンプル レート。正のスカラーで指定します。この引数は、x がベクトルまたは行列の場合にのみ適用されます。

データ型: double | single

ts — サンプル時間

duration スカラー

サンプル時間。duration スカラーで指定します。この引数は、x がベクトルまたは行列の場合にのみ適用されます。

例: seconds(1) は、連続する信号サンプル間の 1 秒間の時間差を表す duration スカラーです。

データ型: duration

名前と値の引数

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

例: Window=hamming(100),OverlapLength=50,FFTLength=128 は、隣接するセグメント間で 50 サンプルがオーバーラップし 128 点の FFT をもつ 100 サンプルのハミング ウィンドウを使用してデータをウィンドウ処理します。

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

例: 'Window',hamming(100),'OverlapLength',50,'FFTLength',128 は、隣接するセグメント間で 50 サンプルがオーバーラップし 128 点の FFT をもつ 100 サンプルのハミング ウィンドウを使用してデータをウィンドウ処理します。

Window — スペクトル ウィンドウ

hann(128,"periodic") (既定値) | ベクトル

スペクトル ウィンドウ。ベクトルとして指定します。ウィンドウを指定しない場合、またはウィンドウを空として指定する場合、関数は長さが 128 のハン ウィンドウを使用します。Window の長さは 2 以上でなければなりません。

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

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

データ型: double | single

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

ウィンドウの長さの 75% (既定値) | 非負の整数

Window の長さより小さな正の整数として指定する、オーバーラップするサンプル数。OverlapLength を省略するか、空として指定する場合、ウィンドウの長さの 75% より少ない最大整数に設定されます。これは、既定のハン ウィンドウの 96 サンプルです。

データ型: double | single

FFTLength — DFT 点の数

128 (既定値) | 正の整数

正の整数として指定する DFT 点の数。この値は、ウィンドウの長さ以上でなければなりません。入力信号の長さが DFT 長より小さい場合、データにはゼロがパディングされます。

データ型: double | single

FrequencyRange — STFT 周波数範囲

"centered" (既定値) | "twosided" | "onesided"

STFT 周波数範囲。"centered"、"twosided"、または "onesided" として指定します。

• "centered" — 中央揃えの両側 STFT を計算します。FFTLength が偶数の場合、s の計算区間は (–π, π] ラジアン/サンプルです。FFTLength が奇数の場合、s の計算区間は (–π, π) ラジアン/サンプルです。時間情報を指定すると、計算区間はそれぞれ (–f_s, f_s/2] サイクル/単位時間、(–f_s, f_s/2) サイクル/単位時間となります。ここで、f_s は有効なサンプル レートです。

• "twosided" — 区間 [0, 2π) ラジアン/サンプルで両側 STFT を計算します。時間情報を指定した場合、計算区間は [0, f_s) サイクル/単位時間となります。

• "onesided" — 片側 STFT を計算します。FFTLength が偶数の場合、s の計算区間は [0, π] ラジアン/サンプルです。FFTLength が奇数の場合、s の計算区間は [0, π) ラジアン/サンプルです。時間情報を指定すると、計算区間はそれぞれ [0, f_s/2] サイクル/単位時間、[0, f_s/2) サイクル/単位時間となります。ここで、f_s は有効なサンプル レートです。このオプションは、実信号でのみ有効です。

  メモ

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

例については、STFT 周波数範囲を参照してください。

データ型: char | string

OutputTimeDimension — 出力時間次元

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

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

出力引数

s — 短時間フーリエ変換

行列 | 3 次元配列

短時間フーリエ変換。行列または 3 次元配列として返されます。時間は s の列方向に、周波数は行方向に下がって増加します。3 番目の次元が存在する場合は、入力チャネルに対応します。

• 信号 x に N_x 個の時間サンプルがある場合、s には k 列があります。ここで、k = ⌊(N_x–L)/(M–L)⌋です。M は Window の長さで、L は OverlapLength であり、⌊ ⌋ 記号は床関数を表します。

• s の行数は FFTLength で指定された値と等しくなります。

データ型: double | single

f — 周波数

ベクトル

STFT が評価される周波数。ベクトルとして返されます。

データ型: double | single

t — 時点

ベクトル

時点。ベクトルとして返されます。t には、短時間のパワー スペクトルの推定の計算に使用されるデータ セグメントの中心に対応する時間値が含まれます。

• サンプル レート fs が指定された場合、ベクトルには時間値 (秒単位) が含まれます。

• サンプル時間 ts が指定された場合、ベクトルは、入力と同じ時間形式をもつ duration 配列です。

• 時間情報が指定されない場合、ベクトルはサンプル数を含みます。

データ型: double | single

詳細

短時間フーリエ変換

短時間フーリエ変換 (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}}.$$

参照

[1] Mitra, Sanjit K. Digital Signal Processing: A Computer-Based Approach. 2nd Ed. New York: McGraw-Hill, 2001.

[2] Sharpe, Bruce. Invertibility of Overlap-Add Processing. https://gauss256.github.io/blog/cola.html, accessed July 2019.

[3] Smith, Julius Orion. Spectral Audio Signal Processing. https://ccrma.stanford.edu/~jos/sasp/, online book, 2011 edition, accessed Nov 2018.

拡張機能

tall 配列

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

stft 関数は tall 配列をサポートしますが、次の使用上の注意と制限事項が伴います。

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

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

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

GPU コード生成

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

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

コード生成では、timetable はサポートされていません。

スレッドベースの環境

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

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

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

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

GPU 配列

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

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

バージョン履歴

R2019a で導入

参考

関数

• dlstft | istft | iscola | pspectrum | spectrogram | stftmag2sig

オブジェクト

• stftLayer

トピック

• Signal Processing Toolbox を使用したスペクトログラムの計算
• 時間-周波数ギャラリー
• 時間-周波数解析の実践的基礎