A quick code tutorial on overlap-save vs single giant IFFT/FFT vs numpy.convolve to achieve linear convolution.

overlapSave.py

"""
Overlap-save approach to fast-convolution.

This script shows how to use

- one long FFT,
- several short FFTs (overlap-save, or overlap-discard),
- and `numpy.convolve`

to achieve the same end: linear convolution. (I don't like overlap-add so that's not included.)

The *one weird trick* I use to remember overlap-save is this: if you take a filter and a
subvector of the same length from inside a longer vector, and do the FFT-based circular
convolution, i.e., if you have a short `h` filter vector and a long `x` data vector and then
```
partial = fft.ifft(fft.fft(x[5:5+h.size]) * np.conj(fft.fft(h)))
```
you'll find that the FIRST element of this result is the same as the element in the true linear
convolution via FFT, that is,
```
full = fft.ifft(fft.fft(x) * np.conj(fft.fft(h, x.size)))
assert full[5] == partial[0]
```

In other words, if you do FFT-circulular-convolution with the shortest data subvector and
the filter, `Nx - Nh + 1 = 1` output sample will be "correct" (equal to the linear convolution
case), and the rest will be invalid (from the linear convolution perspective).

The longer subvectors you use, the more samples of the circular convolution will match the linear
convolution.
"""

import numpy as np
import numpy.fft as fft


def overlapSave(x, h, nfft=None):
    nfft = max(nfft or x.size + h.size - 1, h.size)
    hf = np.conj(fft.fft(h, nfft))

    nvalid = nfft - h.size + 1
    y = np.empty_like(x)
    offset = 0
    chunk = np.zeros(nfft, dtype=x.dtype)
    while offset < x.size:
        thislen = min(offset + nfft, x.size) - offset
        chunk[thislen:] = 0
        chunk[:thislen] = x[offset:offset + thislen]
        chunk = fft.ifft(fft.fft(chunk) * hf)

        thisend = min(offset + nvalid, x.size) - offset
        if np.iscomplexobj(y):
            y[offset:offset + thisend] = chunk[:thisend]
        else:
            y[offset:offset + thisend] = chunk[:thisend].real

        offset += nvalid
    return y


if __name__ == '__main__':
    # Generate some random complex data (works fine for real too)
    x = np.random.randn(111) + 1j * np.random.randn(111)
    h = np.random.randn(7) + 1j * np.random.randn(7)
    nconv = x.size + h.size - 1
    expected = fft.ifft(
        fft.fft(x, nconv) * np.conj(fft.fft(h, nconv)))[:x.size]
    actual = overlapSave(x, h, 16)
    assert np.allclose(actual, expected)

    conv = np.convolve(x, np.conj(h[::-1]))[h.size - 1:]
    assert np.allclose(actual, conv)