I can't find out if it possible to compute the minimum-phase response corresponding to a given magnitude response using a Hilbert transformer. Is that possible?
When I write Hilbert transformer I mean a 90-degree phase shifter.
I know other ways to compute the minimum-phase response but since there are IIR filters that approximately can realize a Hilbert transformer I was wondering if it is possible to use the Hilbert transformer. Not sure if the answer is obvious but it is not a homework question.
Edit:
Implementation of proposed
function y = test_minph(Mag)
Mag = Mag(:);
x = [Mag; Mag(end-1:-1:2)];
len = length(x);
N = (len)/2-1;
wn = [0; -1i*ones(N,1); 0; 1i*ones(N,1)];
xhat = real(fft(log(x)));
y = -ifft(wn.*xhat);
end
But the question is about how to compute the below using a Hilbert transformer (if possible) which is what robert johnson proposed.
function y = minphase(Mag)
Mag = Mag(:);
x = [Mag; Mag(end-1:-1:2)];
len = length(x);
N = (len)/2-1;
wn = [1; 2*ones(N,1) ; 1; zeros(N,1)];
xhat = real(ifft(log(abs(x))));
y = imag(exp(fft(wn.*xhat)));
end
So it seems there are two different Hilbert transforms in play (don't know if they are dual) and I'm not sure how to compute the minimum-phase Hilbert transform using the 90 degree phase shift Hilbert transform. I hope it makes sense.
Answer
as a related aside question i posted this question about minimum-phase filters and the phase-magnitude relationship.
let $N$ be the FFT size you will use. (often $N$ is a power of two, but it doesn't have to be.)
the target magnitude response is
$$ G[k] \qquad \text{for } 0 \le k \le \tfrac{N}{2} $$
$G[0]$ is the magnitude at DC. $G[\tfrac{N}2]$ (if $N$ is even) is the magnitude at Nyquist. $G[k]$ must be purely real and positive. No complex values and no polarity (sign) changes. Then, when converting to log-magnitude, the logarithm function should not give you trouble.
First thing is convert magnitude to nepers using the natural (base-$e$) logarithm.
$$ H[k] = \ln(G[k]) \qquad 0 \le k \le \tfrac{N}2 $$
If you started with gain in dB, it has already been logged. but dB are not nepers. you must multiply each value by $\frac{1}{20}\ln(10)$ = 0.115129255 (that's one of these magic numbers we get to see in DSP).
you need not worry about any constant added to the log-magnitude (which would correspond to a positive gain factor in $G[k]$). the Hilbert transform of a constant is zero so the minimum-phase result will be unchanged no matter how $G[k]$ is scaled.
You must mirror the first half into the latter half (the latter half of the DFT corresponds to negative frequencies or negative times):
$$ H[k] = H[N-k] \qquad \tfrac{N}2 < k \le N-1 $$
Then, to compute the Hilbert transform, there are 3 steps. First, compute
$$\begin{align} h[n] &= \mathcal{DFT}\bigg\{ H[k] \bigg\} \\ \\ &= \sum\limits_{k=0}^{N-1} H[k] e^{-j 2 \pi \frac{nk}N} \\ \end{align}$$
Then, multiply every positive time index ($n<\tfrac{N}2$) with $-j = e^{-j\frac\pi2}$ (or spin those complex values by -90°) and multiply every negative time index ($n>\tfrac{N}2$) with $+j = e^{+j\frac\pi2}$ (or spin those complex values by +90°). $h[0]$ and $h[\tfrac{N}2]$ (if $N$ is even) should be set to 0.
$$ h[n] \leftarrow \begin{cases} 0 & n=0 \\ -j \cdot h[n] \qquad & 1 \le n< \tfrac{N}2 \\ 0 & n=\tfrac{N}2 \\ j \cdot h[n] \qquad & \tfrac{N}2 < n \le N-1 \\ \end{cases} $$
Finally inverse transform that result (and negate)
$$\begin{align} \phi[k] &= -\mathcal{IDFT}\bigg\{ h[n] \bigg\} \\ \\ &= - \tfrac1N \sum\limits_{n=0}^{N-1} h[n] e^{j 2 \pi \frac{nk}N} \\ \end{align}$$
$\phi[k]$ is the phase, in radians, of the minimum-phase system. Your complex transfer function is
$$ G[k] \, e^{j \phi[k]} $$.
No comments:
Post a Comment