dft Archives - GaussianWaves

Window function – figure of merits

Key focus: Window function smooths the observed signal over the edges. Analysis of some important parameters to help select the window for an application.

Spectral leakage

As we know, the DFT operation can be viewed as processing a signal through a set of filter banks with bandwidth Δf centered on the bin (frequency) of interest (Figure 1). To accurately resolve the frequency component in each bin, we desire the resolution bandwidth Δf to be as small as possible. In effect, each filter filter bank is viewed as filter having a brick-wall response with narrow bandwidth. If the signal under study contains several spectral components and broadband noise, the computed amplitude of the DFT output in each bin (filter bank) is significantly affected by the noise contained in the corresponding bandwidth of the filter bank.

*Figure 1: DFT (Discrete Fourier Transform) viewed as processing through a set of filter banks*

*Figure 1: DFT (Discrete Fourier Transform) viewed as processing through a set of filter banks*

In reality, signals are of time-limited nature and nothing can be known about the signal beyond the measured interval. When the time-limited signal slice is subjected to analysis using FFT algorithm (used for computing DFT), the FFT implicitly assumes that the signal essentially repeats itself after the observed interval. This may lead to discontinuities at the edges of each slice, that causes the energy contained in each frequency bin to spill into other bins. This phenomenon is called spectral leakage.

Equivalent Noise Bandwidth:^[2]

To suppress the spectral leakage, the signals are multiplied with a window function so as to smooth the discontinuity at the edges of the FFT slices. As a result, the choice of window function affects the amount of signal and noise that goes inside each filter bank. Hence the amount of noise that gets accumulated into each filter bank is affected by the bandwidth of the chosen window function. In effect, the amplitude estimate in each frequency bin is influenced by the amount of accumulated noise contributed by the window function of choice. In other words, the noise floor displayed at the FFT output, varies according to the window of choice.

Equivalent noise bandwidth (ENBW), often specified in terms of FFT bins, is the bandwidth of a fictitious brick-wall filter such that the amount of noise power accumulated inside the brick-wall filter is same as the noise power accumulated when processing white noise through the selected window.

*Figure 2: Illustrating equivalent noise bandwidth (ENBW) (figure not to scale*)

In discrete domain, the equivalent noise bandwidth B_enbw can be computed using time domain samples as

\[B_{enbw} = L \frac{\sum_{n}|w[n]|^2}{ \left| \sum_{n} w[n] \right|^2}\]

where L is the length of the window.

ENBW is a commonly used metric to characterize the variation in the displayed noise floor if a non rectangular window is used. If noise floor measurements are desired, the ENBW correction factor need to be applied to account for the variation in the noise floor [3]. Table 1 below, lists the ENBW (bins) and ENBW correction factor for some well known window functions.

Following figure illustrates the equivalent noise bandwidth illustrated for Hann window.

*Figure 3: Equivalent noise bandwidth (ENBW) illustrated for Hann window*

Coherent power gain:

In time-domain, windowing a signal reduces the amplitude of the signal at its boundaries (Figure 4). This causes a reduction in the amplitude of the spectral component when FFT used to visualize the signal in frequency domain. This reduction in amplitude in the spectral components is characterized as coherent power gain (which actually is a loss). So, when FFT is computed on a windowed signal, in order to compensate for the loss in the amplitude, we simply add the coherent power gain to the FFT output.

*Figure 4: Windowing a 10 Hz sinusoidal signal with Hann window of length L = 151 (sampling frequency F_s=160 Hz)*

Coherent power gain of a window w[n] of length L is given by

\[\text{Coherent power gain (dB)} = 20 \; log_{10} \left( \frac{\sum_n w[n]}{L} \right)\]

Following plot depicts the coherent power gain (i.e, the reduction in the FFT magnitude when the input signal is processed with a window) of a windowed sinusoidal signal of frequency 10 Hz. The length of the window is L = 151 points and the simulation assumes an oversampling factor of 16 (i.e, F_s=160 Hz). The FFT length is N_FFT =2048.

Table 1 below, lists the coherent power gain for some well known window functions.

*Figure 5: Coherent power gain – illustrated for various window functions*

Scalloping loss

As discussed above, the FFT (equivalently the DFT) operation can be seen as processing a signal through a set of filter banks with bandwidth Δf centered on the successive bin frequencies. The frequency resolution (Δf) of FFT depends on the size of FFT (N_FFT) and the sampling frequency F_s.

\[\Delta f = \frac{F_s}{N_{FFT}}\]

Therefore, the FFT process can measure amplitude of tones that are multiples of the frequency resolution Δf. That is, if the frequency of the observed signal matches a particular bin frequency, the signal amplitude will be maximum at that bin.

However, if the tone of the signal falls between two bins, the signal power is spread between the adjacent bins and hence the displayed signal amplitude/power is reduced.

Scalloping loss quantifies the worst case reduction in the signal level, when the tone of the observed signal falls exactly mid way between two bin frequencies. The worst case scalloping loss is calculated as

