Continuous phase for phase delay calculus in FIR filters

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.

• You can always interpolate more intermediate points (such as by zero-padding before an FFT) to help estimate the unwrapped/continuous phase. Just keep increasing the density of interpolated points until the "jumps" in phase become below whatever bound you think you require to estimate continuous phase. – hotpaw2 Apr 4 '16 at 17:53
• Well the zero padding doesn't really help, because the size on the y/phase axis of the jump is not depending on the FFT size... – IvanC Apr 4 '16 at 18:26
• the phase delay for the linear phase filter should be a constant. that means that the continuous phase has to be a straight line, without any phase jump... I might need something else that the wrapped/unwrapped phase as a basis to get the continuous phase – IvanC Apr 4 '16 at 18:51

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()
• It is true that interpolating noise results in noise. However the sign changes in phase responses for anything above the noise floor are not random. They are usually the result of interpolating incorrectly. – hotpaw2 Apr 4 '16 at 19:57
• @hotpaw2 I don't believe you are correct that it's about interpolating incorrectly. Can you check my additions and let me know what I'm doing wrong? – Peter K. Apr 4 '16 at 21:44
• If the FFT phase alternates in sign when the magnitude is above the noise floor, this is often because the FFT phase is being referenced to a discontinuity. This artifact can often be removed for interpolation purposes by a partial fftshift to the center of an impulse response (e.g. rotate the zero-padded FFT input until the mid or modal high point of the FIR is in the center of the FFT input vector, the result can be un-fft-shifted after the FFT and phase interpolation). Below the noise floor, the sign will not simply alternate, but become more random. – hotpaw2 Apr 4 '16 at 22:10
• Correction: The partial or full fftshift should center the FIR around the first element of the vector fed to the FFT. – hotpaw2 Apr 4 '16 at 22:38

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.

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 (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 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.

• OK. That seems somewhat better; I've updated my answer. Still doesn't show the "ideal" phase, though... Any further thoughts? – Peter K. Apr 4 '16 at 23:32
• It's already done at 200,000 points for a 125 length FIR filter... – Peter K. Apr 5 '16 at 0:03
• Isn't that what complex_magnitude <- polyval(f1[c(63:125,1:62)], exp(-1i*freq)) does ? The indices of f1 should do the same as fftshift. R's start index is 1. – Peter K. Apr 5 '16 at 0:09
• Changing indexes is not the same as an fftshift when using polyval: the phase relationship will be different with respect to all the non-orthogonal-in-aperture exponentials. A shift into negative time (rotation across the 1st element) will work better for all exponentials with phase 0 at t=0. Or try an fftshift with a regular FFT of length (N * padding_ratio). – hotpaw2 Apr 5 '16 at 0:13
• I'm confused (probably too late and too distracted); what else is fftshift BUT changing indices? I'll regroup in the morning and see if things are clearer then. Thanks for your engagement! – Peter K. Apr 5 '16 at 0:19