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

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

• This response you are looking for is very irregular. IIR is for smooth curves, I would go for a FIR filter instead. But it would require a large FIR filter because you want to have a specific response down to very low frequencies, that requires very long impulse response in the time domain. Commented Feb 11, 2022 at 13:41
• The coefficients for the FIR filter is just the time-domain impulse response, that is the frequency response you have, converted to time-domain. Commented Feb 11, 2022 at 13:43
• @LeandroAlsina. FIR is the WORST choice for this. Given the low frequency details, this would require 10s of 1000s of samples. IIR is perfect for this type of thing, they are just hard to design. Commented Feb 11, 2022 at 14:06
• @ScarlettWang: can you add your target data to the post? Commented Feb 11, 2022 at 15:56
• @Hilmar: Hi Hilmar，I will add my data in question hope you can see them, a vector of frequency points and another corresponding target test data. Appreciate for your answer! Commented Feb 12, 2022 at 5:20

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.

• Sorry, that has nothing to do with the order. A 12th order IIR is plenty to fit the target response within a fraction of a dB. The problem is the log spacing, the pole locations, and the extremely steep error surface. Commented Feb 11, 2022 at 13:33
• @Hilmar The OP's target response varies significantly within the whole frequency range, so I think there must be some precision trade off between lower frequencies and higher ones with a relatively low filter order. I used FDLS and fixed-pole parallel filter for HRTF approximation, which design the filter in log space. Both methods give proper results. (I didn't measure how fine the results are but we always need a balance between accuracy and filter order) Commented Feb 11, 2022 at 14:18
• @Hilmar but yes, maybe 12th order can be sufficient but it depends on how accurate you need. Commented Feb 11, 2022 at 14:34
• HRTFs are easy since there is nothing happening below 250Hz and anything below 1kHz is very smooth and gentle. You can fit the whole thing (lows and highs) probably with 8 pole pairs or so. You just need a fitting algorithm that can handle a lot of low frequency details. Commented Feb 11, 2022 at 15:54
• @ZR Han Thank you. I know a little about the equation error method and Quasi-Newton method and DFP algorithm. I guess that is similar to FDLS? I think FDLS may be a good starting point for my problem. Commented Feb 12, 2022 at 5:32

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

• Thank you Hilmar. I accpet your explanation but I cannot totally understand for a while. Because my knowledge on DSP theory is weak and I will read books to get your point. As for the iterative algorithm, I know a little about the equation error method and Quasi-Newton method and DFP algorithm. Do you think that is the right way to solve my problem? Commented Feb 12, 2022 at 5:28
• Thanks for your update. Your fitting result is exactly what I want. Could you give me more thread about your algorithm? Commented Feb 14, 2022 at 6:06
• The outline of what it does is in my answer. Commented Feb 14, 2022 at 18:49