\[\text{Scalloping loss} = \frac{\sum_n w[n] e^{\left( -j \frac{\pi}{N_{FFT}}n \right)}}{\sum_n w[n]}\]

In reality, scalloping loss depends on the type of window and the relationship between the signal tones to the bin frequencies. Table 1 below, lists the scalloping loss for some well known window functions.

*Figure 6: Scalloping loss illustrated for processing a signal with Boxcar and Bartlett window*

Figure of merits for window functions

*Table 1: Figure of merits for various window functions*

List of window functions

Boxcar (Rectangular)

\[w[n] = \begin{cases} 1, \quad 0 \leq n \leq L-1 \\ 0, \quad \text{otherwise} \end{cases}\]

Figure 1: Boxcarr (rectangular) window *(L=51)*

Barthann

\[ w[n] = 0.62 – 0.48 \left| \frac{n}{L – 1} – 0.5 \right| + 0.38 cos \left(2 \pi \left| \frac{n}{L – 1} – 0.5 \right| \right), \quad 0 \leq n \leq L-1\]

Bartlett

\[w[n] = \frac{2}{L-1} \left(\frac{L-1}{2} – \left|n – \frac{L-1}{2}\right| \right),\quad 0 \leq n \leq L-1\]

Blackman

\[w[n] = 0.42 – 0.5 \cos(2\pi n/L) + 0.08 \cos(4\pi n/L), \quad 0 \leq n \leq L-1\]

Blackman-Harris

\[w[n]=a_0 – a_1 \cos \left ( \frac{2 \pi n}{L} \right)+ a_2 \cos \left ( \frac{4 \pi n}{L} \right)- a_3 \cos \left ( \frac{6 \pi n}{L} \right), \quad 0 \leq n \leq L-1\] \[a_0=0.35875;\quad a_1=0.48829;\quad a_2=0.14128;\quad a_3=0.01168\]

Bohman

\[w[n] = \left( 1 – |x|\right) cos \left( \pi |x| \right) + \frac{1}{\pi} sin \left( \pi |x| \right) , \quad -1 \leq x \leq 1 \]

Cosine

\[w[n] = sin \left( \frac{\pi}{ L } \left(n + 0.5 \right) \right),\quad 0 \leq n \leq L-1\]

Flattop

\[w[n]=a_0 – a_1 \cos \left ( \frac{2 \pi n}{L} \right)+ a_2 \cos \left ( \frac{4 \pi n}{L} \right)- a_3 \cos \left ( \frac{6 \pi n}{L} \right), \quad 0 \leq n \leq L-1 \\ a_0=0.21557895;\; a_1=0.41663158 \\ a_2=0.277263158;\; a_3=0.006947368 \]

Hamming

\[w[n] = 0.54 – 0.46 \cos\left(\frac{2\pi{n}}{L-1}\right), \quad 0 \leq n \leq L-1\]

Hann

\[w[n] = 0.5 – 0.5 \cos\left(\frac{2\pi{n}}{L-1}\right), \quad 0 \leq n \leq L-1\]

Nuttall

\[ w[n]=a_0 – a_1 \cos \left ( \frac{2 \pi n}{L} \right)+ a_2 \cos \left ( \frac{4 \pi n}{L} \right)- a_3 \cos \left ( \frac{6 \pi n}{L} \right), \quad 0 \leq n \leq L-1 \\ a_0=0.3635819;\; a_1=0.4891775;\\ a_2=0.1365995;\; a_3=0.0106411\]

Parzen

\[w[n] = \begin{cases} 1- 6 \left( \frac{|n|}{L/2} \right)^2 + 6 \left( \frac{|n|}{L/2} \right)^3, & 0 \leq |n| \leq (L-1)/4 \\ 2 \left( 1 – \frac{|n|}{L/2} \right)^3, & (L-1)/4 < |n| \leq (L-1)/2 \end{cases}\]

Triangular

\[w[n] = 1 – \frac{|n|}{L/2}, \quad n = -L/2, \cdots, -1, 0, 1, \cdots, L/2\]

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

References:

[1] Harris, “On the Use of Windows for Harmonic Analysis with the Discrete Fourier Transform”, proceedings of IEEE, vol. 66, no. 1, January 1978 .↗
[2] Stefan Scholl, “Exact Signal Measurements using FFT Analysis”,Microelectronic Systems Design Research Group, TU Kaiserslautern, Germany.↗

Additional Resources:

[1] If you are particularly interested in calculating the Equivalent Noise Bandwidth of Analog Filters refer – Thomas J. Karras,”Equivalent Noise Bandwidth Analysis from Transfer Functions”,NASA Technical Note,NASA TN D-2842.↗
[2] Strategies for Choosing Windows : Michael Cerna and Audrey F. Harvey,”The Fundamentals of FFT-Based Signal Analysis and Measurement”,National Instruments Application Notes 041.↗

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

Equivalent noise bandwidth (ENBW) of window functions

