# How to estimate a transfer function from a magnitude-only frequency response?

Given an arbitrary frequency response, what signal processing methods might exist that could guess, estimate or determine a transfer function (pole and zero constellation) which gives a "reasonably good" approximation (for some given estimation quality criteria) to that given frequency response? What means exist to estimate the number of poles and zeros required for a given transfer function plus a given approximation error allowance? Or how could one determine these constraints can't be met, if possible?

If the given frequency response was actually produced by a known transfer function, will any of these methods converge on that original transfer function? How about if the given frequency response were subject to (assumed Gaussian) measurement errors?

Assume working in the Z-plane with sampled spectrum, although continuous domain answers might also be interesting.

Added: Are the solution methods any different if only the magnitude of the frequency response is given (e.g. a solution with any phase response is allowed)?

Added: The latter problem is what I'm most interested in, given a known magnitude response around the unit circle, but unknown/unmeasured phase response, can the measured system be estimated, and if so under what conditions?

• Are you attempting to approximate an arbitrary frequency response as a rational spectrum? That is (b[0]+b[1]z^-1...)/(1+a[1]z^-1...)? If so, this is typically referred to as ARMA modeling. It's more difficult than AR modeling because the autocorrelation of a signal tends to be nonlinearly related to the moving average coefficients (the b[]'s, or zeros). If my assumption is correct I can write a more formal response. – Bryan May 23 '12 at 21:22
• @Bryan : Yes. I tried to imply that by stating a "pole and zero" solution (a rational transfer function) was suitable (preferably only if better than an all pole or all zero solution/estimate of the same degree). – hotpaw2 May 23 '12 at 21:34
• What meaning is attached to frequency response? Some people distinguish between the frequency response function $H(\omega)$ or $H(f)$ and the transfer function $H(s)$ and some people don't. See for example, the discussion following this answer to an earlier question. – Dilip Sarwate May 23 '12 at 21:34
• @Dilip Sarwate : Given H(w) only for the unit circle (is that redundant?), solve/estimate a complete z plane representation. Hopefully, that's consistant with my original statement of the question. – hotpaw2 May 23 '12 at 21:49
• You addition changes things. Poles and zeros can change with the magnitude response remaining the same. The most common example of this is when one is designing a minimum phase filter. This typically involves taking an existing system and reflecting the poles and zeros inside the unit circle. This only changes the phase response, not the magnitude response. – Bryan May 23 '12 at 22:30

One approach would be to use the frequency-domain least-squares (FDLS) method. Given a set of (complex) samples of a discrete-time system's frequency response, and a filter order chosen by the designer, the FDLS method uses linear least-squares optimization to solve for the set of coefficients (which map directly to sets of poles and zeros) for the system whose frequency response matches the desired response with minimum total squared error.

The frequency response of an $N$-th order linear discrete-time system can be written as:

$$H(\omega) = H(z)|_{z=e^{j\omega}}$$

where $H(z)$ is the system's transfer function in the $z$ domain. This is typically written in the rational format that follows directly from the system's difference equation:

$$H(z) = \frac{\sum_{k=0}^{N}b_kz^{-k}}{1 + \sum_{k=1}^{N}a_kz^{-k}}$$

The frequency response is therefore:

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

Rearrange the above to get:

$$\sum_{k=0}^{N}b_ke^{-jk\omega} - H(\omega) \left(1 + \sum_{k=1}^{N}a_ke^{-jk\omega}\right) = 0$$

This equation is linear in the $2N+1$ unknown system difference equation 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.

In order to solve for an appropriate set of coefficients using the linear least squares method, we generate an overdetermined system of equations in those unknowns. To generate those equations, choose a collection of frequencies $\omega_m \in [0, 2\pi), m = 0, 1, \ldots , M-1$ (where $M > 2N+1$, and often $M \gg 2N+1)$. For each frequency, substitute the corresponding value of $\omega_k$ into the above equation to yield:

$$\sum_{k=0}^{N}b_ke^{-jk\omega_k} - H(\omega_k) \left(1 + \sum_{k=1}^{N}a_ke^{-jk\omega_k}\right) = 0$$

The values $H(\omega_k)$ are obtained by sampling the desired frequency response at the chosen frequencies $\omega_k$. After generating the system of linear equations, the least-squares solution for the system's coefficients $b_k$ and $a_k$ (and therefore its poles and zeros) is easily obtained. If you substitute those coefficients back into the equation for $H(\omega)$ shown above, it should (hopefully) yield a function that is close to the template frequency response that you started with.

This technique has a few advantages:

• Any arbitrary complex (magnitude and phase) frequency response can be used as the template. If you only have a magnitude constraint, you could just pick a phase response, such as linear phase.

• It can be used to design both FIR and IIR filters; for an FIR realization, just remove the $a_k$ coefficients from the above.

• The technique is very simple to implement and is easily parameterizable based upon the desired system order.

• While there may not be a good way to estimate a priori what the required system order is to meet your design constraints, it is simple to iteratively increase the order $N$ until some selected error metric is met (such as peak error, total squared error, or deviation within a specific band).

You could extend this method a bit to use weighted least-squares optimization if needed; this would allow you to specify regions of the frequency response whose approximation error is weighted more than others. This allows you to more tightly control passband/stopband areas while allowing more slop in "don't-care" areas.

• Excellent answer!! The "art" in doing filter designs with least square error is to properly defined what exactly the "error" is. This is controlled by choosing the right frequency grid, weighting factors at specific frequencies and adding more constraints for out of band behavior and also for keeping your poles inside the unit circle. – Hilmar May 24 '12 at 11:01
• The problem with this potential solution is, if the phase is unknown about an existing transfer function, FDLS may converge on the wrong solution if the wrong phase is assumed, no matter how accurately the order is correctly guessed or the magnitude response is measured. – hotpaw2 Mar 12 '13 at 16:01
• @hotpaw2: That's to be expected. If you don't know anything about the phase response, then there are an infinite number of solutions that are equally valid (i.e. they would have the correct magnitude response). You're going to need some information to steer you toward what you deem to be the most appropriate solution. – Jason R Mar 12 '13 at 16:08
• @JasonR : The only correct solutions should be the permutations of flipping poles/zeros inside/outside, which is a finite number for any (existing) finite order system. – hotpaw2 Mar 12 '13 at 16:50

My colleagues have had great results with vector fitting:

Vector Fitting is a robust numerical method for rational approximation in the frequency domain. It permits to identify state space models directly from measured or computed frequency responses, both for single or multiple input/output systems. The resulting approximation has guaranteed stable poles that are real or come in complex conjugate pairs.

We use it for FIR to IIR conversion.

For less demanding applications, you can just use nonlinear least squares fitting for a fixed number of poles and zeros. This is implemented in Matlab as invfreqs and invfreqz.

Another approach: plot the frequency response and fit a Bode plot to it as best as possible. This could be done very quickly for an approximate solution, or in some elaborate least squares sense for a better fit. GTH