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

Question

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?

Answer 1

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):

Answer 2

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