Key focus: Equivalent noise bandwidth (ENBW), is the bandwidth of a fictitious brick-wall filter that allows same amount of noise as a window function. Learn how to calculate ENBW in applications involving window functions and FFT operation.

FFT and spectral leakage

Here is an article with illustrations on – FFT slices and the phenomenon of spectral leakage,

Hence, to suppress the spectral leakage, the signals are multiplied with a window function so as to smooth the discontinuity at the edges of the FFT slices. As a result, the choice of window function affects the amount of signal and noise that goes inside each filter bank. Hence the amount of noise that gets accumulated into each filter bank is affected by the bandwidth of the chosen window function. In effect, the amplitude estimate in each frequency bin is influenced by the amount of accumulated noise contributed by the window function of choice. In other words, the noise floor displayed at the FFT output, varies according to the window of choice.

Equivalent noise bandwidth

In other words, in frequency domain, the power accumulated in the chosen window is given by

\[P_{window} = \int_{-f}^{f} |W(f)|^2 df\]

Then we wish to find the equivalent noise bandwidth B_enbw, such that

\[P_{brick-wall} = B_{enbw} |W(f_0)|^2 = P_{window}\]

where f₀ is the frequency at which the power peaks and P_brick-wall is the power accumulated in the fictitious rectangular window of bandwidth B_enbw.

Therefore, the equivalent noise bandwidth B_enbw is given by

\[B_{enbw} = \frac{\int_{-f}^{f} |W(f)|^2 df}{|W(f_0)|^2}\]

Translating to discrete domain, the equivalent noise bandwidth can be computed using DFT samples as

\[B_{enbw} =\frac{\sum_{k=0}^{N-1}|W[k]|^2}{|W[k_0]|^2}\]

where, k₀ is the index at which the magnitude of FFT output is maximum and N is the window length. Applying Parseval’s theorem and with a window of length L, B_enbw can also be computed using time domain samples as

\[B_{enbw} = L \frac{\sum_{n}|w[n]|^2}{ \left| \sum_{n} w[n] \right|^2}\]

Following figure illustrates the equivalent noise bandwidth illustrated for Hann window.

Python script

ENBW can be calculated from the time domain samples in a straightforward manner. Following Python 3.0 script calculates the ENBW for some well-known window functions provided as part of Scipy module↗.

#Author : Mathuranathan Viswanathan for gaussianwaves.com
#Date: 13 Sept 2020
#Script to calculate equivalent noise bandwidth (ENBW) for some well known window functions

import numpy as np
import pandas as pd
from scipy import signal

def equivalent_noise_bandwidth(window):
    #Returns the Equivalent Noise BandWidth (ENBW)
    return len(window) * np.sum(window**2) / np.sum(window)**2

def get_enbw_windows():
    #Return ENBW for all the following windows as a dataframe
    window_names = ['boxcar','barthann','bartlett','blackman','blackmanharris','bohman','cosine','exponential','flattop','hamming','hann','nuttall','parzen','triang']
    
    df = pd.DataFrame(columns=['Window','ENBW (bins)','ENBW correction (dB)'])
    for window_name in window_names:
        method_name = window_name
        func_to_run = getattr(signal, method_name) #map window names to window functions in scipy package
        L = 16384 #Number of points in the output window
        window = func_to_run(L) #call the functions
        
        enbw = equivalent_noise_bandwidth(window) #compute ENBW
        
        df = df.append({'Window': window_name.title(),'ENBW (bins)':round(enbw,3),'ENBW correction (dB)': round(10*np.log10(enbw),3)},ignore_index=True)
                       
    return df

df = get_enbw_windows() #call the function
df.head(14) #display dataframe

The resulting output is displayed in the following table (dataframe)

*Table 1: Table of equivalent noise bandwidth for some well known window functions*

As an example, the following plot depicts the difference in the noise floor of FFT output of a noise added sine wave that is processed through Boxcar and Flattop window.

*Figure 4: Noise floor and ENBW for flattop & boxcar window (FFT output) for noise added 10 Hz sinewave (oversampling factor = 16, Fs = 160 Hz, window length L =2048, SNR = 30 dB)*

Continue reading on Window function and figure of merits…

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

Reference

[1] Stefan Scholl, “Exact Signal Measurements using FFT Analysis”,Microelectronic Systems Design Research Group, TU Kaiserslautern, Germany.↗

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

Parseval’s theorem – derivation

The Parseval’s theorem (a.k.a Plancherel theorem) expresses the energy of a signal in time-domain in terms of the average energy in its frequency components.

Suppose if the x[n] is a discrete-time sequence of complex numbers of length N : x_n={x₀,x₁,…,x_N-1}, its N-point discrete Fourier transform (DFT)[1] : X_k={X₀,X₁,…,X_N-1} is given by

\[X[k] = \sum_{n=0}^{N-1} x[n] e^{-j\frac{2 \pi}{N} k n} \]

