# How to implement a 3rd order Butterworth IIR allpass filter with perfectly flat frequency response?

I'm working on an audio plugin which splits audio into frequency bands for individual processing. The filters must have zero latency, so no linear phase filters.

I would like to have perfect reconstruction in the magnitude spectrum, so that when you add together the lowpassed and highpassed signals, you get a perfectly flat frequency response, although gradual phase shifts are ok. Like in crossover filters in speakers.

Working from the analog prototype the math works out perfectly for 1st to 4th order filters. As long as you chose the right gain at the cutoff frequency, and sometimes phase flip one of the signals, you get a perfectly flat magnitude spectrum in the reconstructed signal. It probably works out perfectly for higher orders as well, but I haven't got around to check.

For my digital implementations, it also works perfectly for 1st (6dB/octave), 2nd (12dB/octave) and 4th order (24dB/octave) filters. But for 3rd order filters, the magnitude response is almost perfectly flat up to about 1k crossover frequency, and gets small dips and valleys around the crossover frequency as you increase it towards the nyquist frequency at 22.05 kHz.

From the analog prototype math, I should get perfectly flat reconstruction by cascading a 2nd order filter with a 0dB gain at f0, with a 1st order filter with -3dB gain at the crossover. This should give a Butterworth 3rd order filter, and when adding the hi- and lopassed signals, you should get a perfectly flat response. I have tested and verified this in Python using the ideal analog responses. In the digital implementation I am using the RBJ cookbook filters for the 2nd order stage, setting Q to 1.

For the single pole filters I am using these coeffs:

X = math.exp(-2.0*pi*f0/fs)


Lo-pass:

a0 = (1-X)/2
a1 = (1-X)/2
b1 = X


Hi-pass:

a0 = (1+X)/2
a1 = -(1+X)/2
b1 = X


The lo-pass coeffs give a more correct phase response compared to the more common a0 = 1-X and b1 = X which you find everywhere on the internet.

None of the filters give ideal response towards the nyquist frequency by themselves, but when added together it just magically works to give a perfectly flat response. But only for 1st, 2nd and 4th order.

Here are the responses when adding the 3rd order lopassed and highpassed signals in Python:

Any ideas?

OK, I found the answer. It turns out that the exp-based coefficients in the 1st order filter, referenced in my question, which is taught in engineering textbooks, tutorials, and pretty much everywhere else, gives a -3dB point which is slightly off. It becomes more off the closer you push Fc towards the nyquist frequency.

So, I sat down and derived coefficients myself from the BLT, and now it all sums up to a perfectly flat frequency response.

Here are some code snippets in Python. I haven't tried to simplify the math yet, which may be possible, anyway it works now:

def single_pole_LP(data, f0, fs):
W = math.tan(pi*f0/fs)
a0 = W/(1+W)
a1 = a0
b1 = -(W-1)/(1+W)

result = np.zeros(len(data)+1)
data = np.concatenate((np.zeros(1), data))
for i in range(1,len(data)):
result[i] = a0*data[i]+a1*data[i-1]+b1*result[i-1]
return result[1:]

def single_pole_HP(data, f0, fs):
return data - single_pole_LP(data, f0, fs)

def Butterworth_3rd_LP(data, f0, fs):
temp = rbj_LP(data, f0, fs, 1.0)
return single_pole_LP(temp, f0, fs)

def Butterworth_3rd_HP(data, f0, fs):
temp = rbj_HP(data, f0, fs, 1.0)
return single_pole_HP(temp, f0, fs)

def allPass_1st(data, f0, fs):
lp = single_pole_LP(data, f0, fs)
hp = single_pole_HP(data, f0, fs)
return( lp + hp )

def allPass_2nd(data, f0, fs):
lp = rbj_LP(data, f0, fs, 0.5)
hp = rbj_HP(data, f0, fs, 0.5)
return( lp - hp )

