2
$\begingroup$

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?

enter image description here

  1. 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?

Time-domain coefficients

$\endgroup$
4
  • $\begingroup$ Did you apply fftshift(), or are you using the impulse response exactly as seen in the picture? $\endgroup$ Jul 6, 2022 at 19:42
  • $\begingroup$ Thanks. Yes, I did apply shift to have causal coefficients. $\endgroup$
    – user190055
    Jul 7, 2022 at 15:36
  • $\begingroup$ Without seeing what you actually did (some piece of code), I doubt anyone can answer. What you're saying is very vague (doctor, my foot hurts, what's wrong?). What is the magnitude after which you did the IFFT? How did you apply it: with a built-in function (what program), or manual code (what code)? The 2nd picture doesn't seem to show a symmetric response, so what else is going on in there? $\endgroup$ Jul 7, 2022 at 17:21
  • $\begingroup$ Would you be able to provide your filter coefficients? From that we can determine from that what the expected phase should be. $\endgroup$ Jul 8, 2022 at 23:34

3 Answers 3

3
$\begingroup$

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:

Required phase response

then I get a complex-valued impulse response:

Complex 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:

Phase of complex filter.

The magnitude response looks OK too:

Magnitude response of complex filter.


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.

Phase response of sinc fractional delay filter

Magnitude response of sinc fractional delay filter

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.

$\endgroup$
7
  • $\begingroup$ Thank you very much for your answer. From your script, we've got a compex time-doamin response h. Can we generate real time-domain coefficients using LS method? $\endgroup$
    – user190055
    Jul 9, 2022 at 14:29
  • 1
    $\begingroup$ @user190055 I don't have a good handle on whether firls can be constrained to give real-only coefficients. I suspect if it were to be so constrained that it probably won't give you an answer you like any better. See my update for other suggestions. I'm still not quite clear on what you need or what problem you're trying to solve. $\endgroup$
    – Peter K.
    Jul 9, 2022 at 16:41
  • 2
    $\begingroup$ My favorite approach to doing fractional delay filters (not sure if its best but it's simple straight forward and accurate) is to create the target low pass filter at the interpolated rate using least squares and then decimate that resulting filter into it's polyphase components- each ends up being an all-pass at the fractional delay based on interpolation value chosen (I recently implemented this for fractional correction down to 1/16,384th of a sample!) I then do a binary search for timing correction such that I can reach this accuracy in 14 iterations (when the correction is unknown) $\endgroup$ Jul 9, 2022 at 17:08
  • $\begingroup$ @Peter K. Thank you so much for the suggestion. It's useful to me. My problem is that I am trying to construct a circuit behavior model based on the measurement results, (i.e., sweeping frequency) of a real circuit. I also have a bit concern of how good we can constraint the error to achieve the target response when using the Least-square method. Btw, do you have any Matlab script that implement the above updated part? Thanks. $\endgroup$
    – user190055
    Jul 10, 2022 at 1:26
  • $\begingroup$ @Dan Boschen Thank you. If we use the fractional delay (FD) filter, the FD filter needs to maintain the target amplitude and constant group delay (linear phase) up to FS/2 with limited number of taps (for the sake of low complexity). I think it could be challenging and some tradeoffs may be considered. Other than Least-square method, would you think either Lagrange interpolator or Thrian formular based all-pass filter is another great option? $\endgroup$
    – user190055
    Jul 10, 2022 at 1:31
2
$\begingroup$

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.

$\endgroup$
5
  • $\begingroup$ Thank you so much for the great expalination. I did meadure the phase response by comparing between input and output. I will check that any possible errors indeuced by measurement procedure. $\endgroup$
    – user190055
    Jul 9, 2022 at 14:40
  • $\begingroup$ However, can we explain the ending at a value different from nπ by a fractional delay caused by instrinsic properties of the circuit; such as the electromagnetic wave propagating through the physical length may depend on the parasitic property of transmssion line, any change in the length will lead to a large fractional delay? $\endgroup$
    – user190055
    Jul 9, 2022 at 14:52
  • $\begingroup$ Yes absolutely any delay of the signal in time will result in the slope in phase as given by the Fourier Transform relationship I provided $\endgroup$ Jul 9, 2022 at 15:40
  • $\begingroup$ Thanks. Could you please suggest to me some documents, links to study and implement least squares method for filters design? $\endgroup$
    – user190055
    Jul 9, 2022 at 15:58
  • $\begingroup$ @user190055 Yes I teach an online course on it that will be offered again this October. It's title is "DSP for Wireless Communications" and it will soon be posted for registration at this site: ieeeboston.org/2022-courses $\endgroup$ Jul 9, 2022 at 16:33
1
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ Got it. Thank you for your comment. $\endgroup$
    – user190055
    Jul 12, 2022 at 13:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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