このページの翻訳は最新ではありません。ここをクリックして、英語の最新版を参照してください。

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 の振幅がプロットされます。

例

3D STFT の可視化

10 kHz でサンプリングされる正弦波によって制御された、2 秒間の電圧制御発振器出力を生成します。

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

信号の STFT を計算してプロットします。長さ 256 で形状パラメーター $\beta =5$ のカイザー ウィンドウを使用します。オーバーラップの長さを 220 サンプル、DFT 長を 512 点として指定します。既定のカラーマップとビューをもつ STFT をプロットします。

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

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

view(-45,65)
colormap jet

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);

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 では 200 Hz であり、t = 1 秒では 600 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,200,1,600,'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 として指定します。

メモ

istft を使用して s を反転させて、その結果を x と同じ長さになるようにする場合、(length(x)-noverlap)/(length(window)-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

名前と値のペアの引数

オプションの引数 Name,Value のコンマ区切りペアを指定します。Name は引数名で、Value は対応する値です。Name は引用符で囲まなければなりません。Name1,Value1,...,NameN,ValueN のように、複数の名前と値のペアの引数を、任意の順番で指定できます。

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

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

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

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

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

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

データ型: double | single

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

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

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

データ型: double | single

'FFTLength' — DFT 点の数

128 (既定値) | 正の整数

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

データ型: double | single

'FrequencyRange' — STFT 周波数範囲

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

STFT 周波数範囲。'FrequencyRange' と、'centered'、'twosided' または 'onesided' で構成されるコンマ区切りのペアとして指定します。

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

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

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

  メモ

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

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

データ型: char | string

'OutputTimeDimension' — 出力時間次元

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

出力時間次元。'OutputTimeDimension' と、'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 は、信号上の長さ $$M$$ の "解析ウィンドウ" をスライドして、ウィンドウが適用されたデータの離散フーリエ変換を計算することによって計算されます。ウィンドウは、$$R$$ サンプルの間隔で元の信号を飛び越えます。ほとんどのウィンドウ関数は、スペクトル リンギングを回避するためにエッジで小さくなります。非ゼロのオーバーラップ長 $$L$$ が指定されている場合、ウィンドウが適用されたセグメントのオーバーラップ加算がウィンドウ エッジでの信号の減衰を補正します。ウィンドウが適用された各セグメントの DFT は、時間と周波数の各点の振幅と位相を含む行列に対して追加されます。STFT 配列の行数は DFT 点の数と同じで、列数は

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

によって指定されます。

ここで、$$N_x$$ は元の信号 $$x(n)$$ の長さです。⌊⌋ 記号は床関数を表します。

STFT 行列は、この行列の $$m$$ 番目の要素が

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

である $$X(f)=\left[\begin{array}{ccccc}{X}_1(f)& X_2(f)& X_3(f)& \cdots & X_k(f)\end{array}\right]$$ によって指定されます。

ここで、

• $$g(n)$$ — 長さ $$M$$ のウィンドウ関数。

• $$X_m(f)$$ — 時間 $$mR$$ 付近を中心としたウィンドウが適用されたデータの DFT。

• $$R$$ — 連続する DFT 間のホップ サイズ。ホップ サイズは、ウィンドウの長さ $$M$$ とオーバーラップ長 $$L$$ 間の差異です。

STFT の振幅二乗では、関数のパワー スペクトル密度の spectrogram 表現が得られます。

完全再構成

一般に、入力信号の STFT を計算して逆変換しても、完全再構成は行われません。ISTFT の出力を元の入力信号と可能な限り一致させるには、信号とウィンドウが以下の条件を満たす必要があります。

• 入力サイズ — istft を使用して stft の出力を逆変換して、その結果を入力信号 x と同じ長さにする場合、$$k=\frac{(length(x)-noverlap)}{(length(window)-noverlap)}$$ の値は整数でなければなりません。

• COLA 準拠 — COLA 準拠ウィンドウを使用して、信号の短時間フーリエ変換を変更していないと仮定します。

• パディング — k の値が整数にならないような入力信号の長さの場合、短時間フーリエ変換を計算する前に信号をゼロ パディングします。信号を逆変換した後に、余分なゼロを削除します。

参照

[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 配列

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

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

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

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

C/C++ コード生成

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

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

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

GPU コード生成

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

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

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

GPU 配列

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

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

参考

関数

• iscola | istft | pspectrum

トピック

• 時間-周波数ギャラリー

R2019a で導入

×

例を開く

この例の変更されたバージョンがあります。編集された方の例を開きますか?

いいえ、変更されたバージョンを上書きします。 はい