Power spectral Density Archives

Plot FFT using Python – FFT of sine wave & cosine wave

Key focus: Learn how to plot FFT of sine wave and cosine wave using Python. Understand FFTshift. Plot one-sided, double-sided and normalized spectrum using FFT.

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 Python and how to represent them in frequency domain using FFT. If you are inclined towards Matlab programming, visit here.

This article is part of the book Digital Modulations using Python, ISBN: 978-1712321638 available in ebook (PDF) and Paperback (hardcopy) formats

Sine Wave

In order to generate a sine wave, the first step is to fix the frequency f of the sine wave. For example, we wish to generate a sine wave whose minimum and maximum amplitudes are -1V and +1V respectively. Given the frequency of the sinewave, the next step is to determine the sampling rate.

For baseband signals, the sampling is straight forward. By Nyquist Shannon sampling theorem, for faithful reproduction of a continuous signal in discrete domain, one has to sample the signal at a rate higher than at-least twice the maximum frequency contained in the signal (actually, it is twice the one-sided bandwidth occupied by a real signal. For a baseband signal bandwidth ( to ) and maximum frequency in a given band are equivalent).

For Python implementation, let us write a function to generate a sinusoidal signal using the Python’s Numpy library. Numpy is a fundamental library for scientific computations in Python. In order to use the numpy package, it needs to be imported. Here, we are importing the numpy package and renaming it as a shorter alias np.

import numpy as np

Next, we define a function for generating a sine wave signal with the required parameters.

def sine_wave(f,overSampRate,phase,nCyl):
	"""
	Generate sine wave signal with the following parameters
	Parameters:
		f : frequency of sine wave in Hertz
		overSampRate : oversampling rate (integer)
		phase : desired phase shift in radians
		nCyl : number of cycles of sine wave to generate
	Returns:
		(t,g) : time base (t) and the signal g(t) as tuple
	Example:
		f=10; overSampRate=30;
		phase = 1/3*np.pi;nCyl = 5;
		(t,g) = sine_wave(f,overSampRate,phase,nCyl)
	"""
	fs = overSampRate*f # sampling frequency
	t = np.arange(0,nCyl*1/f-1/fs,1/fs) # time base
	g = np.sin(2*np.pi*f*t+phase) # replace with cos if a cosine wave is desired
	return (t,g) # return time base and signal g(t) as tuple

We note that the function sine wave is defined inside a file named signalgen.py. We will add more such similar functions in the same file. The intent is to hold all the related signal generation functions, in a single file. This approach can be extended to object oriented programming. Now that we have defined the sine wave function in signalgen.py, all we need to do is call it with required parameters and plot the output.

"""
Simulate a sinusoidal signal with given sampling rate
"""
import numpy as np
import matplotlib.pyplot as plt # library for plotting
from signalgen import sine_wave # import the function

f = 10 #frequency = 10 Hz
overSampRate = 30 #oversammpling rate
fs = f*overSampRate #sampling frequency
phase = 1/3*np.pi #phase shift in radians
nCyl = 5 # desired number of cycles of the sine wave

(t,x) = sine_wave(f,overSampRate,phase,nCyl) #function call

plt.plot(t,x) # plot using pyplot library from matplotlib package
plt.title('Sine wave f='+str(f)+' Hz') # plot title
plt.xlabel('Time (s)') # x-axis label
plt.ylabel('Amplitude') # y-axis label
plt.show() # display the figure

Python is an interpreter based software language that processes everything in digital. In order to obtain a smooth sine wave, the sampling rate must be far higher than the prescribed minimum required sampling rate, that is at least twice the frequency – as per Nyquist-Shannon theorem. Hence, we need to sample the input signal at a rate significantly higher than what the Nyquist criterion dictates. Higher oversampling rate requires more memory for signal storage. It is advisable to keep the oversampling factor to an acceptable value.

An oversampling factor of is chosen in the previous function. This is to plot a smooth continuous like sine wave. Thus, the sampling rate becomes . If a phase shift is desired for the sine wave, specify it too.

*Figure 1: A 10Hz sinusoidal wave with 5 cycles and phase shift 1/3π radians*

Sine wave using python — *Figure 1: A 10Hz sinusoidal wave with 5 cycles and phase shift 1/3π radians*

Different representations of FFT:

Since FFT is just a numeric computation of -point DFT, there are many ways to plot the result. The FFT, implemented in Scipy.fftpack package, is an algorithm published in 1965 by J.W.Cooley and
J.W.Tuckey for efficiently calculating the DFT.

The SciPy functions that implement the FFT and IFFT can be invoked as follows

from scipy.fftpack import fft, ifft
X = fft(x,N) #compute X[k]
x = ifft(X,N) #compute x[n]

1. Plotting raw values of DFT:

The x-axis runs from to – representing sample values. Since the DFT values are complex, the magnitude of the DFT is plotted on the y-axis. From this plot we cannot identify the frequency of the sinusoid that was generated.

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft

NFFT=1024 #NFFT-point DFT      
X=fft(x,NFFT) #compute DFT using FFT    

fig1, ax = plt.subplots(nrows=1, ncols=1) #create figure handle
nVals = np.arange(start = 0,stop = NFFT) # raw index for FFT plot
ax.plot(nVals,np.abs(X))      
ax.set_title('Double Sided FFT - without FFTShift')
ax.set_xlabel('Sample points (N-point DFT)')        
ax.set_ylabel('DFT Values')
fig1.show()

*Figure 2: Double sided FFT – without FFTShift*

2. FFT plot – plotting raw values against normalized frequency axis:

In the next version of plot, the frequency axis (x-axis) is normalized to unity. Just divide the sample index on the x-axis by the length of the FFT. This normalizes the x-axis with respect to the sampling rate . Still, we cannot figure out the frequency of the sinusoid from the plot.

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft

NFFT=1024 #NFFT-point DFT  
X=fft(x,NFFT) #compute DFT using FFT     

