10.5 Inverse STFT and spectral processing

import matplotlib
if not hasattr(matplotlib.RcParams, "_get"):
    matplotlib.RcParams._get = dict.get

10.5 Inverse STFT and spectral processing#

We have been computing the STFT; now let us invert it. Is the STFT invertible? We already know the DFT is, since \(x = \texttt{IDFT}(\texttt{DFT}(x))\). So under a rectangular window at 0% overlap, where the frames tile the signal exactly, the STFT is invertible too: applying the inverse DFT to each frame recovers that frame, and overlap-add stitches the frames back together,

\[\texttt{ISTFT}(\texttt{STFT}(x)) = x.\]

Intuitively, the exact invertibility of the DFT implies that the STFT does not change the reconstruction properties of standard frame-based processing. Accordingly, for other windows and overlaps, the same COLA condition from before guarantees perfect reconstruction: as long as the windows overlap-add to a constant, the inverse DFTs stitch the frames back into the original signal, up to a constant gain we divide out. A runnable STFT and inverse STFT are in code/stft.py.

Spectral processing#

The invertibility of the STFT unlocks a whole family of effects. We can transform a sound into the time-frequency domain, edit the spectral coefficients however we like, and transform back, a technique called spectral processing. Now that we have analysis and synthesis in hand, the whole pipeline is a single frame-based flow with an editing step in the middle:

A left-to-right block diagram: the input signal is split into windowed frames, each frame is sent through a DFT, the resulting spectra can be edited, then each is sent through an inverse DFT, and finally the frames are overlap-added back into an output signal. The first half is labeled analysis (STFT) and the second half synthesis (ISTFT).

Fig. 73 The full STFT pipeline. Analysis (the STFT) frames the signal and takes the DFT of each frame; synthesis (the inverse STFT) takes the inverse DFT of each frame and overlap-adds the results. Editing the spectra in between is spectral processing.#

Three quick examples. First, a brick-wall low-pass filter: zero out every bin above a cutoff in every frame, muting the high end. Second, phase randomization: keep each frame’s magnitudes but replace its phases with random values, smearing the sound’s sharp transients into a wash. Third, cross-synthesis: reshape the trio with a human voice so that the trio appears to “speak,” a “talking instrument” effect. The code shows exactly how each is computed.

Brick-wall low-pass (bins above 1 kHz zeroed)

Phase randomized (transients smeared)

The voice used for cross-synthesis

Cross-synthesis (trio made to “speak” by the voice)

Three spectral-processing effects, all computed by editing the STFT and inverting it. For the cross-synthesis, we flatten the trio’s own spectral shape and impose the formants of a spoken clip (Alvin Lucier’s I Am Sitting in a Room) in their place, so the trio keeps its own pitch and rhythm but takes on the shape of the speech.

There is an enormous space of effects to explore here. Try inventing your own by editing the STFT directly:

# hide
import numpy as np
import pyquist as pq


def hann(n):
    return 0.5 * (1 - np.cos(2 * np.pi * np.arange(n) / n))


def stft(x, N_H, N_F, window):
    return np.array([np.fft.rfft(x[s:s + N_F] * window)
                     for s in range(0, len(x) - N_F + 1, N_H)])


def istft(S, N_H, N_F, window):
    out = np.zeros(N_H * (S.shape[0] - 1) + N_F)
    wsum = np.zeros_like(out)
    for k in range(S.shape[0]):
        out[k * N_H:k * N_H + N_F] += np.fft.irfft(S[k], N_F) * window
        wsum[k * N_H:k * N_H + N_F] += window ** 2
    return out / np.maximum(wsum, 1e-8)

# Spectral processing: transform to the time-frequency domain with the STFT,
# EDIT the complex coefficients, and transform back. Try your own edits!
audio = pq.Audio.from_file("./assets/audio-trio.wav")
x = np.asarray(audio.samples).reshape(-1)
sr = audio.sample_rate
N_F, N_H = 2048, 512
w = hann(N_F)

S = stft(x, N_H, N_F, w)          # S is complex, shape (num_frames, num_bins)

# --- edit here --- (this example keeps the magnitudes but randomizes the phases)
magnitude = np.abs(S)
phase = np.random.uniform(-np.pi, np.pi, S.shape)
S = magnitude * np.exp(1j * phase)
# -----------------

y = istft(S, N_H, N_F, w)
pq.play(pq.Audio(y, sr))

The phase vocoder#

We end our discussion of the STFT with a famous spectral-processing algorithm, the phase vocoder, which performs high-quality time stretching without pitch shifting.

Granular synthesis already gave us pitch-preserving time stretching. But we can also frame time stretching as a spectral-processing operation: to slow a sound to half speed, we want an output STFT with twice as many frames, so we simply interpolate between the input frames. Let’s define \(X[j] = \texttt{STFT}_j(x)\), i.e., the DFT of the \(j\)-th frame of \(x\). For output frame \(i\), we blend input frames \(j = \lfloor i/2 \rfloor\) and \(j+1\):

\[Y[i] = (1 - a)\, X[j] + a\, X[j+1], \qquad a = \tfrac{i}{2} - j.\]

This is not unlike the interpolation operation in wavetable synthesis except applied to complex-valued frames instead of wavetable samples. This sounds reasonable, but it has a subtle flaw involving phase. Each STFT bin has a phase, and when we interpolate we are implicitly assuming we know how that phase advances from one frame to the next. But the phase is only known modulo \(2\pi\): if a bin’s phase reads \(\pi/4\) in one frame and \(5\pi/4\) in the next, did it advance by \(\pi\), or by \(3\pi\), or by \(5\pi\)? The sliding-window nature of the STFT makes this ambiguous, and naive interpolation between ambiguous phases produces a smeared, “phasey” artifact.

Two unit circles side by side. The left shows a phasor at angle pi over four; the right shows a phasor at angle five pi over four, half a turn further around. A caption notes the phase could have advanced by pi, or three pi, or five pi.

Fig. 74 The trouble with phase. A bin’s phase jumps from \(\pi/4\) to \(5\pi/4\) between frames, but the true advance could be \(\pi\), \(3\pi\), \(5\pi\), or any of infinitely many possibilities. The STFT alone cannot disambiguate them.#

The phase vocoder resolves this by predicting how the phase should evolve. Bin \(k\) corresponds to an angular frequency of \(\omega_k\) radians per second (as we defined it in Chapter 8), or \(\omega_k / f_s\) radians per sample. So over a single hop of \(N_H\) samples its phase should advance by an expected amount of \(\frac{\omega_k}{f_s} \cdot N_H\) radians. The algorithm compares this expected advance to the observed advance (the actual phase difference between two consecutive frames) and resolves the \(2\pi\) ambiguity by picking whichever multiple lands nearest the expectation. Accumulating these corrected advances frame by frame builds a clean, continuous phase for the output. The details are beyond our scope, but the result is time stretching that exceeds the quality of granular synthesis:

Original

Half speed (pitch preserved)

Double speed (pitch preserved)

Pitched down an octave (phase vocoder + resampling)

The phase vocoder stretches time while holding pitch constant, and combined with resampling it gives independent control over both.

The phase vocoder is the high-quality time-stretching algorithm behind the “playback speed” controls you use every day: the 1.5x and 2x buttons on video and podcast platforms, and the speed sliders in audio software, all rely on it (or a close relative) to speed up or slow down without turning every voice into a chipmunk.