Gaussian Pulse – FFT & PSD in Matlab & Python

by

Key focus: Know how to generate a gaussian pulse, compute its Fourier Transform using FFT and power spectral density (PSD) in Matlab & Python.

Numerous texts are available to explain the basics of Discrete Fourier Transform and its very efficient implementation – Fast Fourier Transform (FFT).  Often we are confronted with the need to generate simple, standard signals (sine, cosine, Gaussian pulse, squarewave, isolated rectangular pulse, exponential decay, chirp signal) for simulation purpose. I intend to show (in a series of articles) how these basic signals can be generated in Matlab and how to represent them in frequency domain using FFT.

This article is part of the following books
● Digital Modulations using Matlab : Build Simulation Models from Scratch, ISBN: 978-1521493885
● Digital Modulations using Python ISBN: 978-1712321638
● Wireless communication systems in Matlab ISBN: 979-8648350779
All books available in ebook (PDF) and Paperback formats

Gaussian Pulse : Mathematical description:

In digital communications, Gaussian Filters are employed in Gaussian Minimum Shift Keying – GMSK (used in GSM technology) and Gaussian Frequency Shift Keying (GFSK). Two dimensional Gaussian Filters are used in Image processing to produce Gaussian blurs. The impulse response of a Gaussian Filter is Gaussian. Gaussian Filters give no overshoot with minimal rise and fall time when excited with a step function. Gaussian Filter has minimum group delay. The impulse response of a Gaussian Filter is written as a Gaussian Function as follows

g(t) = \displaystyle{\frac{1}{\sqrt{2 \pi } \sigma} e^{- \frac{t^2}{2 \sigma^2}}}

The Fourier Transform of a Gaussian pulse preserves its shape.

\begin{aligned} G(f) &=F[g(t)]\\ &= \int_{-\infty }^{\infty } g(t)e^{-j2\pi ft}\, dt\\ &= \frac{1}{\sigma \sqrt{2 \pi } } \int_{-\infty }^{\infty } e^{- \frac{t^2}{2 \sigma^2}}e^{-j2\pi ft}\, dt\\ &=\frac{1}{\sigma \sqrt{2 \pi } } \int_{-\infty }^{\infty } e^{- \frac{1}{2 \sigma^2}\left[t^2+j4 \pi \sigma^2 ft \right]}\, dt\\ &=\frac{1}{\sigma \sqrt{2 \pi } } \int_{-\infty }^{\infty } e^{- \frac{1}{2 \sigma^2}\left[t^2+j4 \pi \sigma^2 ft + (j 2 \pi \sigma^2 f)^2 - (j 2 \pi \sigma^2 f)^2\right]}\, dt\\ &=e^{ \frac{1}{2 \sigma^2}(j 2 \pi \sigma^2 f)^2}\frac{1}{\sigma \sqrt{2 \pi } } \int_{-\infty }^{\infty } e^{- \frac{1}{2 \sigma^2}\left[t+j 2 \pi \sigma^2 f \right]^2}\, dt\\ &=e^{ \frac{1}{2 \sigma^2}(j 2 \pi \sigma^2 f)^2}=e^{ - \frac{1}{2}( 2 \pi \sigma f)^2}\end{aligned}

The above derivation makes use of the following result from complex analysis theory and the property of Gaussian function – total area under Gaussian function integrates to 1.

By change of variable, let (u=t+j 2 \pi \sigma^2 f ). 

\displaystyle{\frac{1}{\sigma \sqrt{2 \pi} }\int_{-\infty }^{\infty }e^{- \frac{1}{2 \sigma^2}\left[t+j 2 \pi \sigma^2 f \right]^2}\, dt =\frac{1}{\sigma \sqrt{2 \pi } }\int_{-\infty }^{\infty }e^{- \frac{1}{2 \sigma^2} u^2}\, du =1}

Thus, the Fourier Transform of a Gaussian pulse is a Gaussian Pulse.

\displaystyle{\frac{1}{\sqrt{2 \pi } \sigma} e^{- \frac{t^2}{2 \sigma^2}} \rightleftharpoons e^{ - \frac{1}{2}( 2 \pi \sigma f)^2}}

Gaussian Pulse – Fourier Transform using FFT (Matlab & Python):

The following code generates a Gaussian Pulse with ( \sigma=0.1s ). The Discrete Fourier Transform of this digitized version of Gaussian Pulse is plotted with the help of (FFT) function in Matlab.

For Python code, please refer the book Digital Modulations using Python 

fs=80; %sampling frequency
sigma=0.1;
t=-0.5:1/fs:0.5; %time base

variance=sigma^2;
x=1/(sqrt(2*pi*variance))*(exp(-t.^2/(2*variance)));
subplot(2,1,1)
plot(t,x,'b');
title(['Gaussian Pulse \sigma=', num2str(sigma),'s']);
xlabel('Time(s)');
ylabel('Amplitude');

