Short-time Fourier transform

Simply, in the continuous-time case, the function to be transformed is multiplied by a window function which is nonzero for only a short period of time. The Fourier transform (a one-dimensional function) of the resulting signal is taken, then the window is slid along the time axis until the end resulting in a two-dimensional representation of the signal. Mathematically, this is written as:

${\displaystyle \mathbf {STFT} \{x(t)\}(\tau ,\omega )\equiv X(\tau ,\omega )=\int _{-\infty }^{\infty }x(t)w(t-\tau )e^{-i\omega t}\,dt}$

where ${\displaystyle w(\tau )}$ is the window function, commonly a Hann window or Gaussian window centered around zero, and ${\displaystyle x(t)}$ is the signal to be transformed (note the difference between the window function ${\displaystyle w}$ and the frequency ${\displaystyle \omega }$). ${\displaystyle X(\tau ,\omega )}$ is essentially the Fourier transform of ${\displaystyle x(t)w(t-\tau )}$, a complex function representing the phase and magnitude of the signal over time and frequency. Often phase unwrapping is employed along either or both the time axis, ${\displaystyle \tau }$, and frequency axis, ${\displaystyle \omega }$, to suppress any jump discontinuity of the phase result of the STFT. The time index ${\displaystyle \tau }$ is normally considered to be "slow" time and usually not expressed in as high resolution as time ${\displaystyle t}$. Given that the STFT is essentially a Fourier transform times a window function, the STFT is also called windowed Fourier transform or time-dependent Fourier transform.

In the discrete time case, the data to be transformed could be broken up into chunks or frames (which usually overlap each other, to reduce artifacts at the boundary). Each chunk is Fourier transformed, and the complex result is added to a matrix, which records magnitude and phase for each point in time and frequency. This can be expressed as:

${\displaystyle \mathbf {STFT} \{x[n]\}(m,\omega )\equiv X(m,\omega )=\sum _{n=-\infty }^{\infty }x[n]w[n-m]e^{-i\omega n}}$

likewise, with signal x[n] and window w[n]. In this case, m is discrete and ω is continuous, but in most typical applications the STFT is performed on a computer using the fast Fourier transform, so both variables are discrete and quantized.

The magnitude squared of the STFT yields the spectrogram representation of the power spectral density of the function:

${\displaystyle \operatorname {spectrogram} \{x(t)\}(\tau ,\omega )\equiv |X(\tau ,\omega )|^{2}}$

See also the modified discrete cosine transform (MDCT), which is also a Fourier-related transform that uses overlapping windows.

If only a small number of ω are desired, or if the STFT is desired to be evaluated for every shift m of the window, then the STFT may be more efficiently evaluated using a sliding DFT algorithm.[2]

Given the width and definition of the window function w(t), we initially require the area of the window function to be scaled so that

${\displaystyle \int _{-\infty }^{\infty }w(\tau )\,d\tau =1.}$

It easily follows that

${\displaystyle \int _{-\infty }^{\infty }w(t-\tau )\,d\tau =1\quad \forall \ t}$

and

${\displaystyle x(t)=x(t)\int _{-\infty }^{\infty }w(t-\tau )\,d\tau =\int _{-\infty }^{\infty }x(t)w(t-\tau )\,d\tau .}$

The continuous Fourier transform is

${\displaystyle X(\omega )=\int _{-\infty }^{\infty }x(t)e^{-i\omega t}\,dt.}$

Substituting x(t) from above:

${\displaystyle X(\omega )=\int _{-\infty }^{\infty }\left[\int _{-\infty }^{\infty }x(t)w(t-\tau )\,d\tau \right]\,e^{-i\omega t}\,dt}$

${\displaystyle =\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }x(t)w(t-\tau )\,e^{-i\omega t}\,d\tau \,dt.}$

Swapping order of integration:

${\displaystyle X(\omega )=\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }x(t)w(t-\tau )\,e^{-i\omega t}\,dt\,d\tau }$

${\displaystyle =\int _{-\infty }^{\infty }\left[\int _{-\infty }^{\infty }x(t)w(t-\tau )\,e^{-i\omega t}\,dt\right]\,d\tau }$

${\displaystyle =\int _{-\infty }^{\infty }X(\tau ,\omega )\,d\tau .}$

So the Fourier transform can be seen as a sort of phase coherent sum of all of the STFTs of x(t). Since the inverse Fourier transform is

${\displaystyle x(t)={\frac {1}{2\pi }}\int _{-\infty }^{\infty }X(\omega )e^{+i\omega t}\,d\omega ,}$

then x(t) can be recovered from X(τ,ω) as

${\displaystyle x(t)={\frac {1}{2\pi }}\int _{-\infty }^{\infty }\int _{-\infty }^{\infty }X(\tau ,\omega )e^{+i\omega t}\,d\tau \,d\omega .}$

or

${\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 .}$

It can be seen, comparing to above that windowed "grain" or "wavelet" of x(t) is

${\displaystyle x(t)w(t-\tau )={\frac {1}{2\pi }}\int _{-\infty }^{\infty }X(\tau ,\omega )e^{+i\omega t}\,d\omega .}$

the inverse Fourier transform of X(τ,ω) for τ fixed.

An alternative definition that is valid only in the vicinity of τ, the inverse transform is:

${\displaystyle x(t)={\frac {1}{w(t-\tau )}}{\frac {1}{2\pi }}\int _{-\infty }^{\infty }X(\tau ,\omega )e^{+i\omega t}\,d\omega .}$

