Interpret FFT results – obtaining magnitude and phase information

by

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
Cosine wave with phase shift

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

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
Noisy phase spectrum

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
Phase spectrum after fixing rounding off errors

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
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: PoorBelow averageAverageGoodExcellent (158 votes, average: 4.37 out of 5)

Topics in this chapter

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

Books by the author

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

33 thoughts on “Interpret FFT results – obtaining magnitude and phase information”

  2. Hi Mathuranathan,
    Thanks for such an amazing article. It is very helpful in interpreting the data and understanding the Fourier Transform. Would you please help me interpreting the same for a 2D Fourier transform? Or can you please share any articles related to the 2D FFT or fft2(). It would be of great help.
    Thanks again for such a vivid explanation of fft function.

    Reply

  3. The phase noise was driving me nuts! Already suspected something fishy with rounding errors but to no avail. The atan2 function did the trick for me! Thanks a lot!

    Reply

  4. While doing the phase calculation from the above equation why i am getting more values than the input? I need that in same length.

    Reply

    • The length of the phase vector (derived from FFT output) depends on the FFT length (N). In the example above the FFT length is set to 256. So it gives 256 phase values corresponding to 256 frequency bins.

      Reply

  5. This example worked well until I changed the signal frequency and found the computed phase only to be correct when the sampling frequency is exactly a power of 2 higher than the signal frequency(x4,8,16,32,64 etc) otherwise it not correct. Is there some other correction that needs to be performed to the result in these cases?

    Reply

    • This is due to spectral leakage phenomenon. You may read about it here
      https://www.gaussianwaves.com/2011/01/fft-and-spectral-leakage-2/

      Due to spectral leakage, you will see non-zero values for frequency components that are near the frequency of the signal. This will also reflect in the phase spectrum.

      Set Fs=100 and N=length(x) (length(x)= 640 and it is not a power of 2). Now you will observe perfect frequency spectrum. However, if you change N=power of 2, spectrum will be affected by spectral leakage.

      Reply

  6. I loved this explanation. Really simple overview of FFT results. Plus, learned about phase, and all the tricks to view it properly. Beautifully done.

    Reply

  7. Anyone knows how to change Matlab code into python? For the part where you apply the threshold thingy. I have tried doing it before but it doesn’t work

    Reply

  11. Hi Mathuranathan, thank you for this article.
    I have a doubt regarding calculation of DFT in matlab using DFT formulae and using inbulit MATLAB function FFT.
    The
    phase of DFT computed using DFT formula and FFT (inbuilt MATLAB
    function) is different at N/2. Why is it so? Attaching sample code
    clc;
    close all;
    clear all;

    x = [2 3 -1 4];
    N = length(x);
    X = zeros(4,1)

    %DFT formulae
    for k = 0:N-1
    for n = 0:N-1
    X(k+1) = X(k+1) + x(n+1)*exp(-j*pi/2*n*k)
    end
    end

    t = 0:N-1
    subplot(311)
    stem(t,x);
    xlabel(‘Time (s)’);
    ylabel(‘Amplitude’);
    title(‘Time domain – Input sequence’)

    %Magniltude Plot of DFT
    magnitude=abs(X)
    subplot(312);
    stem(0:N-1,magnitude);
    xlabel(‘Frequency’);
    ylabel(‘|X(k)|’);
    title(‘Magnitude Response’);

    %Phase Plot of DFT
    phase=atan2(imag(X),real(X))*180/pi
    subplot(313)
    stem(0:N-1,phase)
    xlabel(‘Frequency’);
    ylabel(‘Phase’);
    title(‘Frequency domain – Phase response’)

    %fft
    figure;

    X_k=fft(x,N)

    t=0:N-1;
    subplot(311)
    stem(t,x);
    xlabel(‘Time (s)’);
    ylabel(‘Amplitude’);
    title(‘Input sequence’)

    %Magniltude Plot of DFT
    magnitude=abs(X_k)
    subplot(312);
    stem(0:N-1,magnitude);
    xlabel(‘Frequency’);
    ylabel(‘|X(k)|’);
    title(‘Magnitude Response’);

    %Phase Plot of DFT
    phase=atan2(imag(X_k),real(X_k))*180/pi
    subplot(313);
    stem(0:N-1,phase);
    xlabel(‘Frequency’);
    ylabel(‘Phase’);
    title(‘Phase Response’);

    Reply

  14. Hello Mathuranathan
    Thank you very much for this article, I’ve a naive question, why when I try with a sine wave
    x=A*sin(2*pi*fc*t+phi);
    I get phase spectrum peak at 60 degree (pi/3) instead of 30 degree.
    I’m not sure why this is happening do you have any clue?
    Thank you

    Reply

  18. Hi Mathuranathan, thank you for this article, it’s very helpful. I’d like to extract the phase information from an OFDM signal. How can I do this? Thanks in advance

    Reply

    • In OFDM, the phase estimation should be done after the FFT block in the receiver. Several algorithms/technique available for phase estimation in an OFDM system.

      Consider a coherent detection system.
      1. For each transmitted OFDM symbol, the FFT output contains say N modulated values (say QAM modulation is used)
      2. These values may contain random phase shifts and amplitude variations caused by local oscillator drift, jitter, channel response and other factors.
      3. We usually use a channel estimator (can be a pilot/preamble based one) to estimate the reference phases and amplitudes.
      4. Once the reference phases and amplitudes are learnt, we can apply them to correct phase/amplitude errors in the actual data transmission.

      You can take a look at this patent. I hope it offers more clue on how to use it for your system.
      https://www.google.com/patents/US7423960

      Reply

      • Thank you for your answer. In my system, an OFDM signal is used for extracting the distance between a TX and RX. The signal is composed solely by zadoff-chu pilots. I have extracted a coarse distance measure from the correlation function between the received signal and a refrence signal. Now I’d like to perform a fine estimation exploiting the signal phase estimated in frequency domain but I’m not able to do this.

        Reply

  19. While reconstructing the signal, what happens to the phase information. Since that is not being passed to the ifft(), how is the phase information restored in the time domain?

    Reply

Post your valuable comments !!!