Using MATLAB's invfreqz to fitting a audio frequency response, but the result is bad at low frequencies

Question

I am a DSP beginner. I have some audio frequency response (actually acoustic, especially speaker frequency response by the way) curves and I'm trying to fitting them by IIR filters using MATLAB's invfreqz.

The curves I have are just magnitude data and I use JOS's method to compute the minimum phase in addition to magnitude. Then I use invfreqz to get some the IIR numerator and denominator coefficients. The desired IIR order is 8~12.

I get correct results, but I find the fitting result is good in high frequency and bad in low frequency. Because you know in audio we usually plot curves in a logarithmic x axis, and the curves have more details in low frequency. How can I improve low frequency results?

These two figures are my fitting results in logarithmic x axis and linear x. The IIR orther is 12 both in numerator and denominator. You can see the details in low frequency are a lot, but the IIR does not match well.

I tried increase order to hundreds the filter matches very very good, but that is unapplicable in a real life dsp core.

My target data:
f/Hz 20 21.20000076 22.4 23.6 25 26.5 28 30 31.5 33.5 35.5 37.5 40 42.5 45 47.5 50 53 56 60 63 67 71 75 80 85 90 95 100 106 112 118 125 132 140 150 160 170 180 190 200 212 224 236 250 265 280 300 315 335 355 375 400 425 450 475 500 530 560 600 630 670 710 750 800 850 900 950 1000 1060 1120 1180 1250 1320 1400 1500 1600 1700 1800 1900 2000 2120 2240 2360 2500 2650 2800 3000 3150 3350 3550 3750 4000 4250 4500 4750 5000 5300 5600 6000 6300 6700 7100 7500 8000 8500 9000 9500 10000 10600 11200 11800 12500 13200 14000 15000 16000 17000 18000 19000 20000
data 0.000410959 -0.246575342 -0.01973 -0.05 0.01589 -0.02041 -0.25849 0.80863 1.260959 0.775342 0.458356 -0.03849 -0.45959 -0.54151 -0.54836 -0.53767 -0.49164 -0.35877 -0.18219 -0.13493 -0.18945 -0.35041 -0.6111 -0.7774 -0.86671 -0.82945 -0.71342 -0.64767 -0.5426 -0.38425 -0.17877 0.072329 0.262877 0.30411 0.244247 0.096575 -0.0989 -0.2763 -0.39178 -0.45082 -0.47055 -0.48192 -0.49836 -0.51123 -0.52315 -0.54014 -0.58849 -0.64671 -0.68616 -0.68178 -0.66274 -0.68315 -0.68164 -0.60274 -0.43 -0.21589 -0.14849 -0.22699 -0.32945 -0.41027 -0.48822 -0.58096 -0.66726 -0.73274 -0.76014 -0.73521 -0.69425 -0.63822 -0.59712 -0.57425 -0.57356 -0.58603 -0.58658 -0.58507 -0.54658 -0.42507 -0.32603 -0.23685 -0.23534 -0.27932 -0.29822 -0.30699 -0.33151 -0.39151 -0.41671 -0.36096 -0.25644 -0.16356 -0.13178 -0.18274 -0.28616 -0.43055 -0.58753 -0.74753 -0.8437 -0.80877 -0.7574 -0.75521 -0.81 -0.90603 -0.97795 -0.93178 -0.82123 -0.79397 -0.95356 -1.33466 -1.86589 -1.99534 -1.25178 -0.37795 0.255342 0.595068 0.634521 0.568219 0.452466 0.226438 0.093288 0.194384 -0.13753 -3.46616 -3.14658

Answer 1

Your order is too low to obtain a fine result within whole frequency range. However there may be some ways to go.

---

First one is FDLS, which stands for frequency-domain least-squares.

The transfer function of an IIR filter is given by $$ H(z) = \frac{\sum_{k=0}^{P}b_kz^{-k}}{1 + \sum_{k=1}^{Q}a_kz^{-k}}\tag{1} $$

The frequency response is therefore:

$$ H(\omega) = \frac{\sum_{k=0}^{P}b_ke^{-jk\omega}}{1 + \sum_{k=1}^{Q}a_ke^{-jk\omega}}\tag{2} $$

Rearrange the above equation to get:

$$ \sum_{k=0}^{P}b_ke^{-jk\omega} - H(\omega) \left(\sum_{k=1}^{Q}a_ke^{-jk\omega}\right) =H(\omega)\tag{3} $$

This equation is linear with $P+Q+1$ unknown coefficients $b_k$ and $a_k$. Given a desired frequency response $H(\omega)$, we would like to find coefficients that meet the above equation exactly for all values of $\omega$. For the general case, that's hard. So instead, we will search for a set of coefficients for a system whose frequency response approximates the desired response at a discrete set of frequencies.

Then we choose a collection of frequencies $\omega_m \in [0, \pi], m = 1, 2, \ldots , M$ (where $M > P+Q+1$, and often $M \gg P+Q+1)$. For each frequency, substitute the corresponding value of $\omega_m$ into the above equation to yield:

$$ \sum_{k=0}^{P}b_ke^{-jk\omega_m} - H(\omega_m) \left(\sum_{k=1}^{Q}a_ke^{-jk\omega_m}\right) = H(\omega_m)\tag{4} $$

which yields a set of $M$ linear equations and can be solved by matrix inversion.

In order to have a better approximation at lower frequencies, you may sample more frequency points at low frequency, such as sampling in log space. Also, you may set different weight for each linear equation, bigger weight for lower frequency.

To make sure the filter coefficients are real, you can use the real part of the linear equations.

---

Another way is fixed-pole parallel filter design.

In FDLS we have to solve both numerator and denominator coefficients, while in the second method the poles of the filter is placed logarithmically and slightly inside the unit circle, remaining only numerators to be calculated. Thus the problem is convex in terms of MSE and can be solved by least-squares in either time domain and frequency domain.

---

Recently, I noticed a deep learning method to design IIR filter directly from the desired magnitude response. You can also have a try.

Answer 2

invfreqz() generally doesn't work for "audio type" filters. Due to the log spacing of the frequency axis, most audio filters have the majority of poles extremely close to $z=1$. That means that wiggling a pole or a coefficient by a tiny amount has a massive impact on the frequency response.

Just to see how bad it is you can look at a simple 5th order Butterworth high pass at 40Hz at a sampling rate of 48kHz, i.e. [z,p,k] = butter(5,40*2/48000,'high');

Wiggling a single numerator coefficient by as little as $10^{-14}$ limits your stop band attenuation to $50 dB$ . Even worse if you wiggle a denominator coefficient more than about $5\cdot 10^{-11}$ the thing becomes unstable.

Most fitting algorithms can't deal with this type precision required, so any algorithm that's operating in numerator/denominator representation is bound to give very poor results.

Audio filter IIR fitting can be done, but it typically requires

1. Representation as poles and zeros, not numerator and denominator
2. Some time of frequency warping to get the poles further away from $z=1$
3. An iterative search algorithm (steepest decent, conjugate gradient, etc.)
4. Intelligent seeding of poles and zeros
5. Intelligent management of stability and "maximum Q" constraints

Update

Here is a 12 pole-pair IIR fit based on the OPs data