Main Content

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

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

すべて折りたたむ

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

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

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

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

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

view(-45,65)
colormap jet

1 kHz で 2 秒間サンプリングされた 2 次チャープを生成します。瞬時周波数は、t=0 では 100 Hz であり、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);

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)

信号の片側、両側、および中心の短時間フーリエ変換を計算します。いずれの場合も、形状係数 β=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

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

入力引数

すべて折りたたむ

入力信号。ベクトル、行列、または 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
複素数のサポート: あり

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

データ型: double | single

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

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

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

データ型: double | single

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

データ型: double | single

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

データ型: double | single

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

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

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

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

    メモ

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

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

データ型: char | string

出力時間次元。'OutputTimeDimension' と、'acrosscolumns' または 'downrows' から構成されるコンマ区切りのペアとして指定します。行に沿った s の時間次元と、列に沿った周波数次元が必要な場合は、この値を 'downrows' に設定します。列に沿った s の時間次元と、行に沿った周波数次元が必要な場合は、この値を 'acrosscolumns' に設定します。関数が出力引数なしで呼び出された場合、この入力は無視されます。

出力引数

すべて折りたたむ

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

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

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

データ型: double | single

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

データ型: double | single

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

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

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

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

データ型: double | single

詳細

すべて折りたたむ

短時間フーリエ変換

短時間フーリエ変換 (STFT) を使用して、非定常信号の周波数成分が時間の経過と共に変化する様子を解析します。

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

k=NxLML,

によって指定されます。

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

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

Xm(f)=n=x(n)g(nmR)ej2πfn,

である X(f)=[X1(f)X2(f)X3(f)Xk(f)] によって指定されます。

ここで、

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

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

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

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

完全再構成

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

  • 入力サイズ — istft を使用して stft の出力を逆変換して、その結果を入力信号 x と同じ長さにする場合、k = (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.

拡張機能

R2019a で導入