Custom nearly-linear-phase IIR filter computation gone wrong
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