Derive minimum phase from magnitude

Question

With the desired magnitude of a transfer function in the frequency domain in C++ as described below what is the correct corresponding minimum phase? In general how does one derive the correct minimum phase in the frequency domain if the desired magnitude is known, specifically such that the a complex IFFT of it would return an entirely real part.

// PnCutoff.cpp : Pink noise above a cutoff.
#include <tchar.h>
#include <math.h>
#include <conio.h>

class PnAtFilter
{
    double _fcHz;
    public:
    PnAtFilter(double fcHz) { _fcHz = fcHz; }
    double Mag(double fHz)
    {
        if(fHz < _fcHz) return 1.0;
        else return  1.0 / sqrt(fHz/_fcHz);
    }
    double Phase(double fHz)
    {
        return 0.0; // what is the the correct minimum phase in r?
    }
};

int _tmain(int argc, _TCHAR* argv[])
{
    PnAtFilter* pnf = new PnAtFilter(1000);
    for(int fHz = 10; fHz <= 100000; fHz *= 10)
    {
        double mag = pnf->Mag(fHz); 
        double phase = pnf->Phase(fHz);
        _cprintf_s("%d Mag=%f, Phase=%f\r\n", fHz, mag, phase);
    }
    _getch();
    return 0;
}

/* output:
10 Mag=1.000000, Phase=0.000000
100 Mag=1.000000, Phase=0.000000
1000 Mag=1.000000, Phase=0.000000
10000 Mag=0.316228, Phase=0.000000
100000 Mag=0.100000, Phase=0.000000
*/

```

Answer 1

One method to approximate a minimum phase transfer function from a magnitude-only frequency response is to first find an suitable approximation to the transfer function in the pole-zero Z-plane domain. Then reflect all the zeros to inside the unit circle to get a minimum phase response. The first part of this is sometimes called a "system identification problem" and usually has no closed form solution. Differential evolution algorithms might be useful is finding a "close enough" pole-zero approximation (There appear to be IEEE papers on this method. See the 2005-Jan issue of IEEE Signal Processing Magazine).

Of course if you can do an exact factorization of your frequency response, then the solution is straight forward. (reflect zeros, etc.)

For "smooth magnitude responses" another approximation is to use a property of the complex cepstrum as described here by J.O.S. at CCRMA. The following script for converting an arbitrary FIR filter into a minimum phase approximation requires the frequency response not to have any zeros in order to allow approximate inversion of the real part of the cepstrum of the given magnitude response. Use of FFTs much much longer than your given magnitude response or initial FIR filter approximation may be required to get useful results.

wn = [ones(1,m); 2*ones((n+odd)/2-1,m) ; ones(1-rem(n,2),m); zeros((n+od d)/2-1,m)];
y = real(ifft(exp(fft(wn .*  real(ifft(log(abs(fft(x)))))  ))));

Answer 2

Check out Julius O. Smith III's write up.

There is a Hilbert transform relationship between the magnitude response, $G(\omega)$, and the phase response, $\theta(\omega),$ of the associated minimum phase filter.

If

$$ H(\omega) = G(\omega) \exp(\jmath \theta(\omega)) $$ then $$ \ln( H(\omega) ) = \ln(G(\omega)) + j\theta(\omega) $$ and $$ \theta(\omega) = - {\cal H}\left[ \ln(G(\omega)) \right] $$ where ${\cal H}$ is the Hilbert transform.

Also check out Wikipedia's write up.