0
$\begingroup$

My problem is pretty simple, I've designed a magnitude response and I would like to find the corresponding minimum phase filter. I'm using the code below and unless there is a bug my eyes don't want to see, I'm pretty confident on the methodology.

nbTaps = 100;
fs = 96000;
desiredMagn_dB = [0 0 -100 -100];
f = [0 24000 26000 fs/2];

FIRFreq = (0:nbTaps/2)*fs / nbTaps;

FIRMagn = interp1(f, 10.^(desiredMagn_dB./20), FIRFreq);
FIRMagn = [FIRMagn fliplr(FIRMagn(2:end-1))];

minPhase = -imag(hilbert(log(FIRMagn)));

freqResp = FIRMagn.*exp(1i.*minPhase);

FIRCoefs = ifft(freqResp, 'symmetric');

And yet, my filter has not a minimal phase, since I have zeros outside the unit circle. I am a little bit lost, here, and I would like to understand why is so? or where I am wrong?

EDIT from r b-j:

nbTaps = 16384;
fs = 96000;
desiredMagn_dB = [0 0 -100 -100];
f = [0 24000 26000 fs/2];

FIRFreq = (0:nbTaps/2)*fs / nbTaps;

FIRMagn = interp1(f, 10.^(desiredMagn_dB./20), FIRFreq);
FIRMagn = [FIRMagn fliplr(FIRMagn(2:end-1))];

minPhase = -imag(hilbert(log(FIRMagn)));

freqResp = FIRMagn.*exp(1i.*minPhase);

FIRCoefs = ifft(freqResp, 'symmetric');

figure(1)
plot(20*log10(FIRMagn))
figure(2)
plot(minPhase)
figure(3)
plot(FIRCoefs)

i don't see anything unexpected.

$\endgroup$
  • $\begingroup$ as a quick response, i am not sure about the minus sign on the instruction minPhase = -imag(hilbert(log(FIRMagn))); it appears that Julius Smith has some MATLAB code. i dunno what interp1() does and i would increase nbTaps and i would simply scale your dB representation to nepers and send that directly to the hilbert() operator. $\endgroup$ – robert bristow-johnson May 10 at 18:27
  • $\begingroup$ okay, the minus sign is correct. i just checked. $\endgroup$ – robert bristow-johnson May 10 at 19:01
  • $\begingroup$ Julie, i just ran your script with nbTaps=1024; i don't see what's wrong. $\endgroup$ – robert bristow-johnson May 10 at 19:07
  • $\begingroup$ Thank you Robert for your answer. Two things I don't understand. First, to be sure I have a minimum phase, I run the Matlab function "isminphase" which answers no with 100 taps or 16834. $\endgroup$ – Julie May 13 at 10:23
  • $\begingroup$ Second, if you zoom in FIRCoefs: figure(3) plot(FIRCoefs) axis([1 length(FIRCoefs) -.001 .001]) You will see a "replica" in the middle of the IR. Why? and where is it coming from? I would have expected a smoothed decline until the end of the impulse response. $\endgroup$ – Julie May 13 at 10:28
4
$\begingroup$
  1. More taps. You don't have anywhere near enough taps for a filter that steep. Start large with 8192 or so cut to desired accuracy, if needed
  2. Due to the low number of tabs you are seeing the effect of "circular" hilbert transform. See for example: http://andrewduncan.net/air/
  3. How do you know you have zeros outside the unit circle? Calculating the roots of a high order polynomial is a numerically difficult problem and using a Matlab function like "roots()" will given you often inaccurate or plainly wrong results.
$\endgroup$
  • $\begingroup$ Thank you Hilmar. Yes, I've used the zplane function for Matlab to look at the zeros, which was, for me, a reliable tool until now. I've increased the number of taps to 8192 and I still don't have a minimal phase (I've used 'isminphase'. If this one is not a reliable tool neither, how can I characterize my filter?) I do understand the circular argument but anyway, from my understanding, to each filter having a certain magnitude, there is one filter having a minimal phase. So for any number of taps, I should be able to compute a minimal phase, rigth? $\endgroup$ – Julie May 10 at 9:42
  • $\begingroup$ One more question regarding the circular principle: how predefine the good number of taps to be at the limit of the wrap around due to the circularity? $\endgroup$ – Julie May 10 at 10:29
  • $\begingroup$ Hilmar, that Andrew Duncan paper is excellent. i remember when it first came out. $\endgroup$ – robert bristow-johnson May 10 at 18:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.