# Find filter coefficients to model a device using its measurement

I am trying to realize a digital filter that has the same freq. response of an existing speaker. I have fed an audio sine sweep to the speaker and measured the speaker output, both at 48kHz. Then I perform FFTs of the input(X) and output(Y), divide the absolute values point-wise (absY./absX) to get the transfer function/freq. response (H).

Now I would like to determine the filter coefficients B and A for an FIR/IIR filter, so that I can model the speaker's response digitally. I understand my transfer function is in the domain of w, and filter coefficients are the coefficients of the difference equation of the filter.

This is where all my text-book theory seems to fall apart, and I was hoping someone could clarify:

1. If I understand correctly, since my data is discrete (sampled at 48k), my data is in terms of 'n'. If I performed an ifft of H (which is complex), the resulting vector is basically my filter coefficients?
2. Matlab suggests using invfreqz to find filter coefficients, given a complex frequency response. Does this function map from w to z domain? If so, am I right in understanding that the FFT needs to be converted to Z domain to derive its filter coefficients?
3. What parameters does a auto-regressive filter fit to? for eg: Yule_Walker, or Levinson, or if I want to run a gradient descent to fit a filter? What is the error calculated between?
• The simplest approach would be to truncate the impulse response $h[n]$ to some length with a windowing function and to use it directly as coefficients of an FIR filter. – hulappa Oct 8 '18 at 7:46
• Hi! This is a very exciting topic. Have you searched the existing literature on system modeling (identification) ? Especially those of adaptive systems that rely on LMS, RLS or Kalman filtering? Also you should in the very first place note that a loudspeaker is not strictly an LTI system; i.e., cone displacement is only approximately linear function of excitation amplitude... – Fat32 Oct 8 '18 at 9:13
• @hulappa My concern is along the lines of whether the result of $h[n]$ is the same, whether I use ifft or invfreqz (which I've tested in MATLAB, and doesn't appear to be so), and not to mention, invfreqz fits a desired number of poles and zeros to the given response. So I was wondering how that happens exactly. – Aditya TB Oct 9 '18 at 15:21
• @Fat32 Thanks for your reply. Yes I have looked at System ID, and I understand that the speakers aren't really LTI, and in fact, my speaker is super non-linear. I am making a concession by assuming the speaker is LTI, so as to only affect the magnitude response. There seem to be a lot of NonLinear methods too (NARX, or ARMA type filters), which I am trying to build my knowledge towards. Though my methods maybe primitive, they help me learn the basic steps before getting to the NonLinear ID stage. – Aditya TB Oct 9 '18 at 15:21
• The results of the IFFT-based approach and invfreqz are not the same. invfreqz finds coefficients of an IIR filter using some sort of optimisation procedure which I don't know the details of, whereas the IFFT-based method I described simply uses the impulse response as coefficients of an FIR filter. Filtering with this FIR filter is the same as convolving your input signal with the (truncated) impulse response you computed with the IFFT-based method. – hulappa Oct 9 '18 at 16:23

We can represent a digital filter as:

• A difference equation, relating input and output at every time step $$n$$: $$y[n] = \sum_{k=0}^{K-1} b[k] x[n-k] - \sum_{m=1}^{M} a[m] y[n-m]$$. The filter coefficients are $$b_k$$ and $$a_m$$. In case of an FIR, $$a[m] = 0\,\forall\, m$$. This form is basically how the filter can be implemented (in the discrete time domain).
• By z-Transform of the difference equation, we can also represent the above equation as: $$Y(z) = \frac{\sum_{k=0}^{K-1} b[k] z^{-k}}{1 + \sum_{m=1}^{M} a[m] z^{-m}} X(z) = H(z)X(z)$$.
• The frequency response is simply this z-Transform evaluated at the unit circle $$H(z = \mathrm{e}^{j\omega})$$ (this is often simply written as $$H(\omega)$$).
• The impulse response is the inverse transform of $$H(\mathrm{e}^{j\omega})$$. By convolution with the impulse response, we also can get the output of the filter, as in $$y[n] = h[n] * x[n]$$. However, $$h[n]$$ is only finite for an FIR (hence the name, obviously), where we can easily check that in fact it is simply the coefficients $$b[k]$$.

I hope this helps you with understanding the big picture. So to answer your questions:

1. The IFFT of the frequency response yields the impulse response. This corresponds directly to the filter coefficients only in the case of an FIR filter. In case you want to model your system as an IIR, you can use an algorithm such as the kind implemented in MATLAB's invfreqz, which is essentially an optimization to fit your frequency response to the complex rational function $$H(z)$$.

2. Yes, since $$H(z)$$ is a representation of the filter's difference equation, the coefficients of which are what we call filter coefficients.

In any case, these are important basics you can learn from e.g. the classic book on DSP. Even the MATLAB documentation (which im sure you have already been checking) is a decent resource.

• Thanks @mateC. This helps a lot. If you don't mind, I would like just a little more clarification. – Aditya TB Oct 9 '18 at 15:23
• Pushed enter accidentally. The problem I have is translating the textbook equations, into code. 1. I can see that invfreqz uses a form of curve-fitting approach. Am I right in understanding that minimizing error is the only real way of producing b and a coefficients? Or is there some other approach like Finite Diff. Approximation or Bilinear transform to obtain b and a, while mapping $\omega$ to $z$. 2. Am I right in understanding that the z-domain only offers the power to analyse stability, but for filter implementation, one needs the time-domain representation? – Aditya TB Oct 9 '18 at 16:09
• 1. Since you already have a measurement of the frequency response of your loudspeaker, you can get the coefficients of an approximately equivalent digital filter by estimating from it. This estimation is what the curve-fitting does. Bilinear transform and finite differences are different things, don't know what you mean. 2. $H(z)$ offers a representation of the frequency response, as well as the time-domain difference equation, as mentioned before. Oh btw, you could also implement the filter directly in frequency domain, avoiding the coefficient estimation. Depends on what you wanna do... – mateC Oct 10 '18 at 9:07