The inverse discrete Fourier transform is given by

\[\tilde{x}[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j \frac{2 \pi}{N} kn}\]

Suppose if x[n] and y[n] are two such sequences that follows the above definitions, the Parseval’s theorem is written as

\[\boxed{ \sum_{n=0}^{N-1} x[n] y^{\ast}[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] Y^{\ast}[k]} \]

where, * indicates conjugate operation.

Deriving Parseval’s theorem

\[\begin{aligned} \sum_{n=0}^{N-1} x[n] y^{\ast}[n] &= \sum_{n=0}^{N-1} x[n] \left(\frac{1}{N} \sum_{k=0}^{N-1} Y[k] e^{j\frac{2 \pi}{N} k n} \right )^\ast \\ &= \frac{1}{N}\sum_{n=0}^{N-1} x[n] \sum_{k=0}^{N-1} Y^\ast[k] e^{-j\frac{2 \pi}{N} k n} \\ &= \frac{1}{N} \sum_{k=0}^{N-1} Y^\ast[k] \cdot \sum_{n=0}^{N-1} x[n] e^{-j\frac{2 \pi}{N} k n} \\ &= \frac{1}{N} \sum_{k=0}^{N-1} X[k] Y^\ast[k] \end{aligned}\]

If x[n] = y[n], then

\[ \begin{aligned} \sum_{n=0}^{N-1} x[n] x^{\ast}[n] &= \frac{1}{N} \sum_{k=0}^{N-1} X[k] X^\ast[k] \\ \Rightarrow \sum_{n=0}^{N-1} |x[n]|^2 &= \frac{1}{N} \sum_{k=0}^{N-1} |X[k]|^2\end{aligned} \]

Time-domain and frequency domain representations are equivalent manifestations of the same signal. Therefore, the energy of the signal computed from time domain samples must be equal to the total energy computed from frequency domain representation.

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

Reference

[1] Bertrand Delgutte and Julie Greenberg , “The discrete Fourier Transform”, Biomedical Signal and Image Processing Spring 2005, MIT web

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

Interpret FFT results – obtaining magnitude and phase information

In the previous post, Interpretation of frequency bins, frequency axis arrangement (fftshift/ifftshift) for complex DFT were discussed. In this post, I intend to show you how to interpret FFT results and obtain magnitude and phase information.

Outline

For the discussion here, lets take an arbitrary cosine function of the form $x(t)= A cos \left(2 \pi f_c t + \phi \right)$ and proceed step by step as follows

● Represent the signal $x(t)$ in computer (discrete-time) and plot the signal (time domain)
● Represent the signal in frequency domain using FFT ($ X[k]$)
● Extract amplitude and phase information from the FFT result
● Reconstruct the time domain signal from the frequency domain samples

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

Discrete-time domain representation

Consider a cosine signal of amplitude $A=0.5$, frequency $f_c=10 Hz$ and phase $phi= \pi/6$ radians (or $30^{\circ}$ )

\[x(t) = 0.5 cos \left( 2 \pi 10 t + \pi/6 \right)\]

In order to represent the continuous time signal $x(t)$ in computer memory, we need to sample the signal at sufficiently high rate (according to Nyquist sampling theorem). I have chosen a oversampling factor of $32$ so that the sampling frequency will be $f_s = 32 \times f_c $, and that gives $640$ samples in a $2$ seconds duration of the waveform record.

A = 0.5; %amplitude of the cosine wave
fc=10;%frequency of the cosine wave
phase=30; %desired phase shift of the cosine in degrees
fs=32*fc;%sampling frequency with oversampling factor 32
t=0:1/fs:2-1/fs;%2 seconds duration

phi = phase*pi/180; %convert phase shift in degrees in radians
x=A*cos(2*pi*fc*t+phi);%time domain signal with phase shift

figure; plot(t,x); %plot the signal

Represent the signal in frequency domain using FFT

Lets represent the signal in frequency domain using the FFT function. The FFT function computes $N$-point complex DFT. The length of the transformation $N$ should cover the signal of interest otherwise we will some loose valuable information in the conversion process to frequency domain. However, we can choose a reasonable length if we know about the nature of the signal.

For example, the cosine signal of our interest is periodic in nature and is of length $640$ samples (for 2 seconds duration signal). We can simply use a lower number $N=256$ for computing the FFT. In this case, only the first $256$ time domain samples will be considered for taking FFT. No need to worry about loss of information in this case, as the $256$ samples will have sufficient number of cycles using which we can calculate the frequency information.

N=256; %FFT size
X = 1/N*fftshift(fft(x,N));%N-point complex DFT

In the code above, $fftshift$ is used only for obtaining a nice double-sided frequency spectrum that delineates negative frequencies and positive frequencies in order. This transformation is not necessary. A scaling factor $1/N$ was used to account for the difference between the FFT implementation in Matlab and the text definition of complex DFT.

3a. Extract amplitude of frequency components (amplitude spectrum)

