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 obtain magnitude and phase information from the FFT results.

## Outline

In this discussion, I will 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

## 1. 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.

1 2 3 4 5 6 7 8 9 10 |
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.

1 2 |
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.

1 2 3 4 5 |
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

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

1 2 3 4 5 6 7 |
>> 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.

1 2 3 4 5 6 |
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.

## 4. Reconstruct the time domain signal from the frequency domain samples

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

1 2 3 |
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 (~ 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 article: **

**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 (this article)

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

Constructing the Auto Correlation Matrix using FFT