Generating Basic Signals – Gaussian Pulse and Power Spectral Density using FFT

1 Star2 Stars3 Stars4 Stars5 Stars (3 votes, average: 4.67 out of 5)
Loading...

Articles in this series:
How to Interpret FFT results – complex DFT, frequency bins and FFTShift
How to Interpret FFT results – obtaining Magnitude and Phase information
FFT and Spectral Leakage
How to plot FFT using Matlab – FFT of basic signals : Sine and Cosine waves
Generating Basic signals – Square Wave and Power Spectral Density using FFT
Generating Basic signals – Rectangular Pulse and Power Spectral Density using FFT
Generating Basic Signals – Gaussian Pulse and Power Spectral Density using FFT (this article)
Chirp Signal – Frequency Sweeping – FFT and power spectral density
Constructing the Auto Correlation Matrix using FFT

 For more such examples check this ebook : Simulation of Digital Communications using Matlab – by Mathuranathan Viswanathan

Introduction

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.

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) = \frac{1}{\sqrt{2 \pi } \sigma} e^{- \frac{t^2}{2 \sigma^2}} $$

The Fourier Transform of a Gaussian pulse preserves its shape.

$$ \begin{align*}
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{align*} $$

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 \). Then,
$$ \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.
$$ \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 and its Fourier Transform using FFT:

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.

Gaussian Pulse Matlab FFT

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

Gaussian Pulse Double Sided Power Spectral Density

Gaussian Pulse Single Sided Power Spectral Density

Rate this article: 1 Star2 Stars3 Stars4 Stars5 Stars (3 votes, average: 4.67 out of 5)

Loading...

Articles in this series:
How to Interpret FFT results – complex DFT, frequency bins and FFTShift
How to Interpret FFT results – obtaining Magnitude and Phase information
FFT and Spectral Leakage
How to plot FFT using Matlab – FFT of basic signals : Sine and Cosine waves
Generating Basic signals – Square Wave and Power Spectral Density using FFT
Generating Basic signals – Rectangular Pulse and Power Spectral Density using FFT
Generating Basic Signals – Gaussian Pulse and Power Spectral Density using FFT (this article)
Chirp Signal – Frequency Sweeping – FFT and power spectral density
Constructing the Auto Correlation Matrix using FFT

 For more such examples check this ebook : Simulation of Digital Communications using Matlab – by Mathuranathan Viswanathan

Recommended Signal Processing Books

  • pankaj

    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.

  • Lingaraj kadanavar

    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…

  • KP

    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

    • 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’);

      • KP

        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

        • 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.

          • KP

            Thank you very much Mathuranathan, this will help.

          • Golam Kibria Chowdhury

            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