The FFT function computes the complex DFT and the hence the results in a sequence of complex numbers of form $X_{re} + j X_{im}$. The amplitude spectrum is obtained

\[|X[k]| = \sqrt{X_{re}^2 + X_{im}^2 } \]

For obtaining a double-sided plot, the ordered frequency axis (result of fftshift) is computed based on the sampling frequency and the amplitude spectrum is plotted.

df=fs/N; %frequency resolution
sampleIndex = -N/2:N/2-1; %ordered index for FFT plot
f=sampleIndex*df; %x-axis index converted to ordered frequencies
stem(f,abs(X)); %magnitudes vs frequencies
xlabel('f (Hz)'); ylabel('|X(k)|');

3b. Extract phase of frequency components (phase spectrum)

Extracting the correct phase spectrum is a tricky business. I will show you why it is so. The phase of the spectral components are computed as

\[\angle X[k] = tan^{-1} \left( \frac{X_{im}}{X_{re}} \right)\]

That equation looks naive, but one should be careful when computing the inverse tangents using computers. The obvious choice for implementation seems to be the $atan$ function in Matlab. However, usage of $atan$ function will prove disastrous unless additional precautions are taken. The $atan$ function computes the inverse tangent over two quadrants only, i.e, it will return values only in the $[-\pi/2 , \pi/2]$ interval. Therefore, the phase need to be unwrapped properly. We can simply fix this issue by computing the inverse tangent over all the four quadrants using the $atan2(X_{img},X_{re})$ function.

Lets compute and plot the phase information using $atan2$ function and see how the phase spectrum looks

phase=atan2(imag(X),real(X))*180/pi; %phase information
plot(f,phase); %phase vs frequencies

The phase spectrum is completely noisy. Unexpected !!!. The phase spectrum is noisy due to fact that the inverse tangents are computed from the $ratio$ of imaginary part to real part of the FFT result. Even a small floating rounding off error will amplify the result and manifest incorrectly as useful phase information (read how a computer program approximates very small numbers).

To understand, print the first few samples from the FFT result and observe that they are not absolute zeros (they are very small numbers in the order $10^{-16}$. Computing inverse tangent will result in incorrect results.

>> X(1:5)
ans =
   1.0e-16 *
  -0.7286            -0.3637 - 0.2501i  -0.4809 - 0.1579i  -0.3602 - 0.5579i   0.0261 - 0.4950i
>> atan2(imag(X(1:5)),real(X(1:5)))
ans =
    3.1416   -2.5391   -2.8244   -2.1441   -1.5181

The solution is to define a tolerance threshold and ignore all the computed phase values that are below the threshold.

X2=X;%store the FFT results in another array
%detect noise (very small numbers (eps)) and ignore them
threshold = max(abs(X))/10000; %tolerance threshold
X2(abs(X)<threshold) = 0; %maskout values that are below the threshold
phase=atan2(imag(X2),real(X2))*180/pi; %phase information
plot(f,phase); %phase vs frequencies

The recomputed phase spectrum is plotted below. The phase spectrum has correctly registered the $30^{\circ}$ phase shift at the frequency $f=10 Hz$. The phase spectrum is anti-symmetric ($\phi=-30^{\circ}$ at $f=-10 Hz$ ), which is expected for real-valued signals.

Reconstruct the time domain signal from the frequency domain samples

Reconstruction of the time domain signal from the frequency domain sample is pretty straightforward

x_recon = N*ifft(ifftshift(X),N); %reconstructed signal
t = [0:1:length(x_recon)-1]/fs; %recompute time index 
plot(t,x_recon);%reconstructed signal

The reconstructed signal has preserved the same initial phase shift and the frequency of the original signal. Note: The length of the reconstructed signal is only $256$ sample long ($\approx 0.8$ seconds duration), this is because the size of FFT is considered as $N=256$. Since the signal is periodic it is not a concern. For more complicated signals, appropriate FFT length (better to use a value that is larger than the length of the signal) need to be used.

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

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

Interpret FFT, complex DFT, frequency bins & FFTShift

Key focus: Interpret FFT results, complex DFT, frequency bins, fftshift and ifftshift. Know how to use them in analysis using Matlab and 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

Four types of Fourier Transforms:

Often, one is confronted with the problem of converting a time domain signal to frequency domain and vice-versa. Fourier Transform is an excellent tool to achieve this conversion and is ubiquitously used in many applications. In signal processing, a time domain signal can be continuous or discrete and it can be aperiodic or periodic. This gives rise to four types of Fourier transforms.

From Table 1, we note note that when the signal is discrete in one domain, it will be periodic in other domain. Similarly, if the signal is continuous in one domain, it will be aperiodic (non-periodic) in another domain. For simplicity, let’s not venture into the specific equations for each of the transforms above. We will limit our discussion to DFT, that is widely available as part of software packages like Matlab, Scipy(python) etc.., however we can approximate other transforms using DFT.

Real and Complex versions of transforms

For each of the listed transforms above, there exist a real version and complex version. The real version of the transform, takes in a real numbers and gives two sets of real frequency domain points – one set representing coefficients over $cosine$ basis function and the other set representing the co-efficient over $sine$ basis function. The complex version of the transforms represent positive and negative frequencies in a single array. The complex versions are flexible that it can process both complex valued signals and real valued signals. The following figure captures the difference between real DFT and complex DFT.

Real DFT:

Consider the case of N-point $real$ DFT , it takes in N samples of (real-valued) time domain waveform $x[n]$ and gives two arrays of length $N/2+1$ each set projected on cosine and sine functions respectively.

\[\begin{align} X_{re}[k] &= \phantom{-}\frac{2}{N} \displaystyle{ \sum_{n=0}^{N-1} x[n] \cdot cos\left( \frac{2 \pi k n}{N} \right)} \\ X_{im}[k] &= -\frac{2}{N} \sum_{n=0}^{N-1} \displaystyle{ x[n] \cdot sin\left( \frac{2 \pi k n}{N} \right)} \end{align}\]

Here, the time domain index $n$ runs from $0 \rightarrow N$, the frequency domain index $k$ runs from $0 \rightarrow N/2$

The real-valued time domain signal $x [n]$ can be synthesized from the real DFT pairs as

\[x[n] = \displaystyle{ \sum_{k=0}^{N/2} X_{re}[K] \cdot cos\left( \frac{2 \pi k n}{N} \right) – X_{im}[K] \cdot sin\left( \frac{2 \pi k n}{N} \right)}\]

Caveat: When using the synthesis equation, the values $X_{re}[0]$ and $X_{re}[N/2] $ must be divided by two. This problem is due to the fact that we restrict the analysis to real-values only. These type of problems can be avoided by using complex version of DFT.

Complex DFT:

Consider the case of N-point complex DFT, it takes in N samples of complex-valued time domain waveform $x[n]$ and produces an array $X[k]$ of length N.

\[ X[k]= \displaystyle{\sum_{n=0}^{N-1} x[n] e^{-j2 \pi k n/N}}\]

The arrays values are interpreted as follows

● $X[0]$ represents DC frequency component
● Next $N/2$ terms are positive frequency components with $X[N/2]$ being the Nyquist frequency (which is equal to half of sampling frequency)
● Next $N/2-1$ terms are negative frequency components (note: negative frequency components are the phasors rotating in opposite direction, they can be optionally omitted depending on the application)

The corresponding synthesis equation (reconstruct $x[n]$ from frequency domain samples $X[k]$) is

\[x[n]= \displaystyle{ \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j2 \pi k n/N}} \]