fig2, ax = plt.subplots(nrows=1, ncols=1) #create figure handle
   
nVals=np.arange(start = 0,stop = NFFT)/NFFT #Normalized DFT Sample points         
ax.plot(nVals,np.abs(X))     
ax.set_title('Double Sided FFT - without FFTShift')        
ax.set_xlabel('Normalized Frequency')
ax.set_ylabel('DFT Values')
fig2.show()

*Figure 3: Double sided FFT with normalized x-axis (0 to 1)*

3. FFT plot – plotting raw values against normalized frequency (positive & negative frequencies):

As you know, in the frequency domain, the values take up both positive and negative frequency axis. In order to plot the DFT values on a frequency axis with both positive and negative values, the DFT value at sample index has to be centered at the middle of the array. This is done by using FFTshift function in Scipy Python. The x-axis runs from to where the end points are the normalized ‘folding frequencies’ with respect to the sampling rate .

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft,fftshift

NFFT=1024 #NFFT-point DFT      
X=fftshift(fft(x,NFFT)) #compute DFT using FFT  

fig3, ax = plt.subplots(nrows=1, ncols=1) #create figure handle
    
fVals=np.arange(start = -NFFT/2,stop = NFFT/2)/NFFT #DFT Sample points        
ax.plot(fVals,np.abs(X))
ax.set_title('Double Sided FFT - with FFTShift')
ax.set_xlabel('Normalized Frequency')
ax.set_ylabel('DFT Values');
ax.autoscale(enable=True, axis='x', tight=True)
ax.set_xticks(np.arange(-0.5, 0.5+0.1,0.1))
fig.show()

*Figure 4: Double sided FFT with normalized x-axis (-0.5 to 0.5)*

4. FFT plot – Absolute frequency on the x-axis vs. magnitude on y-axis:

Here, the normalized frequency axis is just multiplied by the sampling rate. From the plot below we can ascertain that the absolute value of FFT peaks at and . Thus the frequency of the generated sinusoid is . The small side-lobes next to the peak values at and are due to spectral leakage.

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft,fftshift

NFFT=1024     
X=fftshift(fft(x,NFFT))

fig4, ax = plt.subplots(nrows=1, ncols=1) #create figure handle

fVals=np.arange(start = -NFFT/2,stop = NFFT/2)*fs/NFFT
ax.plot(fVals,np.abs(X),'b')
ax.set_title('Double Sided FFT - with FFTShift')
ax.set_xlabel('Frequency (Hz)')         
ax.set_ylabel('|DFT Values|')
ax.set_xlim(-50,50)
ax.set_xticks(np.arange(-50, 50+10,10))
fig4.show()

*Figure 5: Double sided FFT – Absolute frequency on the x-axis vs. magnitude on y-axis*

5. Power Spectrum – Absolute frequency on the x-axis vs. power on y-axis:

Where is the frequency domain representation of the signal . In Python, the power has to be calculated with proper scaling terms.

*Figure 6: Power spectral density using FFT*

Plotting the PSD plot with y-axis on log scale, produces the most encountered type of PSD plot in signal processing.

*Figure 7: Power spectral density (y-axis on log scale) using FFT*

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

Books by the author

Wireless Communication Systems in Matlab Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing

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

Generate color noise using Auto-Regressive (AR) model

Key focus: Learn how to generate color noise using auto regressive (AR) model. Apply Yule Walker equations for generating power law noises: pink noise, Brownian noise.

Auto-Regressive (AR) model

An uncorrelated Gaussian random sequence can be transformed into a correlated Gaussian random sequence using an AR time-series model. If a time series random sequence is assumed to be following an auto-regressive model of form,

where is the uncorrelated Gaussian sequence of zero mean and variance , the natural tendency is to estimate the model parameters . Least Squares method can be applied here to find the model parameters, but the computations become cumbersome as the order increases. Fortunately, the AR model coefficients can be solved for using Yule Walker equations.

Yule Walker equations relate auto-regressive model parameters to the auto-correlation of random process . Finding the model parameters using Yule-Walker equations, is a two step process:

1. Given , estimate auto-correlation of the process . If is already specified as a function, utilize it as it is (see auto-correlation equations for Jakes spectrum or Doppler spectrum in section 11.3.2 in the book).

2. Solve Yule Walker equations to find the model parameters and the noise sigma .

Yule-Walker equations

Yule-Walker equations can be compactly written as

Written in matrix form, the Yule-Walker equations that comprises of a set of linear equations and unknown parameters.

Representing equation (3) in a compact form,

The AR model parameters can be found by solving

After solving for the model parameters , the noise variance can be found by applying the estimated values of in equation (2) by setting . The aryule command (in Matlab and Python’s spectrum package) efficiently solves the Yule-Walker equations using Levinson Algorithm [1][2]. Once the model parameters are obtained, the AR model can be implemented as an \emph{infinte impulse response (IIR)} filter of form

Example: power law noise generation

The power law in the power spectrum characterizes the fluctuating observables in many natural systems. Many natural systems exhibit some noise which is a stochastic process with a power spectral density having a power exponent that can take values . Simply put, noise is a colored noise with a power spectral density of over its entire frequency range.

The noise can be classified into different types based on the value of .

Violet noise – = -2, the power spectral density is proportional to .
Blue noise – = -1, the power spectral density is proportional to .
White noise – = 0, the power spectral density is flat across the whole spectrum.
Pink noise – = 1, the power spectral density is proportional to , i.e, it decreases by per octave with increase in frequency.
Brownian noise – = 2, the power spectral density is proportional to , therefore it decreases by per octave with increase in frequency.

The power law noise can be generated by sequencing a zero-mean white noise through an auto-regressive (AR) filter of order :

where, is a zero-mean white noise process. Referring the AR generation method described in [3], the coefficients of the AR filter can be generated as

