In the 1980s I wrote the matlab function invfreqz that implements a fast (FFT-based) equation-error method for recursive filter design given samples of the desired frequency response. This method designs causal filters, fitting both the phase and magnitude of the desired (complex) frequency response. As a result, one must be careful about the specified phase response, especially when it doesn't matter! This article gives some pointers on this topic.

The case I see most often is where there are no phase data at all. For example, a filter response may be specified simply as a desired (real) gain at various frequencies, i.e., only the amplitude response is specified. The easiest choice of phase response in this case is zero, but that would correspond to an impulse response that is symmetric about time zero, and causal filters cannot produce any response before time zero. In this case, phase-sensitive filter-design methods such as invfreqz will generally give their best results for minimum phase in place of zero phase. In other words, we need to convert our desired amplitude response to the corresponding minimum-phase frequency response (whose magnitude equals the original desired amplitude response).

A commonly used method for computing the minimum-phase frequency response from its magnitude is the cepstral method. The steps are as follows:

Interpolate the amplitude response samples from 0 to half the sampling rate, if necessary, and resample to a uniform "FFT frequency axis", if necessary. Denote the real, sampled amplitude response by S(k).

Perform an inverse FFT of log(S) to obtain the real cepstrum of s, denoted by c(n).

Fold the noncausal portion of c(n) onto its causal portion.

Perform a forward FFT, followed by exponentiation to obtain the minimum phase frequency response Sm(k), where now sm(n) is causal, and |Sm(k)|=S(k).

These steps are illustrated by the matlab example below (tested in Octave, version 3.2.4, on Mac OS X). Copy-paste it and give it a try!

Comments:

Hello JOS,
Thanks for posting this very interesting code. I'm going to experiment with it.
It looks like the line:
Ns = length(Gdbfk); if Ns!=Nfft/2+1, error("confusion"); end
should be:
Ns = length(Gdbfk); if Ns~=Nfft/2+1, error('confusion'); end
(I imagine those three "typos" occurred because of some sort of software text-translation problem when when your code was posted as a blog.)
Thanks again JOS,
[-Rick-]