# Derive minimum phase from magnitude

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
*/

-
Are you interested only in obtaining real coefficients for your filter when you apply the IDFT, or is there a significant reason for going after "minimal phase"? –  user2718 Feb 17 '13 at 14:59
it is not clear to me the condition that permits to say that from $$\ln( H(\omega) ) = \ln(G(\omega)) + j\theta(\omega)$$ in case of minimum phase system it can be seen that $$\theta(\omega) = - {\cal H}\left[ \ln(G(\omega)) \right]$$ Could someone please explain in more details? The sentence "Then the phase response $\Theta(\omega)$ can be computed as the Hilbert transform of $\ln G(\omega)$. This can be seen by inspecting the log frequency response" in Julius O. Smith is not clear to me. Why inspecting the ln() of the frequency response should I understand that?. –  user5264 Aug 19 '13 at 17:11

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)))))  ))));

-
Thanks. I marked this one as the answer because it included code and pointed toward IEEE article. The ultimate answer was found in an IEEE paper by Venkata, Evans, and McCaslin: hpl.hp.com/personal/Niranjan_Damera-Venkata/files/…;. –  Speed Coder Feb 26 '13 at 15:14

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.

-
Thanks. Vote up because this was mathematically the correct direction. –  Speed Coder Feb 26 '13 at 15:07