which can be implemented as an infinite impulse response (IIR) filter using the filter transfer function described in equation (6).

The following script implements this method and the sample results are plotted in the next Figure.

Refer the book for the Matlab code

Figure 1: Simulated color noise samples and their PSD estimates: pink noise (α =1) and Brown noise (α =2)

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

References

[1] Gene H. Golub, Charles F. Van Loan, Matrix Computations, ISBN-9780801854149, Johns Hopkins University Press, 1996, p. 143.↗
[2] J. Durbin, The fitting of time series in models, Review of the International Statistical Institute, 28:233-243, 1960.↗
[3] Kasdin, N.J. Discrete Simulation of Colored Noise and Stochastic Processes and Power Law Noise Generation, Proceedings of the IEEE, Vol. 83, No. 5, 1995, pp. 802-827.↗

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

Books by the author

Wireless Communication Systems in Matlab Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing

Topics in this chapter

Random Variables - Simulating Probabilistic Systems
● Introduction
● Plotting the estimated PDF
● Univariate random variables
□ Uniform random variable
□ Bernoulli random variable
□ Binomial random variable
□ Exponential random variable
□ Poisson process
□ Gaussian random variable
□ Chi-squared random variable
□ Non-central Chi-Squared random variable
□ Chi distributed random variable
□ Rayleigh random variable
□ Ricean random variable
□ Nakagami-m distributed random variable
● Central limit theorem - a demonstration
● Generating correlated random variables
□ Generating two sequences of correlated random variables
□ Generating multiple sequences of correlated random variables using Cholesky decomposition
● Generating correlated Gaussian sequences
□ Spectral factorization method
□ Auto-Regressive (AR) model

Generating colored noise with Jakes PSD: Spectral factorization

The aim of this article is to demonstrate the application of spectral factorization method in generating colored noise having Jakes power spectral density. Before continuing, I urge the reader to go through this post: Introduction to generating correlated Gaussian sequences.

This article is part of the book
● Wireless Communication Systems in Matlab (second edition), ISBN: 979-8648350779 available in ebook (PDF) format and Paperback (hardcopy) format.

In spectral factorization method, a filter is designed using the desired frequency domain characteristics (like PSD) to transform an uncorrelated Gaussian sequence into a correlated sequence . In the model shown in Figure 1, the input to the LTI system is a white noise whose amplitude follows Gaussian distribution with zero mean and variance and the power spectral density (PSD) of the white noise is a constant across all frequencies.

The white noise sequence drives the LTI system with frequency response producing the signal of interest . The PSD of the output process is therefore

Figure 1: Relationship among various power spectral densities in a filtering process

If the desired power spectral density of the colored noise sequence is given, assuming , the impulse response of the LTI filter can be found by taking the inverse Fourier transform of the frequency response

Once, the impulse response of the filter is obtained, the colored noise sequence can be produced by driving the filter with a zero-mean white noise sequence of unit variance.

Example: Generating colored noise with Jakes PSD

For example, we wish to generate a Gaussian noise sequence whose power spectral density follows the normalized Jakes power spectral density (see section 11.3.2 in the book) given by

Applying spectral factorization method, the frequency response of the desired filter is

The impulse response of the filter is [1]

where, is the fractional Bessel function of the first kind, is the sampling interval for implementing the digital filter and is a constant. The impulse response of the filter can be normalized by dividing by .

The filter can be implemented as a finite impulse response (FIR) filter structure. However, the FIR implementation requires that the impulse response be truncated to a reasonable length. Such truncation leads to ringing effects due to Gibbs phenomenon. To avoid distortions due to truncation, the filter impulse response is usually windowed using a window function such as Hamming window.

where, the Hamming window is defined as

The function given in the book in section 2.6.1 implements a windowed Jakes filter using the aforementioned equations. The impulse response and the spectral characteristics of the filter are plotted in Figure 2.

Figure 2: Impulse response & spectrum of windowed Jakes filter ( fmax = 10Hz; Ts = 0:01s; N = 512)

A white noise can be transformed into colored noise sequence with Jakes PSD, by processing the white noise through the implemented filter. The script (given in the book in section 2.6.1) illustrates this concept by transforming a white noise sequence into a colored noise sequence. The simulated noise samples and its PSD are plotted in Figure 3.

Figure 3: Simulating colored noise with Jakes PSD

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

Reference

[1] Jeruchim et., al, Simulation of communication systems – modeling, methodology, and techniques, second edition, Kluwer academic publishers, 2002, ISBN: 0306462672.↗

Books by the author

Wireless Communication Systems in Matlab Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing

Topics in this chapter

Generate correlated Gaussian sequence (colored noise)

Key focus: Colored noise sequence (a.k.a correlated Gaussian sequence), is a non-white random sequence, with non-constant power spectral density across frequencies.

Introduction

Speaking of Gaussian random sequences such as Gaussian noise, we generally think that the power spectral density (PSD) of such Gaussian sequences is flat.We should understand that the PSD of a Gausssian sequence need not be flat. This bring out the difference between white and colored random sequences, as captured in Figure 1.

A white noise sequence is defined as any random sequence whose PSD is constant across all frequencies. Gaussian white noise is a Gaussian random sequence, whose amplitude is gaussian distributed and its PSD is a constant. Viewed in another way, a constant PSD in frequency domain implies that the average auto-correlation function in time-domain is an impulse function (Dirac-delta function). That is, the amplitude of noise at any given time instant is correlated only with itself. Therefore, such sequences are also referred as uncorrelated random sequences. White Gaussian noise processes are completely characterized by its mean and variance.

This article is part of the book
● Wireless Communication Systems in Matlab (second edition), ISBN: 979-8648350779 available in ebook (PDF) format and Paperback (hardcopy) format.

Figure 1: Power spectral densities of white noise and colored noise