From these equations we can see that the real DFT is computed by projecting the signal on cosine and sine basis functions. However, the complex DFT projects the input signal on exponential basis functions (Euler’s formula connects these two concepts).

When the input signal in the time domain is real valued, the complex DFT zero-fills the imaginary part during computation (That’s its flexibility and avoids the caveat needed for real DFT). The following figure shows how to interpret the raw FFT results in Matlab that computes complex DFT. The specifics will be discussed next with an example.

Figure 2: Interpretation of frequencies in complex DFT output

Fast Fourier Transform (FFT)

The FFT function in Matlab is an algorithm published in 1965 by J.W.Cooley and J.W.Tuckey for efficiently calculating the DFT. It exploits the special structure of DFT when the signal length is a power of 2, when this happens, the computation complexity is significantly reduced. FFT length is generally considered as power of 2 – this is called $radix-2$ FFT which exploits the twiddle factors. The FFT length can be odd as used in this particular FFT implementation – Prime-factor FFT algorithm where the FFT length factors into two co-primes.

FFT is widely available in software packages like Matlab, Scipy etc.., FFT in Matlab/Scipy implements the complex version of DFT. Matlab’s FFT implementation computes the complex DFT that is very similar to above equations except for the scaling factor. For comparison, the Matlab’s FFT implementation computes the complex DFT and its inverse as

\[ \begin{align} X[k] &= \displaystyle{ \phantom {\frac{1}{N}}\sum_{n=0}^{N-1} x[n] e^{-j2 \pi k n/N}} \\ x[n] &= \displaystyle{ \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j2 \pi k n/N}} \end{align}\]

The Matlab commands that implement the above equations are FFT and IFFT) respectively. The corresponding syntax is as follows

X = fft(x,N) %compute X[k]
x = ifft(X,N) %compute x[n]

Interpreting the FFT results

Lets assume that the $x[n]$ is the time domain cosine signal of frequency $f_c=10Hz$ that is sampled at a frequency $f_s=32*fc$ for representing it in the computer memory.

fc=10;%frequency of the carrier
fs=32*fc;%sampling frequency with oversampling factor=32
t=0:1/fs:2-1/fs;%2 seconds duration
x=cos(2*pi*fc*t);%time domain signal (real number)
subplot(3,1,1); plot(t,x);hold on; %plot the signal
title('x[n]=cos(2 \pi 10 t)'); xlabel('n'); ylabel('x[n]');

