5
$\begingroup$

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:

enter image description here

What have I been doing wrong?

$\endgroup$
2
  • 1
    $\begingroup$ 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. $\endgroup$
    – fibonatic
    Sep 2, 2017 at 14:24
  • $\begingroup$ Is having a mixed bag allowed? If so, use your IIR with an allpass FIR, for equalization. $\endgroup$ Sep 3, 2017 at 7:12

2 Answers 2

8
$\begingroup$

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

enter image description here

$\endgroup$
12
  • $\begingroup$ "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? $\endgroup$
    – filaton
    Aug 31, 2017 at 9:04
  • $\begingroup$ 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. $\endgroup$
    – user28715
    Aug 31, 2017 at 21:37
  • $\begingroup$ Rabiner, L. "Linear program design of finite impulse response (FIR) digital filters." IEEE Transactions on Audio and Electroacoustics 20.4 (1972): 280-288. $\endgroup$
    – user28715
    Aug 31, 2017 at 21:37
  • $\begingroup$ 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. $\endgroup$
    – user28715
    Aug 31, 2017 at 21:38
  • $\begingroup$ 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. $\endgroup$
    – user28715
    Aug 31, 2017 at 21:39
4
$\begingroup$

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
$\endgroup$
10
  • $\begingroup$ 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. $\endgroup$
    – Peter K.
    Aug 31, 2017 at 13:59
  • 2
    $\begingroup$ @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. $\endgroup$
    – Hilmar
    Aug 31, 2017 at 14:07
  • $\begingroup$ I thought about that method but wanted to avoid going through FFT & IFFT, but I might end up doing that anyway $\endgroup$
    – filaton
    Aug 31, 2017 at 14:09
  • $\begingroup$ @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? $\endgroup$
    – filaton
    Sep 2, 2017 at 11:59
  • 1
    $\begingroup$ @robertbristow-johnson: not in this case. The tails of the impulse response are decaying exponentially already and the vast amount of the energy is in the center samples. In this example applying a hanning window increases the error 10e-4 dB to 10e-1 dB. $\endgroup$
    – Hilmar
    Sep 3, 2017 at 11:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.