L=length(x);
NFFT = 1024;
X = fftshift(fft(x,NFFT));
Pxx=X.*conj(X)/(NFFT*NFFT); %computing power with proper scaling
f = fs*(-NFFT/2:NFFT/2-1)/NFFT; %Frequency Vector

subplot(2,1,2)
plot(f,abs(X)/fs,'r');
title('Magnitude of FFT');
xlabel('Frequency (Hz)')
ylabel('Magnitude |X(f)|');
xlim([-10 10])
Figure 1: Gaussian pulse and its FFT (magnitude)
Figure 1: Gaussian pulse and its FFT (magnitude)

Double Sided and Single Power Spectral Density using FFT:

Next, the Power Spectral Density (PSD) of the Gaussian pulse is constructed using the FFT. PSD describes the power contained at each frequency component of the given signal. Double Sided power spectral density is plotted first, followed by single sided power spectral density plot (retaining only the positive frequency side of the spectrum).

Pxx=X.*conj(X)/(L*L); %computing power with proper scaling
figure;
plot(f,10*log10(Pxx),'r');
title('Double Sided - Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('Power Spectral Density- P_{xx} dB/Hz');
Gaussian Pulse Double Sided Power Spectral Density
Figure 2: Double sided power spectral density of Gaussian pulse
X = fft(x,NFFT);
X = X(1:NFFT/2+1);%Throw the samples after NFFT/2 for single sided plot
Pxx=X.*conj(X)/(L*L);
f = fs*(0:NFFT/2)/NFFT; %Frequency Vector
figure;
plot(f,10*log10(Pxx),'r');
title('Single Sided - Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('Power Spectral Density- P_{xx} dB/Hz');
Gaussian Pulse Single Sided Power Spectral Density
Figure 3: One-sided power spectral density of Gaussian pulse

For Python code, please refer the book Digital Modulations using Python

Rate this article: PoorBelow averageAverageGoodExcellent (14 votes, average: 3.71 out of 5)

Topics in this chapter

Essentials of Signal Processing
● Generating standard test signals
 □ Sinusoidal signals
 □ Square wave
 □ Rectangular pulse
 □ Gaussian pulse
 □ Chirp signal
● Interpreting FFT results - complex DFT, frequency bins and FFTShift
 □ Real and complex DFT
 □ Fast Fourier Transform (FFT)
 □ Interpreting the FFT results
 □ FFTShift
 □ IFFTShift
● Obtaining magnitude and phase information from FFT
 □ Discrete-time domain representation
 □ Representing the signal in frequency domain using FFT
 □ Reconstructing the time domain signal from the frequency domain samples
● Power spectral density
● Power and energy of a signal
 □ Energy of a signal
 □ Power of a signal
 □ Classification of signals
 □ Computation of power of a signal - simulation and verification
● Polynomials, convolution and Toeplitz matrices
 □ Polynomial functions
 □ Representing single variable polynomial functions
 □ Multiplication of polynomials and linear convolution
 □ Toeplitz matrix and convolution
● Methods to compute convolution
 □ Method 1: Brute-force method
 □ Method 2: Using Toeplitz matrix
 □ Method 3: Using FFT to compute convolution
 □ Miscellaneous methods
● Analytic signal and its applications
 □ Analytic signal and Fourier transform
 □ Extracting instantaneous amplitude, phase, frequency
 □ Phase demodulation using Hilbert transform
● Choosing a filter : FIR or IIR : understanding the design perspective
 □ Design specification
 □ General considerations in design

Books by the author

Wireless Communication Systems in Matlab
Wireless Communication Systems in Matlab
Second Edition(PDF)
PoorBelow averageAverageGoodExcellent (180 votes, average: 3.62 out of 5)

Digital modulations using Python
Digital Modulations using Python
(PDF ebook)
PoorBelow averageAverageGoodExcellent (134 votes, average: 3.56 out of 5)

digital_modulations_using_matlab_book_cover
Digital Modulations using Matlab
(PDF ebook)
PoorBelow averageAverageGoodExcellent (136 votes, average: 3.63 out of 5)
Hand-picked Best books on Communication Engineering
Best books on Signal Processing

16 thoughts on “Gaussian Pulse – FFT & PSD in Matlab & Python”

  3. Sorry, my previous suggestion is wrong.
    It should be:
    “plot(f,abs(X)/fs,’r’);”
    Then G(0)=1 for all sigmas and all sample frequencies.
    Peter

    Reply

  4. Hi Mathuranathan,
    I think the line “plot(f,abs(X)/(L),’r’);” in the code is not correct. According to the equations, at the frequency f=0 the magnitude of G(0) is equal to 1. But the result of the code is 0.9877.
    I think the correct line is:
    “plot(f,abs(X)/(L-1),’r’);”
    Note: (L) -> (L-1)

    Kind regards from Germany,
    Peter

    Reply

  5. Hi, thanks again for your works here.
    I have a problem, I’m plotting a 3D Gaussian pulse.
    I have been able to create it and to make the fft of it, but the phase it’s not null and quite incomprehensible.
    This is a quick example of what I get:

    close all;
    clear;
    clc;
    N = 64; % Size of a sunction
    D = N/2; % to indicate origin at the center of the function
    a = 8; % radius for cylindrical function
    w = 0.4; % decay rate for exponential function
    y = repmat(1:N,N,1);
    x = y’;
    r = sqrt((D-x).^2+(D-y).^2); % definition of radius
    sig = 5; % Variance for gaussian function
    g = (2*pi*sig^2)*exp(-((r.^2))./(2*sig^2));
    figure;mesh(g); title(‘2D Gaussian Function’);
    G = fft2(g);
    G = fftshift(G);
    figure; mesh(abs(G)); title(‘Fourier Transform of Gaussian function’);
    fzf=atan2(imag(G),real(G));
    figure;
    mesh(fzf);
    title(‘Fourier Phase of Gaussian function’);

    Reply

  6. Dear sir, this code is great for generating the gaussian pulse without using matlab toolboxes.
    I want to ask that from the same code, if i want to generate the train of gaussian pulses, how can i generate (without using the matlab function pulstran)? kindly provide the baisc code for it. I have searched over the internet but not found. suppose I want to generate 12 gaussian pulses wtith sum time duration between them. plzz.

    Reply

  7. Hello this is Lingaraj, doing m-tech in belgaum. I have a problem in finding power spectral density using time domain approach,i have a written code found somewhere but i am not getting it. Please leave a reply as soon as possible…

    Reply

  8. Hi Viswanathan,
    This really awesome article on Gaussian wave generation and analysis. I am trying to generate a Gaussian pulse of 20ns and plot the frequency response of it. When I try to simulate it in matlab, it gives me error saying out of memory. I think I am doing something wrong with the code:

    function [ output_args ] = freqencySpectrum( ~ )
    fs=0.00000001; // sampling at twice the highest frequency (20ns =50MHz, so sampling at 100MHz)
    t=-30:1/fs:30;
    sigma = 20E-9;
    A=1;
    p1 = A/(sqrt(2*pi*sigma^2));
    x=p1*(exp(-(t.^2)/(2*sigma^2)));
    figure(1);
    subplot(4,1,1)
    plot(t,x,’b’);
    title([‘Gaussian Pulse sigma=’, num2str(sigma),’s’]);
    xlabel(‘Time(s)’);
    ylabel(‘Amplitude’);

    L=length(x);
    NFFT = 1024;
    X = fftshift(fft(x,NFFT));
    Pxx=X.*conj(X)/(NFFT*NFFT); %computing power with proper scaling
    f = fs*(-NFFT/2:NFFT/2-1)/NFFT; %Frequency Vector

    subplot(4,1,2)
    plot(f,abs(X)/(L),’r’);
    title(‘Magnitude of FFT’);
    xlabel(‘Frequency (Hz)’)
    ylabel(‘Magnitude |X(f)|’);
    xlim([-10 10])

    Pxx=X.*conj(X)/(L*L); %computing power with proper scaling
    subplot(4,1,3)
    plot(f,10*log10(Pxx),’r’);
    title(‘Double Sided – Power Spectral Density’);
    xlabel(‘Frequency (Hz)’)
    ylabel(‘Power Spectral Density- P_{xx} dB/Hz’);

    X = fft(x,NFFT);
    X = X(1:NFFT/2+1);%Throw the samples after NFFT/2 for single sided plot
    Pxx=X.*conj(X)/(L*L);
    f = fs*(0:NFFT/2)/NFFT; %Frequency Vector
    subplot(4,1,4)
    plot(f,10*log10(Pxx),’r’);
    title(‘Single Sided – Power Spectral Density’);
    xlabel(‘Frequency (Hz)’)
    ylabel(‘Power Spectral Density- P_{xx} dB/Hz’);

    end

    Can you help me with this code?
    Thanks,
    KP

    Reply

    • KP,

      Sigma time is not equal to pulse width. Pulse width should be more than the sigma variation… In the following code, the sigma variation is assumed to be 20ns. The pulse width is 10*20ns. I have fixed the code for you.. Here it is

      sigma = 20e-9%sigma variation of the pulse 20ns
      pulseWidth = 10*sigma %double sided pulse width
      fs = 40*1/pulseWidth %Very very high oversampling factor for smooth curve (40 points), FFT with 1024 points is sufficient to cover this
      t=-(pulseWidth/2):1/fs:(pulseWidth/2)

      A=1;
      p1 = A/(sqrt(2*pi*sigma^2));
      x=p1*(exp(-(t.^2)/(2*sigma^2)));
      figure(1);
      subplot(4,1,1)
      plot(t,x,’b’);
      title([‘Gaussian Pulse sigma=’, num2str(sigma),’s’]);
      xlabel(‘Time(s)’);
      ylabel(‘Amplitude’);

      L=length(x);
      NFFT = 1024;
      X = fftshift(fft(x,NFFT));
      Pxx=X.*conj(X)/(NFFT*NFFT); %computing power with proper scaling
      f = fs*(-NFFT/2:NFFT/2-1)/NFFT; %Frequency Vector
      subplot(4,1,2)
      plot(f,abs(X)/(L),’r’);
      title(‘Magnitude of FFT’);
      xlabel(‘Frequency (Hz)’)
      ylabel(‘Magnitude |X(f)|’);

      Pxx=X.*conj(X)/(L*L); %computing power with proper scaling
      subplot(4,1,3)
      plot(f,10*log10(Pxx),’r’);
      title(‘Double Sided – Power Spectral Density’);
      xlabel(‘Frequency (Hz)’)
      ylabel(‘Power Spectral Density- P_{xx} dB/Hz’);

      X = fft(x,NFFT);
      X = X(1:NFFT/2+1);%Throw the samples after NFFT/2 for single sided plot
      Pxx=X.*conj(X)/(L*L);
      f = fs*(0:NFFT/2)/NFFT; %Frequency Vector
      subplot(4,1,4)

      plot(f,10*log10(Pxx),’r’);
      title(‘Single Sided – Power Spectral Density’);
      xlabel(‘Frequency (Hz)’)
      ylabel(‘Power Spectral Density- P_{xx} dB/Hz’);

      Reply

      • Hi Mathuranathan,
        Thank you very much for the reply. I am little confused about something in the code.

        1. Why did you set the pulse width as “pulseWidth = 10*sigma %double sided pulse width”?
        Why 10?

        2. I have attached the output of the simulation. I see that the pulse width is like 60ns, how can I make it 20ns? Should I even lower the sigma value as the pulse width should be greater than the sigma value?

        And in that case, should I change the sampling frequency (fs = 40*1/pulseWidth %) value as well? or 40 will be fine?

        3. Finally, how should I interpret the output of single sided PSD, in the attached figure, I believe that pulse width is 60ns, so the frequency should like 16.66MHz, but in the plot how can I relate it to 16MHz?
        Thanks
        KP

        Reply

        • 1) Subtle difference exists between the terms sigma and pulse width.

          Sigma is the standard deviation of the pulse. The term pulse width depends on where the measurement is taken. Usually, the width of the pulse at half-the maximum value is called Full-Width at Half maximum pulse duration (FWHM).

          What I was intending is not the FWHM width. It is the full duration of the pulse considered for simulation. Note that the Gaussian pulse gradually tapers nears zero on either ends. As an approximation, we must stop at some point. I chose 10 times the value of sigma for this.

          2) To have a full pulse width of 20ns, the sigma has to be lower.

          3) What we plot is the frequency response of a single pulse. It breaks down the pulse in frequency domain and shows the different frequency components that make-up the pulse in time domain.

          The connotation of frequency (that measures periodicity 60MHz => 16.66ns) is applicable only when the pulse is repeated periodically. You might need to repeat the pulse at a desired regular interval and then plot the frequency response to check for any spike in the frequency domain that is indicative of the frequency of the repeated pattern.

          Reply

          • Hello,

            Your example is absolutely fabulous. This is the first time I have seen theory and practical coding to understand.

            1.) I would like to generate a pulse train of Gaussian pulse in time domain with a certain width (let’s say 20 ns in the above example) with a repetition interval of (let’s say 100 ns). In paper, I have to convolve with a Dirac comb. Again let’s say we limit this Dirac comb to a certain number of cycles (let us say 50 cycles).

            How do I do it in Matlab?

            2.) A further addition to this problem, (I cannot get the first one though) is: let us say we make another pulse offset with the pulse from above. Let us say we keep the pulsewidth the same (20 ns), limit the number of cyles to 50 (like above), but change the repetition interval from every other pulse. For the sake of clearly enquiring it, let me draw a diagram here

            *****——————–*****——————–*****——————–**** (P1)
            *****————————-*****—————*****————————-**** (P2)

            Here P1 = pulse one
            P2 pulse 2

            both pulses are in phase at one time (let us say t = t1), then it shifted in time (let us say 4 ns) with the next pulse, then again it is in phase.

            How could I do it matlab?

            Thank you so much for clear explanation. I deeply appreciate it.

            Kind regards
            Chowdhury

            Reply

Post your valuable comments !!!