Figure 3: A 2 seconds record of 10 Hz cosine wave

Lets consider taking a $N=256$ point FFT, which is the $8^{th}$ power of $2$.

N=256; %FFT size
X = fft(x,N);%N-point complex DFT
%output contains DC at index 1, Nyquist frequency at N/2+1 th index
%positive frequencies from index 2 to N/2
%negative frequencies from index N/2+1 to N

Note: The FFT length should be sufficient to cover the entire length of the input signal. If $N$ is less than the length of the input signal, the input signal will be truncated when computing the FFT. In our case, the cosine wave is of 2 seconds duration and it will have 640 points (a $10Hz$ frequency wave sampled at 32 times oversampling factor will have $2 \times 32 \times 10 = 640$ samples in 2 seconds of the record). Since our input signal is periodic, we can safely use $N=256$ point FFT, anyways the FFT will extend the signal when computing the FFT (see additional topic on spectral leakage that explains this extension concept).

Due to Matlab’s index starting at 1, the DC component of the FFT decomposition is present at index 1.

>>X(1)
 1.1762e-14   (approximately equal to zero)

That’s pretty easy. Note that the index for the raw FFT are integers from $1 \rightarrow N$. We need to process it to convert these integers to $frequencies$. That is where the $sampling$ frequency counts.

Each point/bin in the FFT output array is spaced by the frequency resolution $\Delta f$ that is calculated as

\[ \Delta f = \frac{f_s}{N} \]

where, $f_s$ is the sampling frequency and $N$ is the FFT size that is considered. Thus, for our example, each point in the array is spaced by the frequency resolution

\[\Delta f = \frac{f_s}{N} = \frac{32*f_c}{256} = \frac{320}{256} = 1.25 Hz \]

Now, the $10 Hz$ cosine signal will leave a spike at the 8th sample ($10/1.25=8$), which is located at index 9 (See next figure).

>> abs(X(8:10)) %display samples 7 to 9
ans = 0.0000  128.0000    0.0000

Therefore, from the frequency resolution, the entire frequency axis can be computed as

%calculate frequency bins with FFT
df=fs/N %frequency resolution
sampleIndex = 0:N-1; %raw index for FFT plot
f=sampleIndex*df; %x-axis index converted to frequencies

Now we can plot the absolute value of the FFT against frequencies as

subplot(3,1,2); stem(sampleIndex,abs(X)); %sample values on x-axis
title('X[k]'); xlabel('k'); ylabel('|X(k)|');
subplot(3,1,3); plot(f,abs(X)); %x-axis represent frequencies
title('X[k]'); xlabel('frequencies (f)'); ylabel('|X(k)|');

The following plot shows the frequency axis and the sample index as it is for the complex FFT output.

Figure 4: Magnitude response from FFT output plotted against sample index (top) and computed frequencies (bottom)

After the frequency axis is properly transformed with respect to the sampling frequency, we note that the cosine signal has registered a spike at $10 Hz$. In addition to that, it has also registered a spike at $256-8=248^{th}$ sample that belongs to negative frequency portion. Since we know the nature of the signal, we can optionally ignore the negative frequencies. The sample at the Nyquist frequency ($f_s/2$) mark the boundary between the positive and negative frequencies.

>> nyquistIndex=N/2+1;
>> X(nyquistIndex-2:nyquistIndex+2).'
ans =
1.0e-13 *
  -0.2428 + 0.0404i
  -0.1897 + 0.0999i
  -0.3784          
  -0.1897 - 0.0999i
  -0.2428 - 0.0404i

Note that the complex numbers surrounding the Nyquist index are complex conjugates and they represent positive and negative frequencies respectively.

FFTShift

From the plot we see that the frequency axis starts with DC, followed by positive frequency terms which is in turn followed by the negative frequency terms. To introduce proper order in the x-axis, one can use FFTshift function Matlab, which arranges the frequencies in order: negative frequencies $\rightarrow$ DC $\rightarrow$ positive frequencies. The fftshift function need to be carefully used when $N$ is odd.

For even N, the original order returned by FFT is as follows (note: all indices below corresponds to Matlab’s index)

● $X[1]$ represents DC frequency component
● $X[2]$ to $X[N/2]$ terms are positive frequency components
● $X[N/2+1]$ is the Nyquist frequency ($F_s/2$) that is common to both positive and negative frequencies. We will consider it as part of negative frequencies to have the same equivalence with the fftshift function.
● $X[N/2+1]$ to $X[N]$ terms are considered as negative frequency components

FFTshift shifts the DC component to the center of the spectrum. It is important to remember that the Nyquist frequency at the (N/2+1)th Matlab index is common to both positive and negative frequency sides. FFTshift command puts the Nyquist frequency in the negative frequency side. This is captured in the following illustration.

Figure 5: Role of FFTShift in ordering the frequencies

Therefore, when $N$ is even, ordered frequency axis is set as

