FFT vs DFT Run Time Comparison (Complexity Analysis) in MATLAB

Question

So we were given an assignment to plot the time taken by the FFT algorithm by MATLAB and a DFT algorithm written by me in MATLAB.G

The expected output should have been the DFT algorithm following the O(n^2) complexity, however when I plot the output it was coming out to be linear.

Is there some problem with my implementation of DFT or some other reason that is causing it?

The first graph is time taken to compute dft vs N(N point fft) whereas the other is for my implementation of DFT vs N.

timedft = zeros(1,length(2:1024));
timefft = zeros(1,length(2:1024));
%DFT
sig = rand(1,1024);
for N = 2:1024
    tic;
    opdft = zeros(1,N);
    for i = 0:N-1
        for l = 0:length(sig)-1
            temp = sig(1,l+1).*exp(-2*pi*l*i*j/N);
            opdft(1,i+1)=opdft(1,i+1)+temp;
        end
    end
    timedft(1,N-1) = toc;
    tic;
    fft(rand,N);
    timefft(1,N-1) = toc;
end
figure;
subplot(211);
plot(2:1024,timedft);
title("fftvsdft");
xlabel("N");
ylabel("Time")
subplot(212);
plot(2:1024,timefft);
xlabel("N");
ylabel("Time");

Answer 1

Abhinav Jain, Welcome to DSP Community.

I build for you a proper testing of the run time comparison.

Few tips about timing in MATLAB:

1. Never time in a script. Always call a function to do the heavy lifting. When you run something from script it runs in the global scope which means MATLAB can't optimize it as well as it could if it was in a function.
2. When you time a function, make few iterations of the measurement. I like taking the median of this iterations. Others like the minimum. Usually the mean isn't good as it is sensitive to outliers.
3. MATLAB is a JIT based Script language. It means the first run of anything takes more time. Hence (2) is crucial.
4. Pay attetion that the result of the FFT / DFT is complex. Hence when you allocate memory for a complex array you should use - vArrayName = complex(zeros(arrayLength, 1));.
5. When you need the imaginary number $ i = \sqrt{-1} $ in MATLAB you should use 1i or 1j (No need for multiplication). Of course you can use 5i as well. Yet no need to multiply as it assists the JIT engine to understand you mean the imaginary number and not a variable.

Here are the results:

As you can see, indeed the DFT Implementation has running time which behaves as O(n^2).

The code goes like:

%% Simulation Parameters

vNumSamples     = 2:2:1024;
numIterations   = 6;


%% Generate Data

mDftTime = zeros(numIterations, length(vNumSamples));
mFftTime = zeros(numIterations, length(vNumSamples));

for jj = 1:length(vNumSamples)
    numSamples = vNumSamples(jj);
    vX = randn(numSamples, 1);

    for ii = 1:numIterations
        hDftTimer           = tic();
        vXDft               = ApplyDft(vX, numSamples);
        mDftTime(ii, jj)    = toc(hDftTimer);

        hFftTimer           = tic();
        vXFft               = fft(vX);
        mFftTime(ii, jj)    = toc(hFftTimer);
    end

end


%% Run Time Analysis

vDftMedian = median(mDftTime).';
vFftMedian = median(mFftTime).';

vDftMean = mean(mDftTime).';
vFftMean = mean(mFftTime).';

vDftMax = max(mDftTime).';
vFftMax = max(mFftTime).';

vDftMin = min(mDftTime).';
vFftMin = min(mFftTime).';


%% Display Results

figureIdx = figureIdx + 1;

hFigure = figure('Position', figPosLarge);
hAxes   = subplot(4, 1, 1);
% set(hAxes, 'NextPlot', 'add');
hLineSeries = plot(vNumSamples, [vDftMedian, vFftMedian]);
set(hLineSeries, 'LineWidth', lineWidthNormal);
% set(hLineSeries(2), 'LineStyle', ':');
set(get(hAxes, 'Title'), 'String', {['DFT vs. FFT Run Time - Median']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Input Size']}, ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', {['Run Time [Sec]']}, ...
    'FontSize', fontSizeAxis);
hLegend = ClickableLegend({['DFT'], ['FFT']});

hAxes   = subplot(4, 1, 2);
% set(hAxes, 'NextPlot', 'add');
hLineSeries = plot(vNumSamples, [vDftMean, vFftMean]);
set(hLineSeries, 'LineWidth', lineWidthNormal);
% set(hLineSeries(2), 'LineStyle', ':');
set(get(hAxes, 'Title'), 'String', {['DFT vs. FFT Run Time - Mean']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Input Size']}, ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', {['Run Time [Sec]']}, ...
    'FontSize', fontSizeAxis);
hLegend = ClickableLegend({['DFT'], ['FFT']});

hAxes   = subplot(4, 1, 3);
% set(hAxes, 'NextPlot', 'add');
hLineSeries = plot(vNumSamples, [vDftMax, vFftMax]);
set(hLineSeries, 'LineWidth', lineWidthNormal);
% set(hLineSeries(2), 'LineStyle', ':');
set(get(hAxes, 'Title'), 'String', {['DFT vs. FFT Run Time - Max']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Input Size']}, ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', {['Run Time [Sec]']}, ...
    'FontSize', fontSizeAxis);
hLegend = ClickableLegend({['DFT'], ['FFT']});

hAxes   = subplot(4, 1, 4);
% set(hAxes, 'NextPlot', 'add');
hLineSeries = plot(vNumSamples, [vDftMin, vFftMin]);
set(hLineSeries, 'LineWidth', lineWidthNormal);
% set(hLineSeries(2), 'LineStyle', ':');
set(get(hAxes, 'Title'), 'String', {['DFT vs. FFT Run Time - Min']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Input Size']}, ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', {['Run Time [Sec]']}, ...
    'FontSize', fontSizeAxis);
hLegend = ClickableLegend({['DFT'], ['FFT']});

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

As you can see in the code above I run numIterations = 6 each algorithm for each number of samples.
Then, on the time array mDftTime / mFftTime I can analyze the Median / Mean / Min / Max run time.
This is the proper way to analyze function run time in MATLAB.

The full code is available on my StackExchange Signal Processing Q51516 GitHub Repository.

Remark
You're doing your first steps in DSP world and big part of being modern engineer is about programming code.
You should invest time in learning to properly write code and style in a way it is easy to under stand and manage for you and your partners.
Try learning from other people examples.

Answer 2

Most FFT implementations use a set up step where the complex exponential terms are pre computed and that step is omitted in the legacy complexity calculations.

MATLAB’s interpreter typically runs slower the first time an mfile runs.

MATLAB has been using fftw for its fft routine for a number of years. It probably used fft pack prior to using fftw. fftw has its own prior initializations as well called “wisdom” which should not be part of your bench marking.

You can access the fftw planner with a fftw() call.

So essentially, the initializations need to be eliminated.

you can use dftmtx to precompute the exponentials and eliminate your for loop.

On a deeper level, a modern multithreaded processor’s performance is a lot harder to predict run time based on flops counts. MATLAB originally would return flops but that capability was removed.

There is Lightspeed

https://github.com/tminka/lightspeed

that returns some flops capability. You need to compile it from source. This might help using flops in your benchmarks. MATLAB also has a nice profiler that you can use to tune your code. I'm guessing you will see the majority of time is calculating exp().

Even if you account for the initializations, you may not get the results you expect. I can vouch that on a 8086 theory will agree with measurement as long as you don’t run out of memory.

fft works with both double and single precision numbers. single is faster on a number of processors.