# Custom nearly-linear-phase IIR filter computation gone wrong

I'm using the code in this book: Algorithms for the constrained design of digital filters with arbitrary phase and magnitude responses to compute nearly-linear-phase IIR EQs.

fvtool() shows the correct Eq mag and phase, but when I try to apply the coefficients to some real sound, the output is very very high (silence due to overlimit), it seems very unstable.

Anyone could tell me why fvtool() is ok but filter() is wrong, maybe try the software?

Jeff

EDIT:

Try these coefficients:

a =
1.00
-7057018074799500.00
7419034985051048.00
-2710509977123811.00
3321548070644472.50
19380841977826980.00
-10366336133286466.00
-17395746351478384.00
-241875435346629.44
1537628108108181.70
-7101604243654517.00
7251792090687912.00
23447717670347684.00
-11610491157524712.00
-10852963123569024.00
5000000000000001.00

b =
-65252253992.00
29593281937568.00
31626970863960.00
-7964372882016.00
-32245859702352.00
-25890947971080.00
-10668161380760.00
-4361278103120.00
11507726488128.00
24495364570416.00
6689249360480.00


and the code to write the ouput file:

input= wavread('sine+noiz.wav');
output2=filter(b,a,input);
output2=output2/norm(output2);
wavwrite(output2, 44100,  'e:\sine+noiz2_out.wav');


EDIT 2:

Here is full working code and audio files

EDIT 3:

Okay I figured out where the problem was coming from: to make a lowpass at a very low frequency, you have to add more Zeros : the lower the cutoff frequency, the higher the number of zeros necessary...so then it works But I still don't know why fvtool() shows an ok filter but then the actual computation of an audio file fails...

• Could you post the filter coefficients? How did you implement the filter, with Matlab's filter() function? What's the filter order? – Matt L. Mar 22 '15 at 15:37
• It would help if you added your Matlab code for performing the filtering. – Matt L. Mar 22 '15 at 20:41
• OK, these coefficients are nonsense, so something went wrong in the design routine. Are you sure that a are the denominator coefficients, i.e. you want 45 poles??? Also, these coefficients do not correspond to any of the examples on Vadim's site. If you give me the design specs, I can come up with a reasonable design. The specification of the desired phase requires some experience. – Matt L. Mar 23 '15 at 11:19
• yes I want 45 poles, try fvtool() it works. To obtain linear phase you have to higher the poles and zeros that's the only way you can have linear phase and high order – ionone Mar 23 '15 at 11:21
• this example has a lot of poles/zeros also and it is in the examples page: tau=37;M=40;N=6;r=0.98; om=pi*[linspace(0,0.2,20),linspace(0.23,1,77)]; D=[exp(-1i*om(1:20)*tau),zeros(1,77)]; W=[ones(1,20),1000*ones(1,77)]; [b,a]=mpiir_l2(M,N,om,D,W,r); – ionone Mar 23 '15 at 11:22

I believe that there's something going wrong with the implementation of the filter. The design algorithm lets you prescribe the maximum pole radius, so if you prescribe e.g. $r_{max}=0.98$ the filter is guaranteed to be stable with some reasonable stability margin. So using floating point arithmetic it seems to me that nothing can go wrong. If you use Matlab's filter() function, be careful to call it filter(b,a,x) with b being the numerator coefficients. I've seen people use this function as filter(a,b,x) (which seems more logical to some folks). The latter function call will give you trouble, because filter() will interpret the numerator coefficients b as denominator coefficients, so any zero outside the unit circle will become an unstable pole. And due to the phase approximation of the design algorithm there will almost always be some zeros outside the unit circle.

EDIT:

This is example #7 from this site, with the desired group delay value set to 40 samples (instead of 37). This is what I understood from your comments you're trying to do. This is the code:

tau=40;M=40;N=6;r=0.98;
om=pi*[linspace(0,0.2,20),linspace(0.23,1,77)];
D=[exp(-1i*om(1:20)*tau),zeros(1,77)];
W=[ones(1,20),1000*ones(1,77)];
[b,a]=mpiir_l2(M,N,om,D,W,r);

[H,w] = freqz(b,a,1024);
gd = grpdelay(b',a',1024);
subplot(2,1,1), plot(w/2/pi,20*log10(abs(H)));
grid, title('|H(f)|')
subplot(2,1,2), plot(w/2/pi,gd);
grid, title('passband group delay'), axis([0,.1,35,45])


I ran the design program in Octave. I believe that you have a problem with running the software in the current Matlab version. Note that the software was written about 17 years ago, and I'm not sure if and how you (or someone else) adapted it to the current version. You can find the filter coefficients below and I'm convinced that a filter implemented with these coefficients will work fine. These are the coefficients:

b =

-5.3754e-04
1.4765e-03
-1.7103e-03
1.1605e-03
-3.7906e-04
6.8635e-05
-2.9013e-05
-2.0759e-05
-4.4175e-05
-1.0035e-05
-8.8312e-06
3.8369e-05
3.4096e-05
5.6969e-05
1.6945e-05
5.8498e-06
-5.1954e-05
-5.4722e-05
-7.7327e-05
-2.8043e-05
-1.5210e-07
7.7270e-05
9.5478e-05
1.2200e-04
5.9885e-05
1.6209e-06
-1.1856e-04
-1.8065e-04
-2.3618e-04
-1.7031e-04
-5.9434e-05
1.7262e-04
4.0741e-04
6.8756e-04
8.7908e-04
9.9248e-04
1.3572e-03
-1.4580e-04
3.0311e-03
-1.8688e-03
1.9026e-03

a =

1.00000
-4.65499
9.58904
-11.07801
7.53566
-2.85503
0.47033

• sorry I removed the "completed" thread because some examples work and some don't (in the webpage I provided). So it doesn't come from my code. It seems that the computations are going wrong as the order goes up (small order filters are ok, but higher outputs silence). And the higher a and b, the more likely it is to go awry – ionone Mar 23 '15 at 10:28
• @ionone: Can you post the filter coefficients of one example where things go wrong? I still believe that it is not the actual computation of the filter, because as long as the filter is stable with a reasonable stability margin, everything should be fine in floating point. – Matt L. Mar 23 '15 at 10:31
• I edited the main post with filter coeffs – ionone Mar 23 '15 at 11:15
• I also normalize the output file before writing it to disk – ionone Mar 23 '15 at 11:16