4
$\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$
  • $\begingroup$ Hi! Pure implementation of DFT shall take quadratic time. Probably there's a testing error... $\endgroup$ – Fat32 Aug 26 '18 at 0:17
8
$\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.

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$
  • $\begingroup$ Really nice answer and tips! It's great to have you sharing the code in GitHub. $\endgroup$ – David Aug 26 '18 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 Aug 26 '18 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$ – Abhinav Jain Aug 26 '18 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$ – Marcus Müller Aug 26 '18 at 10:54
  • 1
    $\begingroup$ got it upvoted! $\endgroup$ – Marcus Müller Aug 26 '18 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$

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.