Continuous phase for phase delay calculus in FIR filters

Question

I would like to estimate the phase delay accurately for any random FIR filter.

The definition of the phase delay is the continuous phase divided by the angular frequency (with a sign change). That means that you can't simply use the wrapped phase, between $-\pi$ and $+\pi$ you can get with a FFT to obtain the phase delay estimation. At least, you need to use an unwrapped version of the phase.

For example, with a linear phase FIR filter, the phase delay is supposed to be a constant, so the phase you have to use to do the calculus has to be monophonic and decreasing.

However, by doing some tests in MATLAB, I have been able to see that the unwrapped phase doesn't have this property everywhere. For a halfband lowpass filter for example, there are a few "jumps" in the stopband with a length not related with $\pi$, so they are untouched in the unwrapped phase, which gives me a wrong estimation of the phase delay, featuring jumps as well instead of a constant value.

Indeed, the "continuous phase" I need isn't the same thing than the unwrapped phase, and is associated with the zero-phase amplitude of the FIR filter. I have not found yet a way to get it, and MATLAB with its zerophase function is giving wrong results as well, because it cheats ! If the filter is linear phase, it uses the formula we all know for the phase delay. If it is a FIR linear phase transformed into a minimum-phase filter for example, then the phase delay displayed features again this low size jumps.

I'm thinking about doing something like detecting the jump length, assuming the phase is supposed to be smooth everywhere, and using the derivative at the sample before the jump to do something which might look good...

But I would like to know if someone knows a rigorous way to get the continuous phase.

EDIT : here is an example. You can find there the phase for a FIR linear phase filter, wrapped and unwrapped. As you can see on the unwrapped, some phase jumps are still here...

EDIT 2 : my problem can be redefined this way. Most of the time, any FIR filter can be studied with its frequency response H($\omega$) = |H($\omega$)| $e^{j \phi( \omega)}$ with |H($\omega$)| the magnitude, always positive or zero, and $\phi(\omega)$ the wrapped phase.

There is an alternative notation : H($\omega$) = A($\omega$) $e^{j \theta( \omega)}$ with A($\omega$) the amplitude or zero-phase amplitude and $\theta( \omega)$ the continuous phase representation, still between -$\pi$ and $\pi$. Apparently, it features only 2$\pi$ jumps, and can be fully continuous with a standard unwrap algorithm. That's because A($\omega$) can be positive or negative. Each time there is a zero in the FIR filter, A($\omega$) changes its sign. That's the zeroes which are causing the lowest jumps. Since A($\omega$) contains the sign information, the low jumps are no more present in the new phase representation $\theta( \omega)$.

So, what I need now, is in fact a way to calculate accurately $\theta( \omega)$, maybe by calculating first A($\omega$). I might detect the low jumps easily and change the sign of |H($\omega$)|, then use it to get the remaining $e^{j \phi( \omega)}$. But there is maybe a better way of doing this, and I'd like to know if someone familiar with this notation has already done that.

Answer 1

The problem is that magnitudes close to zero in the stop band change their sign almost arbitrarily / randomly. And that's why you're seeing arbitrary jumps of $\pi$ radians --- it's a sign change.

Because of that, I'd be tempted to just ignore the phase of any part of the magnitude response below some low value. What value to use will depend on your application.

Zero padding will not generally help because that just gives you more "random" sign changes around the frequency axis.

I use "random" in quotes as it's not stochastic, but it is hard to predict.

If I interpolate the phase on 100,000 points for an FIR filter of length 126 taps and unwrap the phase I get the picture below.

The reason the unwrapping algorithm fails is because the phase in the stop band is changing very quickly compared with that in the passband. Comparison below of unwrapped phases. The red lines are $\pm \pi/2$.

---

Changing the complex_magnitude expression to mimic the fftshift idea from @hotpaw, to:

    complex_magnitude <- polyval(f1[c(63:125,1:62)], exp(-1i*freq))

I get:

which is slightly better, but not ideal.

---

R Code Below

#29899

library("signal")

f1 <- remez(125, c(0, 0.3, 0.4, 1), c(1, 1, 0, 0))
plot(f1,type="l")
dev.copy(png, 'Q29899-ImpulseResponse.png')
dev.off()

freq <- seq(0,pi,pi/100000)

complex_magnitude <- polyval(f1, exp(-1i*freq))

plot(freq, unwrap(Arg(complex_magnitude)), type="l");
dev.copy(png, 'Q29899-UnwrappedPhase.png')
dev.off()

plot(freq[60000:100000],Arg(complex_magnitude[60000:100000]), type="l")
lines(freq[c(60000,100000)],c(pi/2,pi/2), col="red")
lines(freq[c(60000,100000)],c(-pi/2,-pi/2), col="red")
title('Stop band phase detail')
dev.copy(png, 'Q29899-PhaseDetailWrapped.png')
dev.off()

plot(freq[1:10000],Arg(complex_magnitude[1:10000]), type="l")
lines(freq[c(1,10000)],c(pi/2,pi/2), col="red")
lines(freq[c(1,10000)],c(-pi/2,-pi/2), col="red")
title('Pass band phase detail')
dev.copy(png, 'Q29899-PhaseDetailWrappedPassband.png')
dev.off()

Answer 2

Just found the real answer. So far it's working well for a variety of different FIR filters.

Here is what I did. I detect the zero-crossings of the magnitude, using the derivative of the phase and the derivative of the magnitude for example. Each time I see one, I change the sign of the magnitude. By doing so, I get the amplitude $A(\omega)$ with values positive or negative.

Then, I divide the $H(\omega)$ by the amplitude I have calculated, to get $e^{j\theta(\omega)}$. Finally I take the phase $\theta(\omega)$ from this new term, unwrap the result, and I have a continuous phase without any extra jump, well suited to calculate the phase delay.

The true power of this method is that it works even when the zeroes are close to be zeroes but are not exactly, which might happen for some minimum-phase FIR filters depending on the algorithm used to convert a linear phase FIR filter. In this particular case, the low jumps are not straight $\pi$ jumps, so all the other methods might fail.

By the way, shifting the filter coefficients is totally useless, because it didn't solve at all the problem of the phase delay calculus.

Answer 3

Interpolating more points can assist the phase unwrapping algorithm. To interpolate phase, do a partial negative circular rotation of the input vector to center the FIR filter around x[0] (first element of the vector being fed to an FFT) in a sufficiently zero-padded window, before the FFT.

What this rotation does is to move the edges of the inherent rectangular window around the impulse response away from the phase measurement point (x[0] for a typical DFT or FFT). Otherwise, you might end up interpolating the implied Gibbs horns around the discontinuity at the first element or last element of the DFT input. Thus a small rotation, plus zero-padding to interpolate more points, might remove some of the zig-zag artifact of the unwrapped phase (for those portions of the passband, or even stop-band sufficiently above the noise floor).

The phase reference point can be moved back to any other point in time, if needed, by the linear phase shift identity for time-shifted data vectors (after interpolation and phase unwrapping, of course).

Note that if the spectrum magnitude passes thru zero at a point where it is being sampled, that point should probably be removed before the phase unwrapping algorithm tries to use it.