Custom nearly-linear-phase IIR filter computation gone wrong

Question

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?

http: vadkudr.org/Algorithms/LangIIR/LangIIR_examples.html

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(-1iom(1:20)tau),zeros(1,77)]; W=[ones(1,20),1000ones(1,77)]; [b,a]=mpiir_l2(M,N,om,D,W,r); — ionone, Mar 23 '15 at 11:22
I think you got it wrong: that example has $M=40$ ZEROS and $N=6$ poles. So I wonder how you get 46 a coeffs and 11 b coeffs. That's the problem. — Matt L., Mar 23 '15 at 11:24
with poles = 6, M = 40 and group delay (tau) = 40 I still get huge coeffs, but they work in fvtool()! that's weird, but I suppose as the coeffs go higher, the errors add up far faster than with lower coeffs...thats should be the problem, because IIR uses the output in the computation of the next output, so each little error adds up and at the end it is huge! — ionone, Mar 23 '15 at 11:29
You can't just change the numbers of poles and zeros without changing the desired frequency response. That will result in a bad design. However, the result will always be stable, unless the routine fails to converge, which may happen for very unrealistic specifications. There is probably another problem in your implementation. — Matt L., Mar 23 '15 at 11:30
I added a design example (#7 from Vadim's site with tau=40) to my answer. It works fine, filter coefficients can be found below. — Matt L., Mar 23 '15 at 14:06
And about the fvtool: if you just look at the magnitude of the frequency response you don't see what's going on. You could e.g. look at the pole-zero plot, then you'd see that there are poles outside the unit circle (for the wrong coefficient set). You can also check the impulse response, and you'll see that it's not decreasing but increasing. So it all depends on which things you actually check in fvtool. — Matt L., Mar 23 '15 at 18:05

Matt L. · Accepted Answer · 2015-03-23T14:05:19.650

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.

enter image description here

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

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

Custom nearly-linear-phase IIR filter computation gone wrong

1 Answers1

Linked