短時間フーリエ変換 ( STFT ) は、信号の時間経過に伴う変化に応じて、信号の局所的な部分の正弦波周波数と位相成分を決定するために使用されるフーリエ変換です。 [ 1 ] 実際には、STFT を計算する手順は、長い時間の信号を 同じ長さの短いセグメントに分割し、各短いセグメントで個別にフーリエ変換を計算することです。これにより、各短いセグメントのフーリエスペクトルが明らかになります。次に、通常、変化するスペクトルを時間の関数としてプロットします。これは、ソフトウェア無線 (SDR) ベースのスペクトル表示でよく使用されるスペクトログラム またはウォーターフォールプロット として知られています。SDR の全範囲をカバーするフル帯域幅表示では、一般的に高速フーリエ変換 (FFT) が使用されます。
「nineteenth century」という単語の短時間フーリエ変換(STFT)の結果を視覚化したスペクトログラム。ここでは、縦軸が周波数の増加、横軸が時間を示しています。右側の凡例は、色の濃さが密度の増加に伴って増すことを示しています。
順方向STFT
離散時間STFT 離散時間の場合、変換対象のデータはチャンクまたはフレーム(境界でのアーティファクトを低減するために、通常は互いに重なり合う)に分割されます。各チャンクはフーリエ変換され 、複素変換された結果が、各時間と周波数における振幅と位相を記録する行列に追加されます。これは次のように表すことができます。
S T F T { x [ n ] } ( m 、 ω ) ≡ X ( m 、 ω ) = ∑ n = − ∞ ∞ x [ n ] w [ n − m ] e − 私 ω n {\displaystyle \mathbf {STFT} \{x[n]\}(m,\omega )\equiv X(m,\omega )=\sum _{n=-\infty }^{\infty }x[n]w[nm]e^{-i\omega n}} 同様に、信号と共にx [ n ] {\displaystyle x[n]} そして窓w [ n ] {\displaystyle w[n]} この場合、mは離散的で ω は連続ですが、ほとんどの典型的なアプリケーションでは、STFT は 高速フーリエ変換を 使用してコンピュータ上で実行されるため、両方の変数は離散的で量子化されます 。
STFTの振幅 の二乗は、関数のパワースペクトル密度のスペクトログラム表現を与える。
スペクトログラム { x ( t ) } ( τ 、 ω ) ≡ | X ( τ 、 ω ) | 2 \displaystyle \operatorname {spectrogram} \{x(t)\}(\tau ,\omega )\equiv |X(\tau ,\omega )|^{2}} 重複するウィンドウを使用するフーリエ変換の一種である修正離散コサイン変換 (MDCT)も参照のこと。
スライディングDFT ωの数が少ない場合、またはSTFTをウィンドウのシフトmごとに評価する必要がある場合は、 スライディングDFT アルゴリズムを使用してSTFTをより効率的に評価できます。[ 2 ]
逆STFT STFTは可逆であり 、つまり、逆STFTによって変換から元の信号を復元できます。STFTを反転する最も広く受け入れられている方法は、オーバーラップ加算(OLA)法を使用することであり、これによりSTFT複素スペクトルの変更も可能になります。これにより 、オーバーラップ加算修正法と 呼ばれる汎用性の高い信号処理 方法が実現します[ 3 ] 。
連続時間STFT ウィンドウ関数w ( t ) の幅と定義を考慮すると、まずウィンドウ関数の領域を次のようにスケーリングする必要があります。
∫ − ∞ ∞ w ( τ ) d τ = 1. {\displaystyle \int _{-\infty }^{\infty }w(\tau )\,d\tau =1.} 容易に次のことが導かれる
∫ − ∞ ∞ w ( t − τ ) d τ = 1 ∀ t {\displaystyle \int _{-\infty }^{\infty }w(t-\tau )\,d\tau =1\quad \forall \ t} そして
x ( t ) = x ( t ) ∫ − ∞ ∞ w ( t − τ ) d τ = ∫ − ∞ ∞ x ( t ) w ( t − τ ) d τ 。 {\displaystyle x(t)=x(t)\int _{-\infty }^{\infty }w(t-\tau )\,d\tau =\int _{-\infty }^{\infty }x(t)w(t-\tau )\,d\tau .} 連続フーリエ変換は
X ( ω ) = ∫ − ∞ ∞ x ( t ) e − 私 ω t d t 。 {\displaystyle X(\omega )=\int _{-\infty }^{\infty }x(t)e^{-i\omega t}\,dt.} 上記のx ( t )を代入すると次のようになります。
X ( ω ) = ∫ − ∞ ∞ [ ∫ − ∞ ∞ x ( t ) w ( t − τ ) d τ ] e − 私 ω t d t {\displaystyle X(\omega )=\int _{-\infty }^{\infty }\left[\int _{-\infty }^{\infty }x(t)w(t-\tau )\,d\tau \right]\,e^{-i\omega t}\,dt} = ∫ − ∞ ∞ ∫ − ∞ ∞ x ( t ) w ( t − τ ) e − 私 ω t d τ d t 。 {\displaystyle =\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }x(t)w(t-\tau )\,e^{-i\omega t}\,d\tau \,dt.} 積分順序の入れ替え:
X ( ω ) = ∫ − ∞ ∞ ∫ − ∞ ∞ x ( t ) w ( t − τ ) e − 私 ω t d t d τ {\displaystyle X(\omega )=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }x(t)w(t-\tau )\,e^{-i\omega t}\,dt\,d\tau } = ∫ − ∞ ∞ [ ∫ − ∞ ∞ x ( t ) w ( t − τ ) e − 私 ω t d t ] d τ {\displaystyle =\int _{-\infty }^{\infty }\left[\int _{-\infty }^{\infty }x(t)w(t-\tau )\,e^{-i\omega t}\,dt\right]\,d\tau } = ∫ − ∞ ∞ X ( τ 、 ω ) d τ 。 {\displaystyle =\int _{-\infty }^{\infty }X(\tau ,\omega )\,d\tau .} したがって、フーリエ変換は、 x ( t )のすべての STFT の位相コヒーレントな和のようなものと見なすことができます。逆フーリエ変換は
x ( t ) = 1 2 π ∫ − ∞ ∞ X ( ω ) e + 私 ω t d ω 、 {\displaystyle x(t)={\frac {1}{2\pi }}\int _{-\infty }^{\infty }X(\omega )e^{+i\omega t}\,d\omega ,} すると、X (τ,ω)からx ( t )を復元できる。
x ( t ) = 1 2 π ∫ − ∞ ∞ ∫ − ∞ ∞ X ( τ 、 ω ) e + 私 ω t d τ d ω 。 {\displaystyle x(t)={\frac {1}{2\pi }}\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }X(\tau ,\omega )e^{+i\omega t}\,d\tau \,d\omega .} または
x ( t ) = ∫ − ∞ ∞ [ 1 2 π ∫ − ∞ ∞ X ( τ 、 ω ) e + 私 ω t d ω ] d τ 。 {\displaystyle x(t)=\int _{-\infty }^{\infty }\left[{\frac {1}{2\pi }}\int _{-\infty }^{\infty }X(\tau ,\omega )e^{+i\omega t}\,d\omega \right]\,d\tau .} 上記と比較すると、x ( t )のウィンドウ化された「粒度」または「ウェーブレット」は
x ( t ) w ( t − τ ) = 1 2 π ∫ − ∞ ∞ X ( τ 、 ω ) e + 私 ω t d ω 。 {\displaystyle x(t)w(t-\tau )={\frac {1}{2\pi }}\int _{-\infty }^{\infty }X(\tau ,\omega )e^{+i\omega t}\,d\omega .} τ を固定した場合のX (τ,ω)の逆フーリエ変換。
τの近傍でのみ有効な別の定義として、逆変換は次のようになる。
x ( t ) = 1 w ( t − τ ) 1 2 π ∫ − ∞ ∞ X ( τ 、 ω ) e + 私 ω t d ω 。 {\displaystyle x(t)={\frac {1}{w(t-\tau )}}{\frac {1}{2\pi }}\int _{-\infty }^{\infty }X(\tau ,\omega )e^{+i\omega t}\,d\omega .} 一般的に、ウィンドウ関数w ( t ) {\displaystyle w(t)} 以下の特性を持つ。
(a)偶数対称性:w ( t ) = w ( − t ) {\displaystyle w(t)=w(-t)} ; (b)非増加(正の時間の場合):w ( t ) ≥ w ( s ) {\displaystyle w(t)\geq w(s)} もし| t | ≤ | s | {\displaystyle |t|\leq |s|} ; (c)コンパクトサポート:w ( t ) {\displaystyle w(t)} |t|が大きい場合、はゼロに等しくなります。
解決に関する問題 STFTの欠点の1つは、分解能が固定されていることです。ウィンドウ関数の幅は、信号の表現方法に関係しており、周波数分解能(近接する周波数成分を分離できるかどうか)が良いか、時間分解能(周波数が変化する時間)が良いかを決定します。ウィンドウ幅が広いほど周波数分解能は高くなりますが、時間分解能は低くなります。ウィンドウ幅が狭いほど時間分解能は高くなりますが、周波数分解能は低くなります。これらはそれぞれ、狭帯域 変換と広帯域変換と呼ばれます。
STFT解像度の比較。左側は時間分解能が高く、右側は周波数分解能が高い。 これがウェーブレット変換 とマルチレゾリューション解析 が開発された理由の一つであり、これらは高周波イベントに対して優れた時間分解能を、低周波イベントに対して優れた周波数分解能を提供することができ、多くの実信号に最適な組み合わせとなっている。
この特性はハイゼンベルクの 不確定性原理 と関連していますが、直接的な関係はありません。詳しくはガボール限界 を参照してください。時間と周波数の標準偏差 の積には制限があります。不確定性原理の境界(両方の同時解像度の最適値)は、ガウス窓関数(またはマスク関数)によって達成されます。これは、ガウス関数がフーリエの不確定性原理を 最小化するためです。これはガボール変換 と呼ばれ(多重解像度用に修正するとモルレーウェーブレット 変換になります)、この原理はガボール変換と呼ばれます。
以下の例に示すように、ウィンドウサイズを変化させた場合の短時間フーリエ変換(STFT)は、2次元領域(時間と周波数)として考えることができ、ウィンドウサイズを変化させることで計算できます。ただし、これは厳密には時間周波数表現ではなく、カーネルは信号全体にわたって一定ではありません。
例 元の関数が次のようになっている場合:
X ( t 、 f ) = ∫ − ∞ ∞ w ( t − τ ) x ( τ ) e − j 2 π f τ d τ {\displaystyle X(t,f)=\int _{-\infty }^{\infty }w(t-\tau )x(\tau )e^{-j2\pi f\tau }d\tau } 簡単な例を挙げてみましょう。
w(t) = 1 (|t|≦Bの場合)
w(t) = 0 それ以外の場合
B = 窓
短時間フーリエ変換の元の関数は次のように変更できます。
X ( t 、 f ) = ∫ t − B t + B x ( τ ) e − j 2 π f τ d τ {\displaystyle X(t,f)=\int _{t-B}^{t+B}x(\tau )e^{-j2\pi f\tau }d\tau } 別の例:
以下のサンプル信号を使用しますx ( t ) {\displaystyle x(t)} これは、4 つの正弦波形が順番に結合されたものです。各波形は、4 つの周波数 (10、25、50、100 Hz ) のうちの 1 つのみで構成されています。x ( t ) {\displaystyle x(t)} は:
x ( t ) = { コス ( 2 π 10 t ) 0 s ≤ t < 5 s コス ( 2 π 25 t ) 5 s ≤ t < 10 s コス ( 2 π 50 t ) 10 s ≤ t < 15 s コス ( 2 π 100 t ) 15 s ≤ t < 20 s {\displaystyle x(t)={\begin{cases}\cos(2\pi 10t)&0\,\mathrm {s} \leq t<5\,\mathrm {s} \\\cos(2\pi 25t)&5\,\mathrm {s} \leq t<10\,\mathrm {s} \\\cos(2\pi 50t)&10\,\mathrm {s} \leq t<15\,\mathrm {s} \\\cos(2\pi 100t)&15\,\mathrm {s} \leq t<20\,\mathrm {s} \\\end{cases}}} 次に、400Hzでサンプリングを行った 。その結果、以下のスペクトログラムが得られた。
25 ミリ秒のウィンドウを用いると、信号が変化する正確な時間を特定できますが、正確な周波数を特定するのは困難です。一方、1000ミリ秒のウィンドウを用いると、周波数を正確に把握できますが、周波数変化間の時間間隔が不明瞭になります。
その他の例:
w ( t ) = e x p ( σ − t 2 ) {\displaystyle w(t)=exp(\sigma -t^{2})} 通常、私たちはe x p ( σ − t 2 ) {\displaystyle exp(\sigma -t^{2})} ガウス関数 またはガボール関数。これを用いる場合、短時間フーリエ変換は「ガボール変換」と呼ばれる。
説明 サンプリングとナイキスト周波数 を参照して説明することもできます。
任意の実数値信号から、サンプリングレートf sで N 個 のサンプルからなるウィンドウを取ります。フーリエ変換を行うと、N 個 の複素係数が得られます。これらの係数のうち、実際に使用できるのは半分だけです(最後のN/2 は、最初のN/2の 複素共役を 逆順にしたものです。これは実数値信号であるためです)。
これらのN/2 個の係数は、周波数 0 からf s /2 (ナイキスト周波数) までを表し、連続する 2 つの係数は f s / N Hzの間隔で配置されます 。
ウィンドウの周波数分解能を高めるには、係数の周波数間隔を狭める必要があります。変数は2つだけですが、f sを小さくする( Nを 一定に保つ)と、単位時間あたりのサンプル数が少なくなるため、ウィンドウサイズが大きくなります。もう1つの方法はN を大きくすることですが、これもまたウィンドウサイズを大きくする原因となります。したがって、周波数分解能を高めようとすると、ウィンドウサイズが大きくなり、結果として時間分解能が低下します。そして、その逆もまた然りです。
実装 元の機能
X ( t 、 f ) = ∫ − ∞ ∞ w ( t − τ ) x ( τ ) e − j 2 π f τ d τ {\displaystyle X(t,f)=\int _{-\infty }^{\infty }w(t-\tau )x(\tau )e^{-j2\pi f\tau }d\tau } 離散形式に変換する:
t = n Δ t 、 f = m Δ f 、 τ = p Δ t {\displaystyle t=n\Delta _{t},f=m\Delta _{f},\tau =p\Delta _{t}} X ( n Δ t 、 m Δ f ) = ∑ − ∞ ∞ w ( ( n − p ) Δ t ) x ( p Δ t ) e − j 2 π p m Δ t Δ f Δ t {\displaystyle X(n\Delta _{t},m\Delta _{f})=\sum _{-\infty }^{\infty }w((n-p)\Delta _{t})x(p\Delta _{t})e^{-j2\pi pm\Delta _{t}\Delta _{f}}\Delta _{t}} 仮に
w ( t ) ≅ 0 のために | t | > B 、 B Δ t = Q {\displaystyle w(t)\cong 0{\text{ for }}|t|>B,{\frac {B}{\Delta _{t}}}=Q} すると、元の関数を次のように書き込むことができます。
X ( n Δ t 、 m Δ f ) = ∑ p = n − Q n + Q w ( ( n − p ) Δ t ) x ( p Δ t ) e − j 2 π p m Δ t Δ f Δ t {\displaystyle X(n\Delta _{t},m\Delta _{f})=\sum _{p=n-Q}^{n+Q}w((n-p)\Delta _{t})x(p\Delta _{t})e^{-j2\pi pm\Delta _{t}\Delta _{f}}\Delta _{t}}
直接実装
制約 a. ナイキスト基準(エイリアシング 効果の回避):
Δ t < 1 2 Ω {\displaystyle \Delta _{t}<{\frac {1}{2\Omega }}} 、 どこΩ {\displaystyle \Omega } 帯域幅はx ( τ ) w ( t − τ ) {\displaystyle x(\tau )w(t-\tau )}
再帰的方法
制約 a.Δ t Δ f = 1 N {\displaystyle \Delta _{t}\Delta _{f}={\tfrac {1}{N}}} 、 どこN {\displaystyle N} 整数です
b.N ≥ 2 Q + 1 {\displaystyle N\geq 2Q+1}
c. ナイキスト基準(エイリアシング効果の回避):
Δ t < 1 2 Ω {\displaystyle \Delta _{t}<{\frac {1}{2\Omega }}} 、Ω {\displaystyle \Omega } 帯域幅はx ( τ ) w ( t − τ ) {\displaystyle x(\tau )w(t-\tau )} d.矩形STFT を実装する場合のみ
長方形ウィンドウは制約を課す
w ( ( n − p ) Δ t ) = 1 {\displaystyle w((n-p)\Delta _{t})=1} 代入すると次のようになる。
X ( n Δ t 、 m Δ f ) = ∑ p = n − Q n + Q w ( ( n − p ) Δ t ) x ( p Δ t ) e − j 2 π p m N Δ t = ∑ p = n − Q n + Q x ( p Δ t ) e − j 2 π p m N Δ t {\displaystyle {\begin{aligned}X(n\Delta _{t},m\Delta _{f})&=\sum _{p=n-Q}^{n+Q}w((n-p)\Delta _{t})&x(p\Delta _{t})e^{-{\frac {j2\pi pm}{N}}}\Delta _{t}\\&=\sum _{p=n-Q}^{n+Q}&x(p\Delta _{t})e^{-{\frac {j2\pi pm}{N}}}\Delta _{t}\\\end{aligned}}} n に対するn -1 の変数変換:
X ( ( n − 1 ) Δ t 、 m Δ f ) = ∑ p = n − 1 − Q n − 1 + Q x ( p Δ t ) e − j 2 π p m N Δ t {\displaystyle X((n-1)\Delta _{t},m\Delta _{f})=\sum _{p=n-1-Q}^{n-1+Q}x(p\Delta _{t})e^{-{\frac {j2\pi pm}{N}}}\Delta _{t}} 計算するX ( ミニ n Δ t 、 m Δ f ) {\displaystyle X(\min {n}\Delta _{t},m\Delta _{f})} N 点FFTによる:
X ( n 0 Δ t 、 m Δ f ) = Δ t e j 2 π ( Q − n 0 ) m N ∑ q = 0 N − 1 x 1 ( q ) e − j 2 π q m N 、 n 0 = ミニ ( n ) {\displaystyle X(n_{0}\Delta _{t},m\Delta _{f})=\Delta _{t}e^{\frac {j2\pi (Q-n_{0})m}{N}}\sum _{q=0}^{N-1}x_{1}(q)e^{-j{\frac {2\pi qm}{N}}},\qquad n_{0}=\min {(n)}} どこ
x 1 ( q ) = { x ( ( n − Q + q ) Δ t ) q ≤ 2 Q 0 q > 2 Q {\displaystyle x_{1}(q)={\begin{cases}x((n-Q+q)\Delta _{t})&q\leq 2Q\\0&q>2Q\end{cases}}} 再帰式を適用して計算するX ( n Δ t 、 m Δ f ) {\displaystyle X(n\Delta _{t},m\Delta _{f})}
X ( n Δ t 、 m Δ f ) = X ( ( n − 1 ) Δ t 、 m Δ f ) − x ( ( n − Q − 1 ) Δ t ) e − j 2 π ( n − Q − 1 ) m N Δ t + x ( ( n + Q ) Δ t ) e − j 2 π ( n + Q ) m N Δ t {\displaystyle X(n\Delta _{t},m\Delta _{f})=X((n-1)\Delta _{t},m\Delta _{f})-x((n-Q-1)\Delta _{t})e^{-{\frac {j2\pi (n-Q-1)m}{N}}}\Delta _{t}+x((n+Q)\Delta _{t})e^{-{\frac {j2\pi (n+Q)m}{N}}}\Delta _{t}}
制約 exp ( − j 2 π p m Δ t Δ f ) = exp ( − j π p 2 Δ t Δ f ) ⋅ exp ( j π ( p − m ) 2 Δ t Δ f ) ⋅ exp ( − j π m 2 Δ t Δ f ) {\displaystyle \exp {(-j2\pi pm\Delta _{t}\Delta _{f})}=\exp {(-j\pi p^{2}\Delta _{t}\Delta _{f})}\cdot \exp {(j\pi (p-m)^{2}\Delta _{t}\Delta _{f})}\cdot \exp {(-j\pi m^{2}\Delta _{t}\Delta _{f})}} それで
X ( n Δ t 、 m Δ f ) = Δ t ∑ p = n − Q n + Q w ( ( n − p ) Δ t ) x ( p Δ t ) e − j 2 π p m Δ t Δ f {\displaystyle X(n\Delta _{t},m\Delta _{f})=\Delta _{t}\sum _{p=n-Q}^{n+Q}w((n-p)\Delta _{t})x(p\Delta _{t})e^{-j2\pi pm\Delta _{t}\Delta _{f}}} X ( n Δ t 、 m Δ f ) = Δ t e − j 2 π m 2 Δ t Δ f ∑ p = n − Q n + Q w ( ( n − p ) Δ t ) x ( p Δ t ) e − j π p 2 Δ t Δ f e j π ( p − m ) 2 Δ t Δ f {\displaystyle X(n\Delta _{t},m\Delta _{f})=\Delta _{t}e^{-j2\pi m^{2}\Delta _{t}\Delta _{f}}\sum _{p=n-Q}^{n+Q}w((n-p)\Delta _{t})x(p\Delta _{t})e^{-j\pi p^{2}\Delta _{t}\Delta _{f}}e^{j\pi (p-m)^{2}\Delta _{t}\Delta _{f}}}
参考文献 ↑ Sejdić E.; Djurović I.; Jiang J. (2009). "エネルギー集中を用いた時間周波数特徴表現:最近の進歩の概要". Digital Signal Processing . 19 (1): 153– 183. Bibcode : 2009DSP....19..153S . doi : 10.1016/j.dsp.2007.12.004 . ↑ E. Jacobsen および R. Lyons、「スライディング DFT」、 Signal Processing Magazine vol. 20、issue 2、pp. 74–80 (2003 年 3 月)。 ↑ Jont B. Allen (1977 年 6 月)「離散フーリエ変換による短時間スペクトル解析、合成、および修正」 IEEE Transactions on Acoustics, Speech, and Signal Processing . ASSP-25 (3): 235–238 . doi : 10.1109/TASSP.1977.1162950 . ↑ Kleinfeld, David; Mitra, Partha P. (2014年3月). "機能的脳イメージングのためのスペクトル法" . Cold Spring Harbor Protocols . 2014 (3): 248– 262. doi : 10.1101/pdb.top081075 . PMID 24591695 . ↑ 「要求された周波数分解能に対してパディングが不十分です」とはどういう意味ですか? – FieldTripツールボックス」 。 ↑ Zeitler M、Fries P 、Gielen S ( 2008)。 「 ガンマ振動の振幅の変動による偏った競合」 。J Comput Neurosci。25 ( 1 ) : 89–107。doi : 10.1007 /s10827-007-0066-2。PMC 2441488。PMID 18293071 。 ↑ ウィンガーデン、マリジン・ヴァン。ヴィンク、マーティン。ジャン・ランケルマ。 Pennartz、Cyriel MA (2010-05-19)。 「報酬期待時の眼窩前頭ニューロンのシータバンド位相ロック」 。 神経科学ジャーナル 。 30 (20): 7078–7087 。 土井 : 10.1523/JNEUROSCI.3860-09.2010 。 ISSN 0270-6474 。 PMC 6632657 。 PMID 20484650 。
外部リンク DiscreteTFDs – 短時間フーリエ変換およびその他の時間周波数分布を計算するためのソフトウェア 特異スペクトル解析 – マルチテーパー法ツールキット– 短くノイズの多い時系列データを解析するための無料ソフトウェアプログラム SpectraWorks社製 Mac OS X 用 kSpectra Toolkit 超広帯域信号の時間周波数解析のための時間伸長短時間フーリエ変換 STFTと逆STFTを実行するBSDライセンスのMatlabクラス LTFAT – 短時間フーリエ変換と時間周波数解析を行うための無料(GPLライセンス)のMatlab/Octaveツールボックス Sonogram visible speech – 短時間フーリエ変換と時間周波数解析のための無料(GPL)フリーウェア 国立台湾大学、時間周波数解析とウェーブレット変換、2021年、教授:丁建俊、電気工学科