\[f = \Delta f \left[ -\frac{N}{2}:1:\frac{N}{2}-1 \right] = \frac{f_s}{N} \left[ -\frac{N}{2}:1:\frac{N}{2}-1 \right]\]

When (N) is odd, the ordered frequency axis should be set as

\[ f = \Delta f \left[ -\frac{N-1}{2}:1:\frac{N+1}{2}-1 \right] = \frac{f_s}{N} \left[ -\frac{N-1}{2}:1:\frac{N-1}{2}-1 \right]\]

The following code snippet, computes the fftshift using both the manual method and using the Matlab’s in-build command. The results are plotted by superimposing them on each other. The plot shows that both the manual method and fftshift method are in good agreement.

%two-sided FFT with negative frequencies on left and positive frequencies on right
%following works only if N is even, for odd N see equation above
X1 = [(X(N/2+1:N)) X(1:N/2)]; %order frequencies without using fftShift
X2 = fftshift(X);%order frequencies by using fftshift

df=fs/N; %frequency resolution
sampleIndex = -N/2:N/2-1; %raw index for FFT plot
f=sampleIndex*df; %x-axis index converted to frequencies
%plot ordered spectrum using the two methods
figure;subplot(2,1,1);stem(sampleIndex,abs(X1));hold on; stem(sampleIndex,abs(X2),'r') %sample index on x-axis
title('Frequency Domain'); xlabel('k'); ylabel('|X(k)|');%x-axis represent sample index
subplot(2,1,2);stem(f,abs(X1)); stem(f,abs(X2),'r') %x-axis represent frequencies
xlabel('frequencies (f)'); ylabel('|X(k)|');

Figure 6: Magnitude response of FFT result after applying FFTShift : plotted against sample index (top) and against computed frequencies (bottom)

Comparing the bottom figures in the Figure 4 and Figure 6, we see that the ordered frequency axis is more meaningful to interpret.

IFFTShift

One can undo the effect of fftshift by employing ifftshift function. The ifftshift function restores the raw frequency order. If the FFT output is ordered using fftshift function, then one must restore the frequency components back to original order BEFORE taking IFFT. Following statements are equivalent.

X = fft(x,N) %compute X[k]
x = ifft(X,N) %compute x[n]

X = fftshift(fft(x,N)); %take FFT and rearrange frequency order (this is mainly done for interpretation)
x = ifft(ifftshift(X),N)% restore raw frequency order and then take IFFT

Some observations on FFTShift and IFFTShift

When $N$ is odd and for an arbitrary sequence, the fftshift and ifftshift functions will produce different results. However, when they are used in tandem, it restores the original sequence.

>>  x=[0,1,2,3,4,5,6,7,8]
0     1     2     3     4     5     6     7     8
>>  fftshift(x)
5     6     7     8     0     1     2     3     4
>>  ifftshift(x)
4     5     6     7     8     0     1     2     3
>>  ifftshift(fftshift(x))
0     1     2     3     4     5     6     7     8
>>  fftshift(ifftshift(x))
0     1     2     3     4     5     6     7     8

When $N$ is even and for an arbitrary sequence, the fftshift and ifftshift functions will produce the same result. When they are used in tandem, it restores the original sequence.

>>  x=[0,1,2,3,4,5,6,7]
0     1     2     3     4     5     6     7
>>  fftshift(x)
4     5     6     7     0     1     2     3
>>  ifftshift(x)
4     5     6     7     0     1     2     3
>>  ifftshift(fftshift(x))
0     1     2     3     4     5     6     7
>>  fftshift(ifftshift(x))
0     1     2     3     4     5     6     7

For Python code, please check the following book: Digital Modulations using Python – by Mathuranathan Viswanathan

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

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.

Spectral leakage

Equivalent Noise Bandwidth:[2]

Coherent power gain:

Scalloping loss

Figure of merits for window functions

List of window functions

Boxcar (Rectangular)

Barthann

Bartlett

Blackman

Blackman-Harris

Bohman

Cosine

Flattop

Hamming

Hann

Nuttall

Parzen

Triangular

References:

Similar articles

Additional Resources:

Books by the author

FFT and spectral leakage

Equivalent noise bandwidth

Python script

Reference

Similar articles

Books by the author

Deriving Parseval’s theorem

Reference

Books by the author

Outline

Discrete-time domain representation

Represent the signal in frequency domain using FFT

3a. Extract amplitude of frequency components (amplitude spectrum)

3b. Extract phase of frequency components (phase spectrum)

Reconstruct the time domain signal from the frequency domain samples

Topics in this chapter

Books by the author

Four types of Fourier Transforms:

Real and Complex versions of transforms

Real DFT:

Complex DFT:

Fast Fourier Transform (FFT)

Interpreting the FFT results

FFTShift

IFFTShift

Some observations on FFTShift and IFFTShift

Topics in this chapter

Books by the author

Equivalent Noise Bandwidth:^[2]