A colored noise sequence is simply a non-white random sequence, whose PSD varies with frequency. For a colored noise, the amplitude of noise at any given time instant is correlated with the amplitude of noise occurring at other instants of time. Hence, colored noise sequences will have an auto-correlation function other than the impulse function. Such sequences are also referred as correlated random sequences. Colored
Gaussian noise processes are completely characterized by its mean and the shaped of power spectral density (or the shape of auto-correlation function).

In mobile channel model simulations, it is often required to generate correlated Gaussian random sequences with specified mean and power spectral density (like Jakes PSD or Gaussian PSD given in section 11.3.2 in the book). An uncorrelated Gaussian sequence can be transformed into a correlated sequence through filtering or linear transformation, that preserves the Gaussian distribution property of amplitudes, but alters only the correlation property (equivalently the power spectral density). We shall see two methods to generate colored Gaussian noise for given mean and PSD shape

• Spectral factorization method
• Auto-regressive (AR) model

Motivation

Let’s say we observe a real world signal that has an arbitrary spectrum . We would like to describe the long sequence of using very few parameters, as in applications like linear predictive coding (LPC). The modeling approach, described here, tries to answer the following two questions:

• Is it possible to model the first order (mean/variance) and second order (correlations, spectrum) statistics of the signal just by shaping a white noise spectrum using a transfer function ? (see Figure 1).
• Does this produce the same statistics (spectrum, correlations, mean and variance) for a white noise input ?

If the answer is yes to the above two questions, we can simply set the modeled parameters of the system and excite the system with white noise, to produce the desired real world signal. This reduces the amount to data we wish to transmit in a communication system application. This approach can be used to transform an uncorrelated white Gaussian noise sequence to a colored Gaussian noise sequence with desired spectral properties.

Linear time invariant (LTI) system model

In the given model, the random signal is observed. Given the observed signal , the goal here is to find a model that best describes the spectral properties of under the following assumptions
• The sequence is WSS (wide sense stationary) and ergodic.
• The input sequence to the LTI system is white noise, whose amplitudes follow Gaussian distribution with zero-mean and variance with flat the power spectral density.
• The LTI system is BIBO (bounded input bounded output) stable.

Read the continuation of this post : Spectral factorization method

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

Reference

[1] Jeruchim et., al, Simulation of communication systems – modeling, methodology, and techniques, second edition, Kluwer academic publishers, 2002, ISBN: 0306462672.↗

Books by the author

Wireless Communication Systems in Matlab Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing

Topics in this chapter

Chirp Signal – FFT & PSD in Matlab & Python

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

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

Introduction

All the signals discussed so far do not change in frequency over time. Obtaining a signal with time-varying frequency is of main focus here. A signal that varies in frequency over time is called “chirp”. The frequency of the chirp signal can vary from low to high frequency (up-chirp) or from high to low frequency (low-chirp).

Observation

Chirp signals/signatures are encountered in many applications ranging from radar, sonar, spread spectrum, optical communication, image processing, doppler effect, motion of a pendulum, as gravitation waves, manifestation as Frequency Modulation (FM), echo location ^[1] etc.

Mathematical Description:

A linear chirp signal sweeps the frequency from low to high frequency (or vice-versa) linearly. One approach to generate a chirp signal is to concatenate a series of segments of sine waves each with increasing(or decreasing) frequency in order. This method introduces discontinuities in the chirp signal due to the mismatch in the phases of each such segments. Modifying the equation of a sinusoid to generate a chirp signal is a better approach.

The equation for generating a sinusoidal (cosine here) signal with amplitude A, angular frequency and initial phase is

This can be written as a function of instantaneous phase

where is the instantaneous phase of the sinusoid and it is linear in time. The time derivative of instantaneous phase is equal to the angular frequency of the sinusoid – which in case is a constant in the above equation.

Instead of having the phase linear in time, let’s change the phase to quadratic form and thus non-linear.

for some constant .

Therefore, the equation for chirp signal takes the following form,

The first derivative of the phase, which is the instantaneous angular frequency becomes a function of time, which is given by

The time-varying frequency in Hertz is given by

In the above equation, the frequency is no longer a constant, rather it is of time-varying nature with initial frequency given by . Thus, from the above equation, given a time duration T, the rate of change of frequency is given by

where, is the starting frequency of the sweep, is the frequency at the end of the duration T.

Substituting (7) & (8) in (6)

From (6) and (8)

where is a constant which will act as the initial phase of the sweep.

Thus the modified equation for generating a chirp signal (from equations (5) and (10)) is given by

where the time-varying frequency function is given by

Generation of Chirp signal, computing its Fourier Transform using FFT and power spectral density (PSD) in Matlab is shown as example, for Python code, please refer the book Digital Modulations using Python.

Generating a chirp signal without using in-built “chirp” Function in Matlab:

Implement a function that describes the chirp using equation (11) and (12). The starting frequency of the sweep is and the frequency at time is . The initial phase forms the final part of the argument in the following function

function x=mychirp(t,f0,t1,f1,phase)
%Y = mychirp(t,f0,t1,f1) generates samples of a linear swept-frequency
%   signal at the time instances defined in timebase array t.  The instantaneous
%   frequency at time 0 is f0 Hertz.  The instantaneous frequency f1
%   is achieved at time t1.
%   The argument 'phase' is optional. It defines the initial phase of the
%   signal degined in radians. By default phase=0 radian
    
if nargin==4
    phase=0;
end
    t0=t(1);
    T=t1-t0;
    k=(f1-f0)/T;
    x=cos(2*pi*(k/2*t+f0).*t+phase);
end

The following wrapper script utilizes the above function and generates a chirp with starting frequency at the start of the time base and at which is the end of the time base. From the PSD plot, it can be ascertained that the signal energy is concentrated only upto 25 Hz

fs=500; %sampling frequency
t=0:1/fs:1; %time base - upto 1 second

