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?
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 &gt; 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.