FIR filter phase response at fs/2

Question

I am currently building a FIR filter based on frequency sweeping data (0 Hz to Fs/2), i.e., transfer curve in frequency with the sampling rate of Fs.

1. The phase response is shown in the following figure. From the measurement results, I saw that the phase did not end at plus/minus multiple of $\pi$. As I know, a normal FIR filter always has it phase end at plus/minus multiple of $\pi$. Does anyone know how to explain this issues or it has problem with my measured response?

2. With that phase response, I used sampling method (iFFT) to generate the time domain coefficients. The result show a large deviation near the main tap as shown in the figure. If I manipulate the phase to let it ends at multiple of pi by adding a fractional delay to the phase response, I will get a very beautify coefficients, with a few taps only (other taps are very small). However, I need to face with fractional delay issues. On the other hand, I can shift, and truncate the coefficients to have smaller number of taps. But the side effect is that there is a drop in amplitude near the Nyquist frequency. And this creates severe disadvantage to my work. How to solve this problem?

Answer 1

If I use the scipy.signal.firls function to design a filter:

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

N = 64
frequencies = np.linspace(0,1,N)
phase = np.linspace(0,-0.6,N)

plt.figure(1)
plt.plot(phase)

h =  signal.firls(9, frequencies, np.exp(1j*phase) )

plt.figure(2)
plt.plot(np.real(h))
plt.plot(np.imag(h))

w, resp = signal.freqz(h)

plt.figure(3)
plt.plot(w,np.angle(resp*np.exp(1j*w*4)))

Where I start with this:

then I get a complex-valued impulse response:

where the blue plot is the real component and the orange plot is the imaginary component.

Then, modulo the delay of 4 samples, this gives a pretty good approximation to the original requirement:

The magnitude response looks OK too:

---

So the problem seems to be that you need a fractional delay with a good response at $f_s/2$, whereas most fractional delay FIR filters end up having a zero magnitude response there.

If I do a standard$^1$ sinc design for a fractional delay filter:

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

HalfT = 32
T = 2*HalfT + 1

t = np.arange(T)-T/2
delay = 0.6/np.pi

h_fractional = np.sin(2*np.pi*0.5*(t-delay))/(np.pi*(t-delay)) #*np.power(-1,t)


w, resp = signal.freqz(h_fractional)

plt.figure(1)
plt.plot(h_fractional)

plt.figure(2)
plt.plot(w,np.angle(resp*np.exp(1j*w*T/2)))
plt.plot(w,np.linspace(0,-0.6,len(w)))
#plt.ylim([-0.8,0])

plt.figure(3)
plt.plot(w,np.abs(resp*np.exp(1j*w*T/2)))

Then I get the following phase and magnitude responses.

Some things I'd suggest:

• Oversample: if you really need the response to be good close to your sampling frequency, you're sampling too slow to start with.

• After oversampling, design a bandpass fractional delay filter. That should avoid the complex-valued problem.

$^1$ $\color{red}{\mbox{Don't do this!}}$ Use these tools instead.

Answer 2

My answer to #1:

A normal FIR filter does not necessarily have its phase response ending at plus or minus a multiple of $\pi$. A linear phase FIR filter specifically would have its phase ending at a multiple of $-\pi/2$ (or positive if the phase isn't unwrapped). This is explained as follows:

As explained here, the delay in samples for linear phase FIR filters is $(N-1)/2$ where $N$ is the number of coefficients in the filter.

One sample delay has a phase response that extends from $0$ at DC to $\pi$ at $f_s/2$ where $f_s$ is the sampling rate. This is consistent with the Fourier Transform of a sample delay $T=1/f_s$ as:

$$\delta(t-T) \rightarrow e^{-j 2 \pi f T}$$

And more directly, the z-Transform of a sample delay where we see the same result:

$$z^{-1} \rightarrow e^{-j2\pi f/f_s}$$

Therefore the phase for any linear phase filter of $N$ taps will go from $0$ to $-(N-1)\pi/2$. Thus for an odd number of taps the phase will "end" at $f_s/2$ at a multiple of $-\pi$, and for ean even number of taps it will be a multiple of $-\pi/2$.

The completely linear phase in the result can only occur with complex conjugate symmetric or antisymmetric coefficients, which would all have the resulting linear phase slope (delay) as I outlined above. Therefore, there appears to be an error in the measurement of what represents phase or the frequency axis is not extending to $f_s/2$.

I note that if the output of the filter was decimated by $D$, the resulting phase slope would also be reduced by a factor $D$. If the OP is measuring the phase response by comparing the input and output, then this could also be an explanation, except in that case we would see the phase extend to $-(N-1)\pi/(2D)$ with integer $D$, which also does not appear to be the case in the result.

The other possibility is a parasitic time delay in the measurement; if we had the actual coefficients and specific details of the measurement approach, we could determine what the error is versus expected which may further refine the different theories.

As for Answer #2, Peter has provided a very good answer and received my upvote - the OP's result is a good example of why not to use the IFFT to determine the impulse response (this is the frequency sampling method of frequency design which is most intuitive but worst choice in most applications for estimating the impulse response for the filter: it will have large error deviations vs other approaches except for the specific frequencies sampled, as I detail further at this post). What Peter chose (least squares) is my favorite approach for an optimum estimate in the least squares sense.

Answer 3

Assuming that the measurements are done correctly, that's most likely aliasing.

In order to sample without aliasing, the energy at and above the Nyquist frequency should be 0. That's probably not the case here, so the periodic repetition of the sampling process creates a discontinuity at Nyquist which generates the time domain ringing.