f0=1;% starting frequency of the chirp
f1=fs/20; %frequency of the chirp at t1=1 second
x = mychirp(t,f0,1,f1); 
subplot(2,2,1)
plot(t,x,'k');
title(['Chirp Signal']);
xlabel('Time(s)');
ylabel('Amplitude');

FFT and power spectral density

As with other signals, describes in the previous posts, let’s plot the FFT of the generated chirp signal and its power spectral density (PSD).

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,2,2)
plot(f,abs(X)/(L),'r');
title('Magnitude of FFT');
xlabel('Frequency (Hz)')
ylabel('Magnitude |X(f)|');
xlim([-50 50])


Pxx=X.*conj(X)/(NFFT*NFFT); %computing power with proper scaling
subplot(2,2,3)
plot(f,10*log10(Pxx),'r');
title('Double Sided - Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('Power Spectral Density- P_{xx} dB/Hz');
xlim([-100 100])

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

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

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

References:

[1] Patrick Flandrin,“Chirps everywhere”,CNRS — Ecole Normale Supérieure de Lyon

Topics in this chapter

Books by the author

Wireless Communication Systems in Matlab Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing

Gaussian Pulse – FFT & PSD in Matlab & Python

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.

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

The Fourier Transform of a Gaussian pulse preserves its shape.

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

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

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

The following code generates a Gaussian Pulse with ( ). 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)

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

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

Figure 3: One-sided power spectral density of Gaussian pulse

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

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

Topics in this chapter

Books by the author

Wireless Communication Systems in Matlab Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing

Generating Basic signals – Rectangular Pulse and Power Spectral Density using FFT

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, square wave, 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 book Digital Modulations using Matlab : Build Simulation Models from Scratch, ISBN: 978-1521493885 available in ebook (PDF) format (click here) and Paperback (hardcopy) format (click here)
Wireless Communication Systems in Matlab, ISBN: 978-1720114352 available in ebook (PDF) format (click here) and Paperback (hardcopy) format (click here).

Rectangular pulse: mathematical description

An isolated rectangular pulse of amplitude A and duration T is represented mathematically as

where

The Fourier transform of isolated rectangular pulse g(t) is

where, the sinc function is given by

Thus, the Fourier Transform pairs are

The Fourier Transform describes the spectral content of the signal at various frequencies. For a given signal g(t), the Fourier Transform is given by

where, the absolute value gives the magnitude of the frequency components (amplitude spectrum) and are their corresponding phase (phase spectrum) . For the rectangular pulse, the amplitude spectrum is given as

The amplitude spectrum peaks at f=0 with value equal to AT. The nulls of the spectrum occur at integral multiples of 1/T, i.e, ( )

Generating an isolated rectangular pulse in Matlab:

An isolated rectangular pulse of unit amplitude and width w (the factor T in equations above ) can be generated easily with the help of in-built function – rectpuls(t,w) command in Matlab. As an example, a unit amplitude rectangular pulse of duration is generated.

fs=500; %sampling frequency
T=0.2; %width of the rectangule pulse in seconds

t=-0.5:1/fs:0.5; %time base

x=rectpuls(t,T); %generating the square wave

plot(t,x,'k');
title(['Rectangular Pulse width=', num2str(T),'s']);
xlabel('Time(s)');
ylabel('Amplitude');

Amplitude spectrum using FFT:

Matlab’s FFT function is utilized for computing the Discrete Fourier Transform (DFT). The magnitude of FFT is plotted. From the following plot, it can be noted that the amplitude of the peak occurs at f=0 with peak value . The nulls in the spectrum are located at ().

L=length(x);
NFFT = 1024;
X = fftshift(fft(x,NFFT)); %FFT with FFTshift for both negative & positive frequencies
f = fs*(-NFFT/2:NFFT/2-1)/NFFT; %Frequency Vector

figure;
plot(f,abs(X)/(L),'r');
title('Magnitude of FFT');
xlabel('Frequency (Hz)')
ylabel('Magnitude |X(f)|');

Power spectral density (PSD) using FFT:

The distribution of power among various frequency components is plotted next. The first plot shows the double-side Power Spectral Density which includes both positive and negative frequency axis. The second plot describes the PSD only for positive frequency axis (as the response is just the mirror image of negative frequency axis).

