# Gaussian filter: Plotting DTFT and DFT (by hand) from the continuous-time impulsive response

I am trying to make an algorithm that plots out the Discrete-Time Fourier Transform (DTFT) and the Discrete Fourier Transform (DFT) of the Gaussian filter. The impulsive response and its transfer function (TF) are respectively given by $$g(t) = \frac{1}{\sigma\sqrt{2\pi}} e^{-\frac{t^2}{2\sigma^2}}$$ and $$G(f) = \mathscr{F}\big\{g(t)\big\} = \int_{-\infty}^{\infty} g(t) \, e^{-j 2 \pi f t} \, \mathrm{d}t = e^{- \frac{(2\pi\sigma f)^2}{2}}$$

Where $$\mathscr{F}\{\cdot\}$$ indicates the Continuous-Time Fourier Transform (CTFT, or just FT). I have studied DSP for a while, but I am still not sure how to proceed in order to derive (and plot) the DTFT and the DFT from this impulsive response. I would like to avoid built-in command that automatically plots it and compute it by myself. Comparing the obtained plot with these built-in results would be interesting, though.

• Are you sure you want the bottom limit of the integral to be $-j\infty$? That's not compatible with the definition of the continuous fourier transform and the integral becomes a little bit non-sensical. $t$ is purely real. Jun 9, 2023 at 3:20
• @robertbristow-johnson That was a typo Jun 18, 2023 at 22:35
• I thought so. I was tempted to correct it, but I had already been dabbling in your question and I want to leave it as your question. Jun 18, 2023 at 22:38

Before anything, let me just rewrite your TF from linear frequency ($$f$$, in Hertz) to angular frequency ($$\Omega = 2\pi f$$, in rad/sec). The notation $$\Omega$$ is usually adopted in the field of DSP to differ it from the discrete angular frequency ($$\omega = \Omega T_s$$, in rad), where $$T_s$$ (in seconds) is the sampling period. Therefore, the TF is given by $$G(j\Omega) = e^{- \frac{(\sigma \Omega)^2}{2}}$$

To plot the DTFT and the DFT, you must first define the parameters of your problem. The Gaussian filter, as you stated, has $$\sigma$$ as parameter, which is the standard deviation. By setting this to, say, $$\sigma = 0.1$$, we can plot this impulsive response and this FT.

## FT

This is a advisable step as it is sensible to take a look at this spectrum before choosing a sampling rate in order to avoid aliasing. In Matlab, you can express a Dynamic System Model in tons of manners. In your case, you already have the analytical expression $$G(j \Omega)$$, so lets make use of it:

%% Continuous-time plot

% parameters

Np1 = 401; % number of points for the plot of g(t) and |G(f)|,i.e., a total of N+1 samples
t = linspace(-2,2,Np1); % time base (from -2s to 2s)
Ts = t(2)-t(1); % sampling period (10 ms)
fs = 1/Ts; % sampling rate (100 Hz)
Omega = linspace(-pi/Ts,pi/Ts,Np1); % frequency base (from -100π rad/s to 100π rad/s)

sigma=0.1;
variance=sigma^2;

g_t = 1/(sqrt(2*pi*variance))*(exp(-t.^2/(2*variance))); % g(t)
G_Omega = exp(-(sigma*Omega).^2 / 2); % G(jΩ)

% plot
subplot(2,1,1)
plot(t,g_t,'b');
grid on;
title('g(t)');
xlabel('Time(s)');
ylabel('Amplitude');

subplot(2,1,2)
plot(Omega, abs(G_Omega),'r');
grid on;
title('|G(j\Omega)|');
ylabel('|G(j\Omega)|');
xlim([Omega(1), Omega(end)]);


You might want to plot $$|G(j\Omega)|$$ in semilogarithm scale. In this case, assuming that $$g(t)$$ express a root-power quantity, such as a voltage or current in a electronic circuit, $$|G(j \Omega)|$$ in dB is given by: $$|G_{dB}(j \Omega)| = 20 \log(|G(j \Omega)|)$$

Anyway, one can see that that most of the effective bandwidth of $$G(j\Omega)$$ is something between $$-30$$ and $$30$$ rad/sec. By denoting the one-sided bandwidth of $$G(f)$$ as $$\Omega_B = 30$$ rad/sec, from the Nyquist criterion, you must select a sampling rate such that $$\Omega_s = \frac{2\pi}{T_s} \geq 2\Omega_B$$ By adopting $$f_s = 1/T_s = 100 Hz \therefore \Omega_s = 2\pi 100$$ rad/s, we can fulfill the Nyquist requirement.

## DTFT

