10
$\begingroup$

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.

enter image description here

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");
$\endgroup$
1
  • $\begingroup$ Hi! Pure implementation of DFT shall take quadratic time. Probably there's a testing error... $\endgroup$
    – Fat32
    Commented Aug 26, 2018 at 0:17

2 Answers 2

9
$\begingroup$

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:

enter image description here

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 (Look at the SignalProcessing\Q51516 folder).

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.

$\endgroup$
10
  • $\begingroup$ Really nice answer and tips! It's great to have you sharing the code in GitHub. $\endgroup$
    – David
    Commented Aug 26, 2018 at 5:29
  • $\begingroup$ Thank You. I really believe that in our days code is almost a "Must" in any "Real World" question. $\endgroup$
    – Royi
    Commented Aug 26, 2018 at 5:46
  • $\begingroup$ Wow, those are some really great points. Thank you for such a comprehensive answer, I'll keep them in mind from now. Yes, planning to look at more MATLAB code and learn. Thanks again. $\endgroup$ Commented Aug 26, 2018 at 9:59
  • $\begingroup$ @Royi he might have too little reputation to upvote yet. I'll do that for him ;) Abhinav: You should probably learn very soon that in Matlab, doing things in loops is incredibly much slower than doing it using vector/matrix operations: thus, your "naive" DFT shouldn't be a for loop, but a single multiplication of your input vector with the DFT matrix (which you can calculate once before). From that, the quadratic behaviour of complexity becomes pretty apparent! $\endgroup$ Commented Aug 26, 2018 at 10:54
  • 1
    $\begingroup$ got it upvoted! $\endgroup$ Commented Aug 26, 2018 at 12:25
2
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ Do you think LightSpeed is still relevant for modern MATLAB versions? $\endgroup$
    – Mark
    Commented May 16, 2021 at 16:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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