Fast convolution

Two routes to one convolution, in a race. Watch the solid FFT curve pass the dashed direct one as the length grows.

Convolution through the FFT

On the page Circular vs linear convolution (13.4) you saw that multiplying two DFTs convolves the two sequences circularly, and that padding both with zeros turns the circle into a line. This page asks a practical question: when is that route faster than flip and slide, and how do you use it on a signal that never ends?

Here is the recipe. Let $x$ have $N_x$ samples and $h$ have $N_h$. The linear convolution has $N_x+N_h-1$ samples (see Discrete convolution, 5.2). Choose the FFT size $N_\text{fft}$: the smallest power of two that is at least $N_x+N_h-1$. Then do five things.

1. Zero-pad $x$ and $h$ to $N_\text{fft}$ samples.
2. Take the FFT of each.
3. Multiply the two spectra, bin by bin.
4. Take the inverse FFT of the product.
5. Keep the first $N_x+N_h-1$ samples.

I call this FFT convolution. In symbols, with both inputs already padded:

Nothing about it is approximate. It gives the same samples as flip and slide, up to rounding in the arithmetic.

When the FFT route wins

To compare the two routes I count real multiplies, because an FFT works on complex numbers and a complex multiply costs 4 real ones, as on the page Complex numbers for signals (3.3). I use no tricks for real inputs. The two-for-one trick of More FFT algorithms (14.2) would roughly halve the FFT part, and I leave it out so the count stays simple.

Take two signals of the same length $N$. Flip and slide multiplies every sample of one by every sample of the other: $N^2$ real multiplies. The FFT route needs three FFTs (two forward, one inverse) of size $N_\text{fft}$, and $N_\text{fft}$ complex products. From The FFT (14.1), one FFT costs $\tfrac{N_\text{fft}}{2}\log_2N_\text{fft}$ complex multiplies. So

The FFT route is like a motorway: slower to get onto, faster for a long trip. For short signals the fixed cost of three FFTs is more than the $N^2$ you were trying to avoid. For long signals $N^2$ grows much faster than $N\log_2N$.

The race at the top of this page runs the formula. At its start, two 4-sample signals, the FFT size is 8, and $4(3\cdot4\cdot3+8)=176$ multiplies against 16 by hand. At $N=256$ it is 65 536 multiplies directly and 29 696 by FFT. At $N=4096$ it is 16.8 million against 671 744, 25 times fewer.

The FFT curve is a staircase, not a smooth line. The size $N_\text{fft}$ has to be a power of two, so it jumps each time $2N-1$ passes one. The jump is why the crossover is not a single point. Here are the numbers around it.

From $N=116$ to $128$ the FFT route wins, because $N_\text{fft}=256$ is already enough. At $N=129$ the length $2N-1=257$ needs 512, the cost more than doubles (13 312 to 29 696), and direct wins again until $N=173$. From 173 onward the FFT route is cheaper at every length up to 4096.

After the sweep you can choose $N$ yourself in that picture. Try 128, then 129, and watch the cheaper trace swap twice.

Overlap-add: a signal that never stops

A microphone never stops, so you cannot wait for the whole signal before you start. Instead, cut $x$ into blocks of $N_\text{blk}$ samples each. Convolution is linear, so by the superposition of The impulse response (5.1) the output is the sum of what each block produces alone.

Each block’s output is $N_\text{blk}+N_h-1$ samples long. That is $N_h-1$ more than the block, so the end of each output, its tail, spills into the time of the next block. The fix is to place each block’s output at the block’s own start and to add wherever tails overlap. This is overlap-add. Think of tiling a floor with tiles that overlap at the edges: the overlapped strips get painted twice, and you add the two coats.

Watch the two hatched zones, where one block’s tail lands on the start of the next: only there do two outputs add.

The input $x$ has three blocks of 8: 2, 1, 3, 0, −1, 2, 1, 1; then 0, 2, −2, 1, 3, 1, 0, −1; then 1, 1, 2, 0, −1, 0, 2, 1. The filter is $h=1,1,1$. Each block of 8 convolved with $h$ gives $8+3-1=10$ samples, so a 16-point FFT is enough.

The block outputs, each computed with a 16-point FFT, are:

• block 1, at $n=0$ to 9: $2,3,6,4,2,1,2,4,2,1$;
• block 2, at $n=8$ to 17: $0,2,0,1,2,5,4,0,-1,-1$;
• block 3, at $n=16$ to 25: $1,2,4,3,1,-1,1,3,3,1$.

