Short-Time Fourier Transform

The short-time Fourier transform (STFT), or alternatively short-term Fourier transform, is a Fourier-related transform used to determine the sinusoidal frequency and phase content of local sections of a signal as it changes over time.In practice, the procedure for computing STFTs is to divide a longer time signal into shorter segments of equal length and then compute the Fourier transform separately on each shorter segment. This reveals the Fourier spectrum on each shorter segment. One then usually plots the changing spectra as a function of time.

Continuous-time STFT

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 as the window is slid along the time axis, resulting in a two-dimensional representation of the signal. Mathematically, this is written as:

STFT[x(t)](τ,ω)≡X(τ,ω)=∫∞−∞x(t)w(t−τ)e−jωtdt

where w(t) is the window function, commonly a Hann window or Gaussian window centered around zero, and x(t) is the signal to be transformed.(Note the difference between w and ω.) X(τ,ω) is essentially the Fourier Transform of x(t)w(t−τ), 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, τ, and frequency axis, ω, to suppress any jump discontinuity of the phase result of the STFT. The time index τ is normally considered to be “slow” time and usually not expressed in as high resolution as time t.

Discrete-time STFT

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:

STFT[x[n]](m,ω)≡X(m,ω)=∑n=−∞∞x[n]w[n−m]e−jωn

ikewise, 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 of the function:

spectrogram[x(t)](τ,ω)≡|X(τ,ω)|2

Inverse STFT

The STFT is invertible, that is, the original signal can be recovered from the transform by the Inverse STFT. The most widely accepted way of inverting the STFT is by using the overlap-add (OLA) method, which also allows for modifications to the STFT complex spectrum. This makes for a versatile signal processing method.

Continuous-time STFT

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

x(t)=∫∞−∞[12π∫∞−∞X(τ,ω)e+jωtdω]dτ

x(t)w(t−τ)=12π∫∞−inftyX(τ,ω)e+jωtdω

PRACTIACL COMPUTATION OF THE STFT

While the definition of the STFT is useful for theoretical work, it is not really a specification of a practical STFT. In practice, the STFT is computed as a succession of FFTs of windowed data frames, where the window slides'' orhops” forward through time. We now derive such an implementation of the STFT from its mathematical definition.

Xm(ω)=∑n=−∞∞x(n+mR)w(n)e−jω(n+mR)

If the window w(n) has the Constant OverLap_Add(COLA) property at hop-size R, i.e., if

∑m=−∞∞w(n−mR)=1,∀n∈Z(w∈COLA(R))

so,

Xm(ω)=e−jωmR∑n=−∞∞x(n+mR)w(n)e−jωn

Xm(ω)=e−jωmRDTFTω(SHIFT−mR(x)⋅w)

In this form, the data centered about time mR are translated to time 0, multiplied by the (let’s assume zero-phase) window w , and then the DTFT is performed. Since the nonzero portion of the windowed data is centered on time zero, the DTFT can be replaced by the DFT (or FFT). This effectively samples the DTFT in frequency. This sampling will not cause (time) aliasing if the number of samples around the unit circle is greater than the width (in samples) of the time interval including all nonzero datapoints. In other words, sampling the frequency axis is information-preserving when the signal is properly time limited. Let M denote the window length (typically an odd number) and N≥M be the DFT length (typically a power of 2). Then at ω=ωk=2πk/N, k=0,1,2,...,N−1, and using the fact that the window w(n) is time-limited to less than N samples centered about time zero, yields

Xm(ωk)=e−jωkmR∑n=−N/2N/2−1x(n+mR)w(n)e−jωkn

Xm(ωk)=e−jωkmRDFTN,ωk(SHIFT−mR(x)⋅w)

Since indexing in the DFT is modulo N , the sum over n can be “rotated” to a sum from 0 to N−1 as is conventionally implemented for the DFT. In practice, this means that the right half of the windowed data frame goes at the beginning of the FFT input buffer, and the left half of the windowed frame goes at the end, with zero-padding in the middle.

STFT IN MATLAB

The following matlab segment illustrates the above processing steps:

    Xtwz = zeros(N,nframes); % pre-allocate STFT output array    M = length(w);           % M = window length, N = FFT length    zp = zeros(N-M,1);       % zero padding (to be inserted)    xoff = 0;                % current offset in input signal x    Mo2 = (M-1)/2;           % Assume M odd for simplicity here    for m=1:nframes    xt = x(xoff+1:xoff+M); % extract frame of input data    xtw = w .* xt;         % apply window to current frame    xtwz = [xtw(Mo2+1:M); zp; xtw(1:Mo2)]; % windowed, zero padded    Xtwz(:,m) = fft(xtwz); % STFT for frame m    xoff = xoff + R;       % advance in-pointer by hop-size R    end

This article comes from https://en.wikipedia.org/wiki/Short-time_Fourier_transform and https://www.dsprelated.com/freebooks/sasp/Short_Time_Fourier_Transform.html