Spectrum Analysis with Discrete Fourier Transform

Question

Here is a simple implementation of the Discrete Fourier Transform:

myFourierTransform.m

function [ft] = myFourierTransform (X, n)
% Objective: 
%   Apply the Discrete Fourier Transform on X.
% Input:
%   X - 1xM - complex vector - data points (signal discretisation). 
%   n - 1x1 - integer scalar - number of discrete frequencies (spectrum discretisation). 
% Output:
%   ft - 1xn - complex vector - "quantity of" i-th frequency in X.
% Complexity: O(n * M).
% Author:  
% Date:  

  ft = zeros(1, n);
  freq = @(point_index, freq_index, total_freqs)...
  exp( -2 .* pi .* j  .*(point_index - 1) .* (freq_index - 1) / total_freqs );  

  for i = 1 : n

    f_i = freq(1: numel(X), i, n);

    ft(i) = dot(X, f_i);

  end

end  

Here is a script file applying the above and comparing it to the built-in function fft(X, n):

applyMyFourierTransform.m

clc; clear; close all;

% Signal parameters
A = 1;                                 % signal amplitude.
f = 10;                                % frequency of the signal 10 Hz.
theta = 0.25 * pi;                     % initial phase of the signal.     
measurement_rate = f * 30;             % Nyquist-Shannon Th.-rate of measurement >=2*f. 

nPeriods = 7;                          % number of periods to be displayed.
tBeg = 0;                                                               
tEnd = nPeriods * 1 / f;                                       
tStep = 1 / measurement_rate ;  

t = tBeg : tStep : tEnd;                         
x = A * sin(2 * pi * f * t + theta);   % sine wave A= 1 V,f= 10 Hz, phi_0= pi/4 rad.
P = A^2 / 2;                           % signal power=\frac{1}{T}\int_0^T {x{^2}(t)dt}.
point_num = numel(x);

% Plot the test signal.
subplot(3, 1, 1)
plot(t, x);
xlabel('Time [sec]');
ylabel('Amplitude');
title('Test Signal.');
legend( sprintf( 'f = %.2f [Hz] \n A = %.2f [V]', f, A ) );

% Apply FFT and obtain spectrum of test signal.
freq_num = 2^10;                       % Discretise the spectrum in 1024 frequencies.
bfft = fft(x, freq_num);               % built-in Discrete Fourier Transform.

bfft = fftshift(bfft);                 % shift f(0) at the middle of the vector.            

% Calculate signal power \frac{X_k * X_k^{*}}{Frequency_Discretisation * Data_size}.                  
power = bfft .* conj(bfft) / (freq_num * point_num); 

% Frequency axis starts at: -f/2 ends at: f/2 and includes: n frequencies.  
frequencies = measurement_rate * (-freq_num / 2 : freq_num/ 2 -1) / freq_num; 

subplot(3, 1, 2)
plot(frequencies, power);
xlabel('Frequency [Hz]');
axis([-10*f 10*f])
set(gca, 'XTick', [-10*f : f : 10*f]);
ylabel('Power [Watt]');
title('Power Spectral Density.')
legend( sprintf( 'f = %.2f [Hz]\n P = %.2f [W].', f ,P ) )

% Apply custom implementation of Discrete Fourier Transform.
cdft = myFourierTransform(x, freq_num);
cdft = fftshift(cdft);      
cpower = cdft .* conj(cdft) / (freq_num * point_num); 

subplot(3, 1, 3)
plot(frequencies, cpower);
xlabel('Frequency [Hz]');
axis([-10*f 10*f])
set(gca, 'XTick', [-10*f : f : 10*f]);
ylabel('Power [Watt]');
title('Power Spectral Density.')
legend( sprintf( 'f = %.2f [Hz]\n P = %.2f [W].', f, P ) )

error = abs( sum(cdft) - sum(bfft) )

Command Line:

>> applyMyFourierTransform

Result:

>> error =   1.7346e-011

Questions

How can I estimate the validity of the values produced by the custom implementation? Is this a reasonable estimate?

error = abs( sum(cdft) - sum(bfft) )

Answer 1

A few things:

When you have the start, endpoint and the number of steps you want in a vector then linspace is better than colon:

t = linspace(tBeg, tEnd, measurement_rate);

It's cleaner, since you don't have to calculate the length of each step.

---

exp( -2 .* pi .* j  .*(point_index - 1) .* (freq_index - 1) / total_freqs );

j looks like a variable, especially since you use i as a variable later on. The recommended way to denote the imaginary unit is 1i or 1j. That way it can't be confused with a variable, and you don't risk overwriting it.

You don't need .* when multiplying by scalars. It's only needed between matrices or vectors. If you only use it between matrices then it's easy to see which variables are vectors/matrices, and which are just scalars. It's now hard to tell if point_index, freq_index and/or total_freqs are scalars or matrices.

There's only one vector in that statement (point_index), so no dots are needed.

The for loop can be vectorized, which is the MATLAB way to do something like that. I don't have MATLAB at hand, so I can't do it for you, unfortunately. If you have a version newer than R2015B, then you can do implicit matrix expansion by adding two vectors of different orientation:

[1; 2; 3] + [4, 5, 6] = [5, 6, 7;6, 7, 8;7, 8, 9]

You can do this with point_index and freq_index, then do matrix multiplication to get the final ft.

---

Don't use error as a variable name. It's a useful function that you can use to, well, throw errors.

---

You don't need a space after \n in a string, unless you actually want a space in front of the first letter in the next line:

f = %.2f [Hz] \nA = %.2f [V]'

---

• Your variable names are very good
• Your coding style is good
• I don't like clear and close all in the beginning of the script, but that's up to you.

Answer 2

Two things to add to Stewie's answer:

1. The second statement in your DFT function wraps using .... I would recommend that you indent the second line of this statement to make that more obvious.

2. Regarding your computation of the error: positive and negative errors will cancel out the way you do this. It would be better to compute the mean of absolute differences, or the mean of their squares:

   err = mean(abs(cdft - bfft));
   err = mean(abs(cdft - bfft).^2);
   

   These are the "mean absolute error" and "mean square error" metrics, very common in all fields of engineering.

   I've also seen people compare the power spectra to evaluate FFT algorithms.