0
$\begingroup$

I have a first order linear filter $$ H(s) = \frac{1}{1+\tau\,s} $$

With the impulse response function, found using the inverse Laplace transform as $$ h(t) = \frac{1}{\tau}\,\exp\left({-\frac{t}{\tau}}\right) $$

I'm trying to match the output from the frequency response to the impulse response of the filter. For that I'm applying the discrete Fourier transform on the impulse response, evaluating $t$ as $t/f_s$. However, the phase I'm getting from the impulse response is significantly different from that of the frequency response (figures below). Any idea why is this happening and how to get an identical magnitude and phase from the impulse response function.

magnitude plot of the filter

phase plot of the filter

To reproduce the figures in R:


library(pracma)

# First order filter definition
H1st <- function(w, tau, descritize = FALSE, Ts) {
    s <- 1i*w
    z  <- exp(1i*w*Ts)

    if (descritize) {
        s <- 2/Ts * (z-1)/(z+1)
    }

    1 /(1 + tau * s)
}

# impulse response function: inverse laplace transform of the function H1st
h1st <- function(t, tau) {
    exp(-t/tau) / tau
}

sampling_freq      <- 20
N                  <- 2000
f_seq              <- seq(0, sampling_freq/2, length.out = N)
omega_seq          <- 2*pi*f_seq
Ts                 <- 1/sampling_freq
tau1               <- 0.8



h1st_calc <- h1st(1:(N*2)*Ts, tau = tau1)

# Normalize
h1st_calc <- h1st_calc/sum(h1st_calc)

# Magnitude plot
plot((abs(H1st(omega_seq, tau = tau1, 
            descritize = FALSE,  Ts = Ts))) ~ f_seq, ty ='l', log = 'xy', 
     ylab = 'magnitude, |H|', xlab = 'frequency [Hz]')
lines((abs(H1st(omega_seq, tau = tau1, 
            descritize = TRUE,  Ts = Ts))) ~ f_seq,
      ty ='l', lty = 2, col = 'blue')
lines(abs(fft(h1st_calc)[1:N]) ~ f_seq, ty ='l', col = 'red')
legend("bottomleft", c("evaluate s directly",
                     'bilinear transform', 'abs of fft of h(t)'), lty = 1,
       col = c('black', 'blue', 'red'))

# Evaluate phase
plot((angle(H1st(omega_seq, tau = tau1, 
            descritize = FALSE,  Ts = Ts))) ~ f_seq, ty ='l', log = 'x', 
     ylab = 'angle', xlab = 'frequency [Hz]')
lines((angle(H1st(omega_seq, tau = tau1, 
            descritize = TRUE,  Ts = Ts))) ~ f_seq,
      ty ='l', lty = 2, col = 'blue')
lines(angle(fft(h1st_calc))[1:(N)] ~ f_seq, ty ='l', log = 'x', col = 'red')
legend("bottomleft", c("evaluate s directly",
                     'bilinear transform', 'angle of fft of h(t)'), lty = 1,
       col = c('black', 'blue', 'red'))


$\endgroup$

1 Answer 1

2
$\begingroup$

I'm applying the discrete Fourier transform

That implies that the signal was sampled and it was sampled with a finite sample rate. Since your signal is not band-limited, this will result in aliasing, which is exactly what you are seeing.

Specifically, sampling in the time domain makes the spectrum periodic and the spectrum must be real at the Nyquist frequency. Hence the phase at Nyquist must be an integer multiple of $\pi$ and the magnitude around Nyquist tend to "flatten out" so the periodic repetition is reasonably continuous.

Additional questions in the comments

could you clarify why must the spectrum be real at Nyquist

That's simply a consequence of the math of the DFT

$$X(k) = \sum_{n=0}^{N-1} x[n] e^{-j2\pi\frac{nk}{N}}$$

Nyquist is at $k = N/2$ so we get

$$X(N/2) = \sum_{n=0}^{N-1} x[n] e^{-j2\pi\frac{nN/2}{N}} = \sum_{n=0}^{N-1} x[n] e^{-j\pi n} = \sum_{n=0}^{N-1} x[n] (-1)^n $$

which is obviously real.

Physically speaking that's a consequence of the spectrum being periodic. Periodicity means that $X(-N/2) = X(N/2)$. For a real input the DFT is conjugate symmetric. A number can only be it's own conjugate symmetric if the imaginary part is zero.

and is there a way we can minimize this aliasing effect if

It's not trivial:

  1. Make the sample rate "high enough" so that the problems occur outside your frequency area of interest
  2. Try fitting a discrete FIR or IIR to the analog target using a least square error method.
  3. Often you can better result doing an "in-between" of sampling the impulse response and the bilinear transform. In this case you can move the zero from $z = -1$ to something a little higher $z = .9$ or or thereabouts

The best method depends on the details and requirements of your specific application.

$\endgroup$
4
  • $\begingroup$ Thank you, could you clarify why must the spectrum be real at Nyquist? and is there a way we can minimize this aliasing effect if we want to use the impulse response function e. g. for convolving a signal. $\endgroup$
    – Crataegus
    Oct 29, 2021 at 14:16
  • 1
    $\begingroup$ @Crataegus: You don't need to use the sampled impulse response if you want to filter a signal. You could just use the discrete-time filter obtained via the bilinear transform. Note that that filter is an IIR filter, unlike the FIR filter given by the sampled impulse response. $\endgroup$
    – Matt L.
    Oct 30, 2021 at 10:33
  • $\begingroup$ @MattL. you certainly can, but that doesn't match the analog filter either $\endgroup$
    – Hilmar
    Oct 30, 2021 at 13:03
  • 1
    $\begingroup$ @Hilmar: No, but the phase is much closer to the analog filter's phase. Note that the phase at Nyquist doesn't need to be an integer multiple of $\pi$ in discrete time if the transfer function has a zero at Nyquist, as is the case with the filter design via the bilinear transform. At Nyquist the phase jumps from $-\pi/2$ to $\pi/2$ and for this reason it approximates the analog filter's phase pretty closely, even close to Nyquist. $\endgroup$
    – Matt L.
    Oct 30, 2021 at 19:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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