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]for 0≤k≤N2
G[0] is the magnitude at DC. G[N2] (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])0≤k≤N2
If you started with gain in dB, it has already been logged. but dB are not nepers. you must multiply each value by 120ln(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]N2<k≤N−1
Then, to compute the Hilbert transform, there are 3 steps. First, compute
h[n]=DFT{H[k]}=N−1∑k=0H[k]e−j2πnkN
Then, multiply every positive time index (n<N2) with −j=e−jπ2 (or spin those complex values by -90°) and multiply every negative time index (n>N2) with +j=e+jπ2 (or spin those complex values by +90°). h[0] and h[N2] (if N is even) should be set to 0.
h[n]←{0n=0−j⋅h[n]1≤n<N20n=N2j⋅h[n]N2<n≤N−1
Finally inverse transform that result (and negate)
ϕ[k]=−IDFT{h[n]}=−1NN−1∑n=0h[n]ej2πnkN
ϕ[k] is the phase, in radians, of the minimum-phase system. Your complex transfer function is
G[k]ejϕ[k]
No comments:
Post a Comment