figure;
Pxx=X.*conj(X)/(L*L); %computing power with proper scaling
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
plot(f,10*log10(Pxx),'r');
title('Single Sided - Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('Power Spectral Density- P_{xx} dB/Hz');

Magnitude and phase spectrum:

The phase spectrum of the rectangular pulse manifests as series of pulse trains bounded between 0 and , provided the rectangular pulse is symmetrically centered around sample zero. This is explained in the reference here and the demo below.

clearvars;
x = [ones(1,7) zeros(1,127-13) ones(1,6)];
subplot(3,1,1); plot(x,'k');
title('Rectangular Pulse'); xlabel('Sample#'); ylabel('Amplitude');

NFFT = 127;
X = fftshift(fft(x,NFFT)); %FFT with FFTshift for both negative & positive frequencies
f = (-NFFT/2:NFFT/2-1)/NFFT; %Frequency Vector

subplot(3,1,2); plot(f,abs(X),'r');
title('Magnitude Spectrum'); xlabel('Frequency (Hz)'); ylabel('|X(f)|');

subplot(3,1,3); plot(f,atan2(imag(X),real(X)),'r');
title('Phase Spectrum'); xlabel('Frequency (Hz)'); ylabel('\angle X(f)');

*Magnitude and phase spectrum of rectangular pulse*

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

Topics in this chapter

Books by the author

Wireless Communication Systems in Matlab Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing

Generating Basic signals – Square Wave and Power Spectral Density using FFT

Note: There is a rating embedded within this post, please visit this post to rate it.

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.

Significance of Square Waves

The most logical way of transmitting information across a communication channel is through a stream of square pulse – a distinct pulse for ‘0‘ and another for ‘1‘. Digital signals are graphically represented as square waves with certain symbol/bit period. Square waves are also used universally in switching circuits, as clock signals synchronizing various blocks of digital circuits, as reference clock for a given system domain and so on.

Square wave manifests itself as a wide range of harmonics in frequency domain and therefore can cause electromagnetic interference. Square waves are periodic and contain odd harmonics when expanded as Fourier Series (where as signals like saw-tooth and other real word signals contain harmonics at all integer frequencies). Since a square wave literally expands to infinite number of odd harmonic terms in frequency domain, approximation of square wave is another area of interest. The number of terms of its Fourier Series expansion, taken for approximating the square wave is often seen as Gibbs Phenomenon, which manifests as ringing effect at the corners of the square wave in time domain (visual explanation here).

True Square waves are a special class of rectangular waves with 50% duty cycle. Varying the duty cycle of a rectangular wave leads to pulse width modulation, where the information is conveyed by changing the duty-cycle of each transmitted rectangular wave.

How to generate a square wave in Matlab

If you know the trick of generating a sine wave in Matlab, the task is pretty much simple. Square wave is generated using “square” function in Matlab. The command sytax – square(t,dutyCycle) – generates a square wave with period for the given time base. The command behaves similar to “sin” command (used for generating sine waves), but in this case it generates a square wave instead of a sine wave. The argument – dutyCycle is optional and it defines the desired duty cycle of the square wave. By default (when the dutyCycle argument is not supplied) the square wave is generated with (50%) duty cycle.

f=10; %frequency of sine wave in Hz
overSampRate=30; %oversampling rate
fs=overSampRate*f; %sampling frequency
duty_cycle=50; % Square wave with 50% Duty cycle (default)
nCyl = 5; %to generate five cycles of sine wave

t=0:1/fs:nCyl*1/f; %time base

x=square(2*pi*f*t,duty_cycle); %generating the square wave

plot(t,x,'k');
title(['Square Wave f=', num2str(f), 'Hz']);
xlabel('Time(s)');
ylabel('Amplitude');

Power Spectral Density using FFT

Let’s check out how the generated square wave will look in frequency domain. The Fast Fourier Transform (FFT) is utilized here. As discussed in the article here, there are numerous ways to plot the response of FFT. Single Sided power spectral density is plotted first, followed by the Double-sided power spectral density.

Single Sided Power Spectral Density

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

plot(f,10*log10(Pxx),'r');
title('Single Sided Power Spectral Density');
xlabel('Frequency (Hz)')
ylabel('Power Spectral Density- P_{xx} dB/Hz');
ylim([-45 -5])

Double Sided Power Spectral Density

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

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

Note: There is a rating embedded within this post, please visit this post to rate it.

Topics in this chapter

Books by the author

Wireless Communication Systems in Matlab Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing

Plot FFT using Matlab – FFT of sine wave & cosine wave

Key focus: Learn how to plot FFT of sine wave and cosine wave using Matlab. Understand FFTshift. Plot one-sided, double-sided and normalized spectrum.

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. If you are inclined towards Python programming, visit here.

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

Sine Wave

In order to generate a sine wave in Matlab, the first step is to fix the frequency of the sine wave. For example, I intend to generate a f=10 Hz sine wave whose minimum and maximum amplitudes are and respectively. Now that you have determined the frequency of the sinewave, the next step is to determine the sampling rate. Matlab is a software that processes everything in digital. In order to generate/plot a smooth sine wave, the sampling rate must be far higher than the prescribed minimum required sampling rate which is at least twice the frequency – as per Nyquist Shannon Theorem. A oversampling factor of is chosen here – this is to plot a smooth continuous-like sine wave (If this is not the requirement, reduce the oversampling factor to desired level). Thus the sampling rate becomes . If a phase shift is desired for the sine wave, specify it too.

f=10; %frequency of sine wave
overSampRate=30; %oversampling rate
fs=overSampRate*f; %sampling frequency
phase = 1/3*pi; %desired phase shift in radians
nCyl = 5; %to generate five cycles of sine wave

t=0:1/fs:nCyl*1/f; %time base

x=sin(2*pi*f*t+phase); %replace with cos if a cosine wave is desired
plot(t,x);
title(['Sine Wave f=', num2str(f), 'Hz']);
xlabel('Time(s)');
ylabel('Amplitude');

Representing in Frequency Domain

Representing the given signal in frequency domain is done via Fast Fourier Transform (FFT) which implements Discrete Fourier Transform (DFT) in an efficient manner. Usually, power spectrum is desired for analysis in frequency domain. In a power spectrum, power of each frequency component of the given signal is plotted against their respective frequency. The command computes the -point DFT. The number of points – – in the DFT computation is taken as power of (2) for facilitating efficient computation with FFT. A value of is chosen here. It can also be chosen as next power of 2 of the length of the signal.

Different representations of FFT:

Since FFT is just a numeric computation of -point DFT, there are many ways to plot the result.

1. Plotting raw values of DFT:

NFFT=1024; %NFFT-point DFT      
X=fft(x,NFFT); %compute DFT using FFT        
nVals=0:NFFT-1; %DFT Sample points       
plot(nVals,abs(X));      
title('Double Sided FFT - without FFTShift');        
xlabel('Sample points (N-point DFT)')        
ylabel('DFT Values');

2. FFT plot – plotting raw values against Normalized Frequency axis:

NFFT=1024; %NFFT-point DFT      
X=fft(x,NFFT); %compute DFT using FFT        
nVals=(0:NFFT-1)/NFFT; %Normalized DFT Sample points         
plot(nVals,abs(X));      
title('Double Sided FFT - without FFTShift');        
xlabel('Normalized Frequency')       
ylabel('DFT Values');

3. FFT plot – plotting raw values against normalized frequency (positive & negative frequencies):

As you know, in the frequency domain, the values take up both positive and negative frequency axis. In order to plot the DFT values on a frequency axis with both positive and negative values, the DFT value at sample index has to be centered at the middle of the array. This is done by using function in Matlab. The x-axis runs from to where the end points are the normalized ‘folding frequencies’ with respect to the sampling rate .

NFFT=1024; %NFFT-point DFT      
X=fftshift(fft(x,NFFT)); %compute DFT using FFT      
fVals=(-NFFT/2:NFFT/2-1)/NFFT; %DFT Sample points        
plot(fVals,abs(X));      
title('Double Sided FFT - with FFTShift');       
xlabel('Normalized Frequency')       
ylabel('DFT Values');

4. FFT plot – Absolute frequency on the x-axis Vs Magnitude on Y-axis:

NFFT=1024;      
X=fftshift(fft(x,NFFT));         
fVals=fs*(-NFFT/2:NFFT/2-1)/NFFT;        
plot(fVals,abs(X),'b');      
title('Double Sided FFT - with FFTShift');       
xlabel('Frequency (Hz)')         
ylabel('|DFT Values|');

5. Power Spectrum – Absolute frequency on the x-axis Vs Power on Y-axis:

The following is the most important representation of FFT. It plots the power of each frequency component on the y-axis and the frequency on the x-axis. The power can be plotted in linear scale or in log scale. The power of each frequency component is calculated as

Where is the frequency domain representation of the signal . In Matlab, the power has to be calculated with proper scaling terms (since the length of the signal and transform length of FFT may differ from case to case).

NFFT=1024;
L=length(x);         
X=fftshift(fft(x,NFFT));         
Px=X.*conj(X)/(NFFT*L); %Power of each freq components       
fVals=fs*(-NFFT/2:NFFT/2-1)/NFFT;        
plot(fVals,Px,'b');      
title('Power Spectral Density');         
xlabel('Frequency (Hz)')         
ylabel('Power');

If you wish to verify the total power of the signal from time domain and frequency domain plots, follow this link.
Plotting the power spectral density (PSD) plot with y-axis on log scale, produces the most encountered type of PSD plot in signal processing.

NFFT=1024;      
L=length(x);         
X=fftshift(fft(x,NFFT));         
Px=X.*conj(X)/(NFFT*L); %Power of each freq components       
fVals=fs*(-NFFT/2:NFFT/2-1)/NFFT;        
plot(fVals,10*log10(Px),'b');        
title('Power Spectral Density');         
xlabel('Frequency (Hz)')         
ylabel('Power');

6. Power Spectrum – One-Sided frequencies

In this type of plot, the negative frequency part of x-axis is omitted. Only the FFT values corresponding to to sample points of -point DFT are plotted. Correspondingly, the normalized frequency axis runs between to . The absolute frequency (x-axis) runs from to .

L=length(x);        
NFFT=1024;       
X=fft(x,NFFT);       
Px=X.*conj(X)/(NFFT*L); %Power of each freq components       
fVals=fs*(0:NFFT/2-1)/NFFT;      
plot(fVals,Px(1:NFFT/2),'b','LineSmoothing','on','LineWidth',1);         
title('One Sided Power Spectral Density');       
xlabel('Frequency (Hz)')         
ylabel('PSD');

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

For further reading

[1] Power spectral density – MIT opencourse ware↗

Topics in this chapter

Books by the author

Wireless Communication Systems in Matlab Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing

Compute signal power in Matlab

Calculating the energy and power of a signal was discussed in one of the previous posts. Here, we will verify the calculation of signal power using Discrete Fourier Transform (DFT) in Matlab. Check here to know more on the concept of power and energy.

The total power of a signal can be computed using the following equation

For other forms of equations: refer here.

Case Study:

Let x(t) be a sine wave of amplitude A and frequency f_c represented by the following equation.

When represented in frequency domain, it will look like the one on the right side plot in the following figure. This is evident from the fact that the sinewave can be mathematically represented by applying Euler’s formula.

Taking the Fourier transform of x(t) to represent it in frequency domain,

When considering the amplitude part,the above decomposition gives two spikes of amplitude A/2 on either side of the frequency domain at f_c and -f_c.

*Figure 1: Time domain representation and frequency spectrum of a pure sinewave*

Squaring the amplitudes gives the magnitude of power of the individual spikes/frequency components. The power spectrum is plotted below.

*Figure 2: Power spectrum of a pure sinewave*

Thus if the pure sinewave is of amplitude A=1V and frequency=100Hz, the power spectrum will have two spikes of value at 100 Hz and -100 Hz frequencies. The total power will be

Let’s verify this through simulation.

Simulation and Verification

A sine wave of 100 Hz frequency and amplitude 1V is taken for the experiment.

A=1; %Amplitude of sine wave
Fc=100; %Frequency of sine wave
Fs=1000; %Sampling rate - oversampled by the rate of 10
Ts=1/Fs; %Sampling period
nCycles=200; %Number of cycles of the sinewave

subplot(2,1,1);
t=0:Ts:nCycles/Fc-Ts; %Time base
x=A*sin(2*pi*Fc*t); %Sinusoidal function
stem(t,x); %Plot command

A sinusoidal wave of 10 cycles is plotted here

*Figure 3: Sine wave generated in Matlab*

Matlab’s Norm function:

Matlab’s basic installation comes with “norm” function. The p-norm in Matlab is computed as

By default, the single argument norm function computed 2-norm given as

To compute the total power of the signal x[n] (as in equation (1) above), all we have to do is – compute norm(x), square it and divide by the length of the signal.

L=length(x);
P=(norm(x)^2)/L;
sprintf('Power of the Signal from Time domain %f',P);

The above given piece of code will result in the following output

Power of the Signal from Time domain 0.500000

Verifying the total Power by DFT : Frequency Domain

Here, the total power is verified by applying DFT on the sinusoidal sequence. The sinusoidal sequence x[n] is represented in frequency domain X[f] using Matlab’s FFT function. The power associated with each frequency point is computed as

Finally, the total power is calculated as the sum of all the points in the frequency domain representation.

X = fft(x);
Px=sum(X.*conj(X))/(L^2); %Compute power with proper scaling.
subplot(2,1,2)
% Plot single-sided amplitude spectrum.
stem(Px);
sprintf('Total Power of the Signal from DFT %f',P);

Result:

Total Power of the Signal from DFT 0.500000

*Figure 4: Power spectrum of a pure sinewave simulated in Matlab*

A word on Matlab’s FFT: Matlab’s FFT is optimized for faster performance if the transform length is a power of 2. The following snippet of code simply calls “fft” without the transform length. In this case, the window length and the transform length are the same. This is to simplify the calculation of power. You can re-write the call to the FFT routine with transform length set to next power of two which is greater than or equal to the window length (sequence length). Then the step to compute total power will be differing slightly in the denominator. But that will not improve resolution (Remember : zero padding to compute FFT will not improve resolution).

Also note that in the above simulation we are using a pure sinusoid. The entire sequence of sinusoid defined all the cycles completely. There are no discontinuities in the sequence. If you call FFT with the transform length set to next power of 2 (as given in Matlab manuals), it will pad additional zeros to the sequence and creates discontinuities in the FFT computation. This will lead to spectral leakage. FFT and spectral leakage is discussed here.

Rate this article: Note: There is a rating embedded within this post, please visit this post to rate it.

References:

[1] Sanjay Lall,”Norm and Vector spaces”,Information Systems Laboratory,Stanford University.↗

Books by the author

Wireless Communication Systems in Matlab Second Edition(PDF) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Python (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart	Digital Modulations using Matlab (PDF ebook) Note: There is a rating embedded within this post, please visit this post to rate it. Checkout Added to cart
Hand-picked Best books on Communication Engineering Best books on Signal Processing

Cookie	Duration	Description
cookielawinfo-checbox-analytics	11 months	This cookie is set by GDPR Cookie Consent plugin. The cookie is used to store the user consent for the cookies in the category "Analytics".
cookielawinfo-checbox-analytics	11 months	This cookie is set by GDPR Cookie Consent plugin. The cookie is used to store the user consent for the cookies in the category "Analytics".
cookielawinfo-checbox-functional	11 months	The cookie is set by GDPR cookie consent to record the user consent for the cookies in the category "Functional".
cookielawinfo-checbox-functional	11 months	The cookie is set by GDPR cookie consent to record the user consent for the cookies in the category "Functional".
cookielawinfo-checbox-others	11 months	This cookie is set by GDPR Cookie Consent plugin. The cookie is used to store the user consent for the cookies in the category "Other.
cookielawinfo-checbox-others	11 months	This cookie is set by GDPR Cookie Consent plugin. The cookie is used to store the user consent for the cookies in the category "Other.
cookielawinfo-checkbox-necessary	11 months	This cookie is set by GDPR Cookie Consent plugin. The cookies is used to store the user consent for the cookies in the category "Necessary".
cookielawinfo-checkbox-performance	11 months	This cookie is set by GDPR Cookie Consent plugin. The cookie is used to store the user consent for the cookies in the category "Performance".
viewed_cookie_policy	11 months	The cookie is set by the GDPR Cookie Consent plugin and is used to store whether or not user has consented to the use of cookies. It does not store any personal data.

Introduction

Sine Wave

Different representations of FFT:

1. Plotting raw values of DFT:

2. FFT plot – plotting raw values against normalized frequency axis:

3. FFT plot – plotting raw values against normalized frequency (positive & negative frequencies):

4. FFT plot – Absolute frequency on the x-axis vs. magnitude on y-axis:

5. Power Spectrum – Absolute frequency on the x-axis vs. power on y-axis:

Books by the author

Topics in this chapter

Auto-Regressive (AR) model

Yule-Walker equations

Example: <img decoding="async" src="https://s0.wp.com/latex.php?latex=1%2Ff%5E%5Calpha&amp;bg=ffffff&amp;fg=000&amp;s=3&amp;c=20201002" alt="1/f^\alpha" width="54" height="29"> power law noise generation

References

Books by the author

Topics in this chapter

Example: Generating colored noise with Jakes PSD

Reference

Books by the author

Topics in this chapter

Introduction

Motivation

Linear time invariant (LTI) system model

Reference

Books by the author

Topics in this chapter

Introduction

Observation

Mathematical Description:

Generating a chirp signal without using in-built “chirp” Function in Matlab:

FFT and power spectral density

References:

Topics in this chapter

Books by the author

Gaussian Pulse : Mathematical description:

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

Double Sided and Single Power Spectral Density using FFT:

Topics in this chapter

Books by the author

Rectangular pulse: mathematical description

Generating an isolated rectangular pulse in Matlab:

Amplitude spectrum using FFT:

Power spectral density (PSD) using FFT:

Magnitude and phase spectrum:

Topics in this chapter

Books by the author

Significance of Square Waves

How to generate a square wave in Matlab

Power Spectral Density using FFT

Single Sided Power Spectral Density

Double Sided Power Spectral Density

Topics in this chapter

Books by the author

Introduction

Sine Wave

Representing in Frequency Domain

Different representations of FFT:

1. Plotting raw values of DFT:

2. FFT plot – plotting raw values against Normalized Frequency axis:

3. FFT plot – plotting raw values against normalized frequency (positive & negative frequencies):

4. FFT plot – Absolute frequency on the x-axis Vs Magnitude on Y-axis:

5. Power Spectrum – Absolute frequency on the x-axis Vs Power on Y-axis:

6. Power Spectrum – One-Sided frequencies

For further reading

Topics in this chapter

Books by the author

Case Study:

Simulation and Verification

Matlab’s Norm function:

Verifying the total Power by DFT : Frequency Domain

Result:

References:

Books by the author

Example: power law noise generation