Computing filter(b,a,x,zi) using FFTs

Question

I would like to try to compute y=filter(b,a,x,zi) and dy[i]/dx[j] using FFTs rather than in the time domain for possible speedup in a GPU implementation.

I am not sure it's possible, particularly when zi is non-zero. I looked through how scipy.signal.lfilter in scipy and filter in octave are implemented. They are both done directly in the time domain, with scipy using direct form 2 and octave direct form 1 (from looking through code in DLD-FUNCTIONS/filter.cc). I haven't seen anywhere an FFT implementation analogous to fftfilt for FIR filters in MATLAB (i.e. a = [1.]).

I tried doing y = ifft(fft(b) / fft(a) * fft(x)) but this seems to be conceptually wrong. Also, I am not sure how to handle the initial transient zi. Any references, pointer to existing implementation, would be appreciated.

Example code,

import numpy as np
import scipy.signal as sg
import matplotlib.pyplot as plt

# create an IRR lowpass filter
N = 5
b, a = sg.butter(N, .4)
MN = max(len(a), len(b))

# create a random signal to be filtered
T = 100
P = T + MN - 1
x = np.random.randn(T)
zi = np.zeros(MN-1)

# time domain filter
ylf, zo = sg.lfilter(b, a, x, zi=zi)

# frequency domain filter
af = sg.fft(a, P)
bf = sg.fft(b, P)
xf = sg.fft(x, P)
yfft = np.real(sg.ifft(bf/af * xf))[:T]

# error
print np.linalg.norm(yfft - ylf)

# plot, note error is larger at beginning and with larger N
plt.figure(1)
plt.clf()
plt.plot(ylf)
plt.plot(yfft)

Answer 1

You can reduce the error in your existing implementation by replacing P = T + MN - 1 with P = T + 2*MN - 1. This is purely intuitive, but it seems to me that the division of bf and af will require 2*MN terms, due to wraparound.

C.S. Burrus has a pretty terse writeup of how to regard filtering, whether FIR or IIR, in a block oriented way, here. I haven't read it in detail, but I think it gives you the equations you need to implement IIR filtering by convolution, including intermediate states.