def allPass_3rd(data, f0, fs):
lp = Butterworth_3rd_LP(data, f0, fs)
hp = Butterworth_3rd_HP(data, f0, fs)
return( lp + hp )

def allPass_4th(data, f0, fs):
q = 1.0/math.sqrt(2.0)
lp = rbj_LP(data, f0, fs, q)
lp = rbj_LP(lp, f0, fs, q)
hp = rbj_HP(data, f0, fs, q)
hp = rbj_HP(hp, f0, fs, q)
return( lp + hp )


The rjb_* algorithms can be found here. The fourth parameter is the Q.

Here are the new frequency responses for the 3rd order, with crossover frequencies at 1kHz, 5kHz and 10kHz

EDIT: I have now found the coefficients for the single pole high pass filter. they are:

a0 = 1/(1+W)
a1 = -a0


b1 stays the same. The results are exactly identical to subracting the lo-passed signal from the original signal, except for roundoff errors, but if you only need the hi-passed signal you might save a few nanoseconds of performance... :)

EDIT 2: Different authors use different conventions for naming the coefficients, and for the signs of the coefficients. The convention I'm using for the 1st order LPF is taken from the DSP guide here, where the recursion algorithm is:

y[n] = a0*x[n]+a1*x[n-1]+b1*y[n-1]


In the RBJ cookbook, the convention is that the recursion algorithm is:

y[n] = b0*x[n]+b1*x[n-1]-a1*y[n-1]


With this convention, the LPF coefficients become:

W = math.tan(pi*f0/fs)


LPF:

b0 = W/(1+W)
b1 = b0
a1 = (W-1)/(1+W)


HPF:

b0 = 1/(1+W)
b1 = -b0
a1 = (W-1)/(1+W)


In a comment is suggested another set of coefficients using sine and cosine instead of tan. This set of coefficients give exactly the same results. The tan function should give infinite if you set f0 to the nyquist frequency fs/2, which might indicate the sine/cosine based coeffs are better, however testing shows that my set of coeffs work as expected with f0 set to nyquist (maybe due to inaccuracies in the pi constant), while the other suggested coeffs produce an overflow. It probably means none of them are safe with f0 == nyquist, so you should set an upper limit on f0 just below nyquist anyway.

• Butterworth ist order LPF: H(s) = 1/(s+1): b0 = sin(w0) b1 = sin(w0) a0 = cos(w0) + sin(w0) + 1 a1 = sin(w0) - cos(w0) - 1 and 1st order HPF: H(s) = s/(s+1): b0 = cos(w0) + 1 b1 = -(cos(w0) + 1) a0 = cos(w0) + sin(w0) + 1 a1 = sin(w0) - cos(w0) - 1 Commented Dec 3, 2022 at 15:43
• Thank you. Those coefficients give me completely wrong results, but it doesn't matter, as my own coefficients are correct, give exact results and are (slightly) simpler to implement. Commented Dec 4, 2022 at 10:28
• PS: I am assuming that w0 is 2.0*pi*f0/fs (where f0 is cutoff and fs is the sampling rate in the same unit), and that you have switched a's and b's compared to what I am using in my code above (otherwise we just get a flat response with a dB boost). I am also assuming all coefficients must be scaled by 1/a0 (your a0), otherwise we just get an overflow. Then the filter just gives a lowpass response with a completely wrong cutoff frequency. But as stated, it doesn't matter. The problem is solved, and I am getting the correct results with my own coeffs. Commented Dec 4, 2022 at 10:47
• PS: The convention I am using for naming a's and b's is taken from dspguide.com, where y[n] = a0*x[n]+a1*x[n-1]+b1*y[n-1] Commented Dec 4, 2022 at 10:57
• OK, I got it to work now, by negating b1. There are different conventions for using negative b's, and adding them in the recursion algorithm, or positive b's and subtracting them. AND for calling them either b or a. Now my and your sets of coeffs give exactly the same response, except for miniscule rounding errors. Commented Dec 4, 2022 at 11:14