# Spectrum Analysis with Discrete Fourier Transform

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


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.

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.
err = mean(abs(cdft - bfft));