Chirp Signal – Frequency Sweeping – FFT and power spectral density

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

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, cosineGaussian pulsesquarewaveisolated 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.

Checkout this ebook : Simulation of Digital Communications using Matlab – by Mathuranathan Viswanathan

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 \(\omega_0\) and initial phase \(\phi\) is
$$x(t)=A cos(\omega_0 t + \phi ) \;\;\;\;\;\;\;\;\;\;\;\;\; (1) $$

This can be written as a function of instantaneous phase
$$x(t)=A cos(\theta (t)) \;\;\;\;\;\;\;\;\;\;\;\;\; (2)$$

where \(\theta (t) = \omega_0(t)+ \phi \) is the instantaneous phase of the sinusoid and it is linear in time. The time derivative of instantaneous phase \(\theta(t)\) is equal to the angular frequency \(\omega\) of the sinusoid – which in case is a constant in the above equation.
$$ \omega (t) = \frac{d}{dt} \theta (t) \;\;\;\;\;\;\;\;\;\;\;\;\; (3)$$

Instead of having the phase linear in time, let’s change the phase to quadratic form and thus non-linear.
$$\theta(t) = 2 \pi \alpha t^2 + 2 \pi f_0 t + \phi \;\;\;\;\;\;\;\;\;\;\;\;\; (4)$$

for some constant \(\alpha\). The first derivative the phase, which is the instantaneous angular frequency becomes
$$ \omega_i(t)=\frac{d}{dt} \theta (t) = 4 \pi \alpha t + 2 \pi f_0 \;\;\;\;\;\;\;\;\;\;\;\;\; (5)$$

The time-varying frequency in \(Hz\) is given by
$$ f_i (t) = 2 \alpha t + f_0 \;\;\;\;\;\;\;\;\;\;\;\;\; (6)$$

In the above equation, the frequency is no longer a constant, rather it is of time-varying nature with initial frequency given by \(f_0\). Thus, from the above equation, given a time duration \(T\), the rate of change of frequency is given by
$$ k = 2 \alpha = \frac{f_1-f_0}{T} \;\;\;\;\;\;\;\;\;\;\;\;\; (7)$$
where, \(f_0\) is the starting frequency of the sweep, f_1 is the frequency at the end of the duration \(T\).

The time varying frequency function for a phase derivative equal to \( \omega_i(t)= 2 \pi (\alpha t + f_0 ) \) is
$$f_i (t) = \alpha t + f_0  \;\;\;\;\;\;\;\;\;\;\;\;\; (8)$$

Substituting \((6)\) & \((7)\) in \((5)\)
$$ \omega_i(t)=\frac{d}{dt} \theta (t) = 2 \pi (kt+f_0) \;\;\;\;\;\;\;\;\;\;\;\;\; (9) $$

From \((3)\) and \((9)\)
$$ \begin{align*}
\theta (t) &= \int \omega_i(t) \, dt\\
&= 2 \pi \int (kt+f_0) \, dt = 2 \pi (k\frac{t^2}{2}+f_0 t) + \phi_0\\
&= 2 \pi (k\frac{t^2}{2}+f_0 t) + \phi_0 = 2 \pi (\frac{k}{2} t +f_0 ) t + \phi_0
\end{align*} \;\;\;\;\;\;\;\;\;\;\;\;\; (10)$$

where \( \phi_0 \) is a constant which will act as the initial phase of the sweep.

Thus the modified equation for generating a chirp signal is given by
$$ x(t) = A cos (2 \pi f(t) t + \phi_0) \;\;\;\;\;\;\;\;\;\;\;\;\; (11)$$
where the time-varying frequency function is given by
$$ f(t)= \frac{k}{2} t +f_0 \;\;\;\;\;\;\;\;\;\;\;\;\; (12)$$

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 \(f_0\) and the frequency at time \(t_1\) is \(f_1\). The initial phase \(\phi_0\) forms the final part of the argument in the following function

The following wrapper script utilizes the above function and generates a chirp with starting frequency \(f_0=1Hz\) at the start of the time base and \(f_1=25 Hz \) at \(t_1=1s\) 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\)

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

Chirp signal FFT and power spectral density in Matlab

References:

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

Rate this article: 1 Star2 Stars3 Stars4 Stars5 Stars (13 votes, average: 4.54 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
Chirp Signal – Frequency Sweeping – FFT and power spectral density (this article)
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

  • raphael

    Hey nice article. I have a question: why is the amplitude of the second plot (magnitude of fft) not one? Because all the frequencys in your chirp signal have the amplitude of one? Am I missing something

    • The frequency of the chirp signal varies with time, hence the amplitude of each frequency component does not stay the same. To have a same amplitude for all frequencies, the signal needs to have 1 complete cycle for each frequency component. For chirp, the frequency continuously changes from one time instant to the next, you cannot pin-point a cycle.

      White noise is an excellent example that has a flat magnitude for all frequencies. See here: http://www.gaussianwaves.com/2013/11/simulation-and-analysis-of-white-noise-in-matlab/

  • Lloyd Blythen

    Hi Mathuranathan. Thanks for the detailed explanation. Is there a typo? Formula (2) correctly states “x(t) = Acos(Ɵ(t))”. But immediately after that you state “where Ɵ(t) = ω0 + Ф is the instantaneous phase and it is linear in time.” Looks like there’s a ‘t’ missing: I think it should read “where Ɵ(t) = ω0t + Ф is the instantaneous phase and it is linear in time.”

    • Dear Lloyd,
      Thanks for spotting the typo. It is corrected now.

  • Clark Kent

    Does anyone know any way to produce a chirp in Matlab that follows a custom polynomial rather than the pre-built in functions provided in the chirp() function?

    Does the chirp() function allow for more customization than I am aware of?

    I want to produce a chirp based on the following polynomial (determined using polyfit() on a set of data points)…

    y2 = p(1).*x.^10 + p(2).*x.^9 + p(3).*x.^8 + p(4).*x.^7 + p(5).*x.^6 + p(6).*x.^5 + p(7).*x.^4 + p(8).*x.^3 + p(9).*x.^2 + p(10).*x + p(11) ;