Block 2 starts at $n=8$, on top of block 1’s tail: at $n=8$ and 9 the two add, $2+0=2$ and $1+2=3$. At $n=16$ and $17$ the tails of block 2 meet the start of block 3: $-1+1=0$ and $-1+2=1$. Only those four samples are sums; every other sample belongs to one block alone.

Added up, the 26 samples of $y$ are 2, 3, 6, 4, 2, 1, 2, 4, 2, 3, 0, 1, 2, 5, 4, 0, 0, 1, 4, 3, 1, −1, 1, 3, 3, 1, and they equal the output of flip and slide sample for sample. Each block needed only a 16-point FFT, however long $x$ is.

Overlap-save: keep what did not wrap

There is a second way to cut a signal, and it works the other way round. Let the FFT wrap, then throw away what wrapped. This is overlap-save.

Take overlapping segments of $N_\text{fft}$ input samples. Each segment starts $N_\text{fft}-N_h+1$ samples after the last one, so it repeats the end of the one before. Convolve each segment with $h$ circularly, as on page 13.4. The circular convolution wraps the tail of the segment back onto its first $N_h-1$ outputs, so those are wrong. The remaining $N_\text{fft}-N_h+1$ outputs have nothing wrapped onto them and are correct. Keep those and drop the rest.

The figure uses the same $x$ and $h$ as before. The segment holds two zeros of history, then $x[0]$ to $x[5]$. Of the eight circular outputs, the first two ($1,2$) are wrapped. The other six equal $y[0]$ to $y[5]$. Each segment yields $8-3+1=6$ good samples, so the next one starts 6 samples later. The two methods differ in where the repeated work sits: overlap-add overlaps the outputs and adds, overlap-save overlaps the inputs and discards.

Bigger blocks: cheaper, but later

Overlap-add leaves one choice: the block length $N_\text{blk}$. To see what it costs, take a real job, a reverb. The echo of a room is its impulse response, as on page 5.1. Record it once, and any dry sound can be placed in that room by convolution with it.

A hall with a 1-second echo, sampled at $f_s=16\,000$ Hz, has $N_h=16\,000$ taps. Directly, that is 16 000 multiplies for every output sample. With overlap-add each block costs two FFTs (the FFT of $h$ is done once and kept) and $N_\text{fft}$ products, where $N_\text{fft}$ is the first power of two at least $N_\text{blk}+N_h-1$. Spread over the $N_\text{blk}$ outputs of the block, the cost per output sample is

There is a price. The first output of a block can only appear after the whole block has arrived, so every output is late by about $N_\text{blk}/f_s$. This waiting time is the latency. It is like a bus that waits to fill: cheaper per passenger, later to leave.

Watch the cost fall as the blocks double, while the delay doubles with them, until past 16 384 the cost climbs again.

Blocks of 64 give only 4 ms of delay, but 15 360 multiplies per output, about as many as direct. At 512 the cost rises a little, from 3840 to 4096, because the FFT jumps to 32 768.

Blocks of 16 384 cost 128 multiplies per output, 125 times fewer than direct, but each block waits 1024 ms to fill. This is the cheapest block length of the sweep.

Past that point the FFT grows faster than the block, and the cost climbs again: 144 at 65 536. The choice is delay against cost. The whole table, from the same formula:

Press “Hear the room” in the picture. It plays a 0.3 s dry pluck, a 440 Hz tone with a 60 ms decay, and then the same pluck convolved with a synthetic 1 s room. The room is seeded white noise times $e^{-6.91t/0.8}$, which has decayed to $-60$ dB at 0.8 s. The output is the same for every block length. It sounds the same; only the wait changes.

Real-time reverbs split $h$ into pieces too. This is partitioned convolution: small blocks for the start of $h$ give low delay, and big blocks for its tail give low cost.

Worked example

Use the reverb above with $N_\text{blk}=1024$. The first power of two at least $1024+16\,000-1=17\,023$ is $N_\text{fft}=32\,768$. The cost per output is $4(32768\cdot15+32768)/1024=2048$ real multiplies, against 16 000 directly: 7.8 times fewer. The delay is $1024/16\,000=0.064$ s, which is 64 ms.

Is that a good choice? If 64 ms of delay is acceptable, 7.8 times fewer multiplies is a good return. If you need under 10 ms, blocks of 64 give 4.0 ms but save nothing. That tension is why real systems partition $h$.

Where you’ll meet this

Long FIR filters are applied by FFT in this way, as on the pages of Chapter 19. The short-time Fourier transform of Chapter 15 also cuts a signal into overlapping frames (15.5 and 15.6). Multirate filter banks use the same blocking (Chapter 22).

Reference card