matlab - 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?

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...

Answer

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