3
\$\begingroup\$

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

Sinusoid signal with power spectral density plots

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) )
\$\endgroup\$
2
\$\begingroup\$

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.
\$\endgroup\$
1
\$\begingroup\$

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.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.