In general, the window function ${\displaystyle w(t)}$ has the following properties:

(a) even symmetry: ${\displaystyle w(t)=w(-t)}$;

(b) non-increasing (for positive time): ${\displaystyle w(t)\geq w(s)}$ if ${\displaystyle |t|\leq |s|}$;

(c) compact support: ${\displaystyle w(t)}$ is equal to zero when |t| is large.

When the original function is:

${\displaystyle X(t,f)=\int _{-\infty }^{\infty }w(t-\tau )x(\tau )e^{-j2\pi f\tau }d\tau }$

We can have a simple example:

w(t) = 1 for |t| smaller than or equal B

w(t) = 0 otherwise

B = window

Now the original function of the Short-time Fourier transform can be changed as

${\displaystyle X(t,f)=\int _{t-B}^{t+B}x(\tau )e^{-j2\pi f\tau }d\tau }$

Another example:

Using the following sample signal ${\displaystyle x(t)}$ that is composed of a set of four sinusoidal waveforms joined together in sequence. Each waveform is only composed of one of four frequencies (10, 25, 50, 100 Hz). The definition of ${\displaystyle x(t)}$ is:

${\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}}}$

Then it is sampled at 400 Hz. The following spectrograms were produced:

25 ms window	125 ms window
375 ms window	1000 ms window

The 25 ms window allows us to identify a precise time at which the signals change but the precise frequencies are difficult to identify. At the other end of the scale, the 1000 ms window allows the frequencies to be precisely seen but the time between frequency changes is blurred.

Other examples:

${\displaystyle w(t)=exp(\sigma -t^{2})}$

Normally we call ${\displaystyle exp(\sigma -t^{2})}$ a Gaussian function or Gabor function. When we use it, the short-time Fourier transform is called the "Gabor transform".

It can also be explained with reference to the sampling and Nyquist frequency.

Take a window of N samples from an arbitrary real-valued signal at sampling rate f_s . Taking the Fourier transform produces N complex coefficients. Of these coefficients only half are useful (the last N/2 being the complex conjugate of the first N/2 in reverse order, as this is a real valued signal).

These N/2 coefficients represent the frequencies 0 to f_s/2 (Nyquist) and two consecutive coefficients are spaced apart by f_s/N Hz.

To increase the frequency resolution of the window the frequency spacing of the coefficients needs to be reduced. There are only two variables, but decreasing f_s (and keeping N constant) will cause the window size to increase — since there are now fewer samples per unit time. The other alternative is to increase N, but this again causes the window size to increase. So any attempt to increase the frequency resolution causes a larger window size and therefore a reduction in time resolution—and vice versa.

a. Nyquist criterion (avoiding the aliasing effect):

${\displaystyle \Delta _{t}<{\frac {1}{2\Omega }}}$, where ${\displaystyle \Omega }$ is the bandwidth of ${\displaystyle x(\tau )w(t-\tau )}$

a. ${\displaystyle \Delta _{t}\Delta _{f}={\tfrac {1}{N}}}$, where ${\displaystyle N}$ is an integer

b. ${\displaystyle N\geq 2Q+1}$

c. Nyquist criterion (avoiding the aliasing effect):

${\displaystyle \Delta _{t}<{\frac {1}{2\Omega }}}$, ${\displaystyle \Omega }$ is the bandwidth of ${\displaystyle x(\tau )w(t-\tau )}$

${\displaystyle X(n\Delta _{t},m\Delta _{f})=\sum _{p=n-Q}^{n+Q}w((n-p)\Delta _{t})x(p\Delta _{t})e^{-{\frac {2\pi jpm}{N}}}\Delta _{t}}$

${\displaystyle {\text{if we have }}q=p-(n-Q),{\text{ then }}p=(n-Q)+q}$

${\displaystyle X(n\Delta _{t},m\Delta _{f})=\Delta _{t}e^{\frac {2\pi j(Q-n)m}{N}}\sum _{q=0}^{N-1}x_{1}(q)e^{-{\frac {2\pi jqm}{N}}}}$

${\displaystyle {\text{where }}x_{1}(q)={\begin{cases}w((Q-q)\Delta _{t})x((n-Q+q)\Delta _{t})&0\leq q\leq 2Q\\0&2Q<q<N\end{cases}}}$

a. ${\displaystyle \Delta _{t}\Delta _{f}={\tfrac {1}{N}}}$, where ${\displaystyle N}$ is an integer

b. ${\displaystyle N\geq 2Q+1}$

c. Nyquist criterion (avoiding the aliasing effect):

${\displaystyle \Delta _{t}<{\frac {1}{2\Omega }}}$, ${\displaystyle \Omega }$ is the bandwidth of ${\displaystyle x(\tau )w(t-\tau )}$

d. Only for implementing the rectangular-STFT

Rectangular window imposes the constraint

${\displaystyle w((n-p)\Delta _{t})=1}$

Substitution gives:

${\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}}}$

Change of variable n-1 for n:

${\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}}$

Calculate ${\displaystyle X(\min {n}\Delta _{t},m\Delta _{f})}$ by the N-point FFT:

${\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)}}$

where

${\displaystyle x_{1}(q)={\begin{cases}x((n-Q+q)\Delta _{t})&q\leq 2Q\\0&q>2Q\end{cases}}}$

Applying the recursive formula to calculate ${\displaystyle X(n\Delta _{t},m\Delta _{f})}$

${\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}}$

${\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})}}$

so

${\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}}}$

${\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}}}$