6
$\begingroup$

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?

$\endgroup$

1 Answer 1

8
$\begingroup$

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:

enter image description here

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:

enter image description here

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:

enter image description here

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):

enter image description here

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:

enter image description here

enter image description here

"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.

$\endgroup$
3
  • 2
    $\begingroup$ Great explanation of two different interpretations of the STFT (sliding windowed transform vs complex filterbank). I especially appreciate the often-missed fact that the hop size from the sliding window view can be mapped to a downsampling factor in the complex filterbank view. And thank you for the great CWT animation! $\endgroup$
    – Jdip
    Mar 3, 2023 at 17:23
  • 1
    $\begingroup$ @Jdip Well it's good to see someone other than I acknowledge this, I've never seen it elsewhere and keep (re)inventing stuff. Only saw it in a synchrosqueezing paper, and all they say is, "it's convenient for our purposes"... like, ok?? More like "absolutely necessary and for good reasons". $\endgroup$ Mar 4, 2023 at 13:13
  • 1
    $\begingroup$ Hay, I'm in agreement with you about the "Improved STFT". $\endgroup$ Mar 6, 2023 at 20:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.