mpiir_l2(): Results are not IIR but FIR filters - why and how to avoid this?

I'm using mpiir_l2() from user Matt L's PhD thesis to design IIR filters. I set the number of numerator and denominator coefficients both to the same value (between $$100$$ and $$300$$). The maximum pole raduis is $$0.98$$. My observation is that the result is effectively a FIR filter, with the denominator coefficients allways being $$[1, 0, 0, ...]$$. Basically, the result is equal to that of the lslevin() function for least squares FIR design from the same source.

Clearly something in my filter spec makes the FIR solution have better errors than any IIR filter. The specification is the result of numerical optimization and looks pretty random to a humans eye, so I'm not sure what properties to look for.

Is there a way to get a true IIR filter from this function? I want to compare IIR anf FIR filters for the same specification and see what's better. Right now, I can't really compare when the results are basically the same.

EDIT: Here is a sample script and corresponding desired response. It seems like I have to use a third party hosting site for this, sorry.

close all;
clear all;

% loads the specification, weights and frequencies as three vectors
% the sepcification has a relevant band that is padded with zeros to both
% sides. There is a "don't care" gap between the relevant band and the
% zeros.

% search area
allowedOrder = 1:1:3;       % numerator and denominator order
allowedShifts = 0:0.125:2;      % linear phase term

% other variables
Fs = 48000;
normalFreqs = freqs / Fs * 2 * pi; % normalized frequencies

% preallcoate for the search
lowestError = inf;
bestNum = [];
bestDenom = [];
errors = zeros(length(allowedOrder), length(allowedShifts));
orders = zeros(length(allowedOrder), length(allowedShifts));
shifts = zeros(length(allowedOrder), length(allowedShifts));

% search for best aproximation
for nl = 1:length(allowedOrder)
ord = allowedOrder(nl);
for ns = 1:length(allowedShifts)
orders(nl, ns) = ord;
shifts(nl, ns) = allowedShifts(ns);
% add the linear phase shift
phaseShiftTerm = exp(-1i * shifts(nl, ns) * (ord-1)/2/Fs * 2 * pi * freqs);
D = desiredResponse .* phaseShiftTerm;

% design the filter (change denominator order here, if you like)
[num,denom] = mpiir_l2(ord, ord, normalFreqs, D, Ws, poleRadius);

% calculate filter response
H = freqz(num, denom, freqs, Fs);
% calculate the weighted error
errors(nl, ns) = sum(Ws .* abs(H - D).^2);

if errors(nl, ns) < lowestError
lowestError = errors(nl, ns);
bestNum = num;
bestDenom = denom;
end
end
end

• I can't tell you exactly what is happening for your specs, but I'm actually positively surprised that the result is at least not completely useless. An order of several 100 is totally unrealistic for an IIR filter, and even if it could be designed (which I don't think is the case with the current algorithms), it is very likely that it can't be implemented in a numerically stable way. – Matt L. Nov 19 '18 at 11:42
• In my experience, the best results are achieved with just a few non-zero denominator coefficients, say 4 to 8, and as many numerator coefficients as necessary to achieve an acceptable approximation error. What is crucial is the determination of the best average delay of the desired phase response. It must match the chosen filter order, and it can only be determined experimentally. – Matt L. Nov 19 '18 at 12:00
• Interessting. I had planned to split the filter into biquads with tf2sos (I'm reading just now that this is unstable for high orders as well). With double precision floating point math I wouldn't have expected the design algorithm to be numerically problematic. But I don't know much about the internal workings of it, so it's interesting to hear I'm doing something unexpected with this algorithm. Re delay: I added an aritifcial delay to my desired response, and chose the best delay out of several options in the vincinity of $(N-1)/2$ – TheSlowGrowth Nov 19 '18 at 12:11
• A delay of $(N-1)/2$ samples works well for FIR filters, but it might not be optimal for IIR filters. – Matt L. Nov 19 '18 at 12:14
• I edited your question by adding your code. I hope thats okay for you. – Irreducible Nov 20 '18 at 6:54

I had a look at your specs, and I designed a few filters to see what is going on. First of all, we shouldn't expect that an IIR filter should perform much better than an FIR filter for the given specifications, because poles are mainly useful if sharp transitions from pass bands to stop bands must be realized, as is the case for frequency selective filters. In your case, we don't have sharp transitions, and adding poles is unlikely to improve the design substantially. Of course, theoretically, the approximation error should be decreased by adding poles (if the numerator order is not changed compared to the FIR solution), but the improvement over the FIR filter may be small. In practice, however, the resulting IIR filter could even be worse than the FIR filter with the same (numerator) order. There are two reasons for that: first, the design problem is much more difficult for the IIR filter and we might only find a local optimum, and, secondly, we need to restrict the pole radius to some value $$r<1$$ in order to get a stable filter that can be implemented. That restriction may increase the possible approximation error.
It indeed turns out that using an IIR filter does not improve things very much in your case. I used an FIR filter of length $$N=100$$ and an IIR filter with the same number of numerator coefficients, and $$6$$ additional poles away from the origin of the $$z$$-plane. The magnitudes of the frequency responses are shown in the top figure below (blue: desired, green: FIR, red: IIR). The IIR approximation is slightly better, but the difference might not be worth it. The bottom figure zooms in the lower frequencies, where the IIR filter approximates the desired response a bit better. The phase responses (not shown) are comparable for both filters. Note that for the FIR and the IIR designs I individually optimized the necessary linear phase term that was added to the given phase response.
The pole radii of the IIR filter are between $$0.89$$ and $$0.98$$ ($$r=0.98$$ was chosen as the maximum pole radius in the design process). So the designed IIR filter is really IIR, not just FIR with all its poles (approximately) at the origin. I do not know why you only obtain IIR designs that are actually FIR.