One way to implement a digital filter is by sampling the impulse response, $$g(t)$$, in order the obtain the same impulse response of the equivalent analog filter. Such an approach is usually called impulse invariance method. However, we usually normalize the signal samples by $$T_s$$ so that we cancel the $$T_s$$ is the DTFT and get the very same transfer function of the analog filter (vide Oppenheim, eq (4.20)). Hence, by defining $$g[n] = T_s g(nT_s)$$ its DTFT is given by $$G(e^{j\omega}) = \sum_{k = - \infty}^{\infty} G\left[j\left(\frac{\omega}{T_s} - \frac{2\pi k}{T_s} \right)\right],$$

In the particular case of no aliasing, we have that

\begin{align} G (e^{j\omega}) = G\left( j\frac{\omega}{T_s}\right) & \textrm{ , for }|\omega| \leq \pi \end{align}

It is obvious that $$G_s (e^{j\omega})$$ is actually defined in $$(-\infty, \infty)$$, but we usually represent this signal only in $$[-\pi, \pi]$$ as it is periodic in $$2\pi$$.

%% DTFT
G_DTFT = G_Omega;
omega = Ts .* Omega;

plot(omega, abs(G_DTFT), 'b');
title('|G(e^{j\omega})|');
ylabel('Amplitude');
grid on;
xlim([-pi,pi]);
xticks([-pi,-pi/2, 0, pi/2, pi]);
xticklabels(["-\pi","-\pi/2", "0", "\pi/2", "\pi"]);


As you can see, plotting the DTFT is ridiculously easy (as long as you ensure that aliasing does not occur!) since it just involves frequency scaling (if you don't normalize the samples by $$T_s$$, it involves a scaling in the amplitude as well).

## DFT

For the DFT, you actually have two way outs. For the first one, note that (vide Oppenhneim, again)

$$\left.G[k] = G(e^{j\omega})\right\vert_{\omega = \frac{2\pi k}{N+1}}$$ where $$N+1$$ is the number of samples of the sequence $$\{G[k]\}$$, which is also equal to the number of samples of the sequence $$\{g[n]\}$$. Since what you have on computer is samples of $$G(e^{j\omega})$$ instead of itself, strictly speaking, you already made it. The plot() function makes the illusion, by interpolating the samples, that you have a continuous curve instead of a vector of values. Since we have plotted $$G(e^{j\omega})$$ from $$-\pi$$ to $$\pi$$, let us define $$k \in \left\{-\frac{N}{2}, -\frac{N}{2}+1, \cdots, \frac{N}{2}-1, \frac{N}{2}\right\}$$.

%% DFT 1th method - by hand
G_DFT = G_DTFT;
k = -(Np1-1)/2:(Np1-1)/2;
stem(k, abs(G_DFT), 'b');
title('|G[k]|');
xlabel('k');
ylabel('Amplitude');
xticks([-(Np1-1)/2, -(Np1-1)/4, 0, (Np1-1)/4, (Np1-1)/2]);
xticklabels(["-N/2", "-N/4", "0", "N/4" "N/2"]);
grid on;


The second way out, which kinds of serves as a comparison with the Matlab built-in function, is directly applying the FFT on the samples $$\{g[n]\}$$ (remember that the FFT is a efficient implementation of the DFT computation, and the final result is exactly the same). The computation of the FFT is given by: $$G[k] = \frac{1}{N+1} \sum_{n = 0}^{N} g[n] W_{N+1}^{nk},$$ where $$W_{N+1} = e^{-j\frac{2\pi}{N+1}}$$ is called twiddle factor.

The magic is that, by computing $$G[k]$$ from the previous equation, it shall exactly be the samples of $$X(e^{j\omega})$$ in the points $$\omega = \frac{2\pi k}{N}$$ as the relationship between DTFT and DFT must match. Just pay attention that your vector of samples must begin with $$g[0]$$. Likewise, the sequence $$\{G[k]\}$$ returned by fft() also begins with $$G[0]$$. Shifting it in order to obtain the very same plot is easily done by fftshift().

%% DFT 2th method - by using the built-in FFT function
g_k = Ts * g_t; % g[k]
G_k = fft(fftshift(g_k)); % G[k] (the DFT of g[k])

stem(k, abs(fftshift(G_k)), 'b');
title('|G[k]|');
xlabel('k');
ylabel('Amplitude');
xticks([-(Np1-1)/2, -(Np1-1)/4, 0, (Np1-1)/4, (Np1-1)/2]);
xticklabels(["-N/2", "-N/4", "0", "N/4" "N/2"]);
ylim([0, 1]);
grid on;


• In your explanation, what your call the total number of samples, $N+1$, other authors call it just $N$. That is right? Why did you choose to denote it differently? May 17, 2023 at 18:23
• Yes, that is right. I denoted it differently just to make the number of samples around $0$ symmetric. May 17, 2023 at 18:25