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.
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 &lt;- 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 &lt;- 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 &lt;- seq(0,pi,pi/100000)

complex_magnitude &lt;- 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()