Phase plot from impulse response function does not match that of the frequency response function

Question

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.

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'))

Answer 1

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.