# Design a linear-phase FIR filter approximating the magnitude of a given IIR filter

I have a biquad IIR filter from which I want to get a linear-phase FIR. I saw that related question but the OP does not care much about phase.

I have troubles getting the exact same magnitude curve when extracting an FIR from my IIR.

So far, my method has been:

1. Apply my IIR filter on a 1024 samples buffer which contains nothing but a Dirac as first sample.
2. Shift the impulse response obtained in step 1 and make it symmetric (to get phase linearity)

If I stop after step 1, I get the exact same magnitude curve for my IIR and my FIR but I also get the exact same phase curve (which is non linear, and hence not interesting).
If I stop after step 2, I get an linear phase but not the exact same magnitude.

On the following picture, IIR's and FIR's (after step 2) transfer functions are (respectively) in red and blue: What have I been doing wrong?

• From the image it seems that you are interested in a constant/zero phase filter instead of linear phase. You could obtain this by filtering a signal forwards and backwards in time by the same filter. However this operation is noncasual, so can't be done in real time. You can delay this filter in time to make it casual, but this makes the filter linear in phase. – fibonatic Sep 2 '17 at 14:24
• Is having a mixed bag allowed? If so, use your IIR with an allpass FIR, for equalization. – a concerned citizen Sep 3 '17 at 7:12

What you do in step 1 is simply truncate the infinite impulse response to approximate it by an FIR filter. If you use sufficiently many filter taps, the approximation becomes arbitrarily accurate. This means that the resulting FIR filter approximates the magnitude and the phase characteristic of the original IIR filter. So with this approach the phase will never become linear.

Making the impulse response symmetric to obtain phase linearity, as you do in step 2, will of course change the magnitude response.

What you should do is use the magnitude of the IIR filter as a desired response in a (linear-phase) FIR filter design routine. In that case you will get an FIR filter with an exactly linear phase and with a certain magnitude approximation error. That magnitude error can be made sufficiently small by choosing an appropriate filter order. The simplest approach is probably to use a least squares approximation, which just involves solving a system of linear equations.

Example: I use a peaking EQ filter as the IIR prototype. The coefficients are (b are the numerator coefficients, a are the denominator coefficients):

b = [1.2223e+00, 0, 7.7775e-01];
a = [1.1250e+00, 0, 8.7502e-01];


You can use the magnitude of the IIR filter's frequency response and combine it with a linear phase to obtain the desired response for the FIR filter design routine (N is the filter length). The code is Matlab/Octave syntax:

[H,w] = freqz(b,a,256);
N = 61;
D = abs(H).*exp(-1i*w*(N-1)/2);


You can use a least squares FIR filter design routine called lslevin.m, which you can find here.

h = lslevin(N,w,D,ones(length(w),1));
Hh = freqz(h,1,256);


The figure below shows the magnitudes of the two frequency responses (IIR and FIR): • "a (linear-phase) FIR filter design routine" Can you elaborate on that please? Do you know some websites/papers describing precisely a method to do that? – filaton Aug 31 '17 at 9:04
• McClellan, J., T. W. Parks, and L. Rabiner. "A computer program for designing optimum FIR linear phase digital filters." IEEE Transactions on Audio and Electroacoustics 21.6 (1973): 506-526. – user28715 Aug 31 '17 at 21:37
• Rabiner, L. "Linear program design of finite impulse response (FIR) digital filters." IEEE Transactions on Audio and Electroacoustics 20.4 (1972): 280-288. – user28715 Aug 31 '17 at 21:37
• Wu, Shao-Po, Stephen Boyd, and Lieven Vandenberghe. "FIR filter design via semidefinite programming and spectral factorization." Decision and Control, 1996., Proceedings of the 35th IEEE Conference on. Vol. 1. IEEE, 1996. – user28715 Aug 31 '17 at 21:38
• Ababneh, Jehad I., and Mohammad H. Bataineh. "Linear phase FIR filter design using particle swarm optimization and genetic algorithms." Digital Signal Processing 18.4 (2008): 657-668. – user28715 Aug 31 '17 at 21:39

Simple solution:

1. Sample the impulse response of the IIR with sufficient length, 8192 or so should be plenty in this case
2. FFT
3. Set phase to zero
4. Inverse FFT
5. Time shift and truncate to desired accuracy/filter-length

EDIT: here is the code how to do it

%% get a filter target
%sos = audioEQ(6,5000,sqrt(.5),'para')
fs = 44100;
% paramtric: 6 dB, 5 kHz, Q = 1
a = [1.000000000000000  -1.216444449798070   0.533294672146362];
b = [1.232247112503961  -1.216444449798070   0.301047559642400];

%% Go through it step by step
% 1. Sample the impulse response
nx = 8192;
delta = zeros(nx,1); delta(1) = 1;
hiir = filter(b,a,delta);
% 2. FFT
fh = fft(hiir);
% 3. Set phase to zero
fhZeroPhase = abs(fh);
% 4. inverse fft
hfir = ifft(fhZeroPhase);
% 5. cut and shift to desired size. Let's go with 63 tabs
nFinal = 63;
hFinal = circshift(hfir,(nFinal-1)/2);
hFinal = hFinal(1:nFinal,:);

%% Plot the difference between the two spectra
freqAxis = (0:nx/2)'/nx*fs;
fDiff = fft(hFinal,nx)./fh;
semilogx(freqAxis,20*log10(abs(fDiff(1:nx/2+1))));
ylabel('Error in dB');
xlabel('Frequency in Hz');
set(gca,'xlim',[20 20000]);
% Note the scale: they magtnitude of the filters matches better than 1.2e-4
% dB. Depending on how good your match needs to be, you can probably get
% away with a much shorter filter

• A POCS on that approach! :D (POCS = Projections Onto Convex Sets). Does this work? You might need to repeat 2->3->4 several times before it converges. – Peter K. Aug 31 '17 at 13:59
• @PeterK.: This works just fine on the first attempt. No convergence involved. The only source of error is the truncation or wrap-around of the impulse response. – Hilmar Aug 31 '17 at 14:07
• I thought about that method but wanted to avoid going through FFT & IFFT, but I might end up doing that anyway – filaton Aug 31 '17 at 14:09
• @Hilmar I tried that method and had some problem. I asked a new question (see there) where people say that it is not the correct method. Do you have an idea? – filaton Sep 2 '17 at 11:59
• @filaton: example code added – Hilmar Sep 2 '17 at 19:15