Equivalence between "windowed Fourier transform" and STFT as convolutions/filtering

Question

I've heard, that "windowed Fourier transform" is but one perspective on STFT, and that STFT is fundamentally convolutions of windowed complex sinusoids with the input, i.e. bandpass filtering. Yet the standard implementation and perspective on STFT is, take a segment of input x_seg = x[start:end], window it x_seg * window, and take its FFT (DFT).

Is that really equivalent to conv(x, window * cisoid)? How so?

Answer 1

This answer isn't fully developed (but is correct), I may update in the future.

Indeed, STFT is strided convolutions, and there's many practical caveats, discussed in this answer.

Standard vs improved STFT

Standard:

$$ \text{STFT}_x(\tau, \omega) = \int_{-\infty}^{\infty} x(t) w(t - \tau) e^{-j\omega t} dt $$

It's what most libraries (including MATLAB, scipy, librosa) implement. In short, its complex-valued result is suboptimal, or even flawed - it lacks some important properties of time-frequency transforms, e.g. one-integral inverse and zero phase (even if window is zero phase), making it ill-suited for many kinds of operations on complex values - though the modulus (spectrogram) is unaffected.

The better variant, implemented in ssqueezepy, writes

$$ \text{STFT}_x^\text{mod}(\tau, \omega) = \int_{-\infty}^{\infty} x(t) w(t - \tau) e^{-j\omega (t - \tau)} dt $$

I'm unsure whether MATLAB is really implementing "standard"; see "Demos" below.

---

STFT as convolutions

The simplicity of STFT is obfuscated in the "windowed Fourier" implementation. Now let's visualize:

Usually, we take DFT (FFT) or DTFT of this and call this "the frequency response" and do lobe analysis. But that's wrong! (Before I get pitchforked, see "caveats" below.) The true operators of STFT are convolutions, hence the true frequency response is that of its filters padded to the convolution length. Below is that, on the $\texttt{STFT}^{\texttt{mod}}$ variant:

That's the real part, with imaginary zero: zero-phase, as is ideal. For complex input, the filters would also completely tile the negative frequencies. Note, despite center frequencies being positive-only, there's still leaks into negatives - hence loss of analyticity. Brief detour:

Convolution <=> Windowed Fourier equivalence

"Windowed Fourier transform" takes DFT of a windowed segment of x. Output of DFT at e.g. index 3, is the inner product of a complex sinusoid at DFT frequency 3, with input to DFT, which here is the windowed segment. I.e., it's sum(x_seg * window * cisoid(f=3)), and that's assigned to STFT[3, seg_idx]. More concretely,

• Let y = fft(x), and N=len(x). y[3] is equal to sum(x * np.exp(-1j * 2*np.pi * 3 * np.arange(N) / N)) - or, sum(x * cisoid(f=3)).
• Hence, fft(x * window) just replaces x with x * window above: sum(x * window * cisoid(f=3))

This is repeated for all hops and f to produce the full STFT. Below are window * cisoid(f) for two f:

Above is exactly describing convolution of x with window * cisoid(f), where len(x_seg) < len(x) is permitted since window is time-limited. Yet, for reasons, it's a flawed implementation of said convolution - the modulus is correct, but phase messy; separate topic. We assume $\texttt{STFT}^{\texttt{mod}}$ throughout.

STFT filterbank (cont'd)

out[4, :] is hence computed by convolving with filter4 = window * cisoid(f=4). Multiplication in time is convolution in frequency - hence, the frequency response of filter4 is the frequency response of window convolved with frequency response of cisoid(f=4) - latter is just a unit impulse, and the result is shifting the frequency response of window to f=4. Repeat for all f, and we get the fbank_f0 plot.

Now, let's check n_fft < N:

With all said so far, it makes sense that if N / n_fft is integer, then fbank_f1 is a subset of fbank_f0; you'll find in code, np.allclose(fbank_f1, fbank_f0[::4]) == True.

Now, for N / n_fft that's a fraction, superimposed with above filterbank (orange):

No longer a subset, but some frequencies may still coincide.

Demos

Let $x(t)$ be sum of pure sines of frequencies $4, 20, 40, 60$, with $N=128$ samples. We have:

"Standard" introduces behavioral dependence on n_fft while "improved" resamples at higher resolution as intended; "standard" indicates modulations along frequency for a pure tone, which is nonsense. These have reaching implications which amount to largely invalidating the 2D time-frequency structure of the complex-valued STFT.

MATLAB's stft, from my limited experimentation, has yet additional distortions (ignore edge effects):

This shouldn't be possible in either variant, "standard" or "improved", I've no clue where it's coming from. Likely, it's still "standard", but with extra steps that don't fit the standard's mathematical formulation - I welcome others' input here. (Note, here it looks closer to "improved", but other examples I've tried rule it out.)

t = 0:1/128:(1 - 1/128);
x = cos(2*pi*30*t);
xp = cat(2, zeros(1, 10), x, zeros(1, 9));
Sxx = real(stft(xp, Window=hamming(20), OverlapLength=19, FFTLength=128));

imagesc(Sxx(1:65, :)); figure
plot(Sxx(30, 45:70)); hold on; plot(Sxx(32, 45:70)); plot(Sxx(34, 45:70));

Caveats / observations

1. Above when I say convolution, I actually mean cross-correlation.
2. n_fft=N corresponds to (actually) convolution of fft(window) with fft(x) for the overall frequency response of STFT
3. Even if original window is zero-phase, it may not be so in the true response. In fact, in majority of cases, it's not - though it's often close. Zero-phase requires DFT-symmetry - and, for example, zero-padding kills that for even-length DFT-symmetric windows with no zeros in window. But zero-padding isn't causing the problem, only exposing it.
4. "But that's wrong!" - it's correct columnwise. But almost any time we speak of filtering, we speak of conv(x, filter), and fft(window) doesn't reflect that. More importantly, zero-padding in one domain is interpolation in another - meaning, the magnitude response of fft(window) remains accurate; it's phase that most likely isn't.
5. hop_size = 4 is just hop_size = 1 with convolution stride 4, or equivalently, subsampling STFT(x, hop_size=1) by a factor of 4, like so: Sx4 = Sx1[:, ::4]. The effects on frequency are understood via subsampling in time <=> folding in frequency.
6. STFT pads and unpads, and unpadding is aliasing - in this case exact frequency representation is (likely) impossible, but in most cases it's very close.
7. scipy's sym=False for odd lengths is incorrect

CWT animation

1. Slide wavelet psi_i by some time shift tau
2. Compute similarity (inner product, sum(x * psi))
3. Repeat for all tau and psi_i

STFT animation

It's exactly "CWT animation", where psi_i = fft(window * cisoid_i) (i is index of a given frequency, and window is full length), except

1. Filter bandwidth is constant
2. Filters are spaced uniformly (linearly distributed)

There's also zero-mean enforcement in CWT, but that's negligible for our purposes.

What's in code

I provide a limited row-wise implementation of STFT (mod) for exploration purposes, that also returns the exact filterbank used to compute STFT. Main caveats:

1. The filterbanks aren't actually zero-phase, though it's very close - per technicalities, I've not had time to account for it, though you could. Easiest workaround is scipy.signal.windows with input length N.
2. There's no padding, but it's equivalently 'wrap' (periodic) padding. This bypasses the unpadding caveat and represents the filters perfectly. It can also be worked around by passing padded x as input, and unpadding the output.

How to translate between standard and "mod" implems?

It's a simple conversion - code at ssqueezepy, discussion here (also see comments above the one linked).

Code

Available at Github.