0

I have implemented a (direct) DFT in MatLab (following script) and compared it to the built-in FFT routine. The magnitude response seems to be identical (excluding some possibly round-off errors), but in the phase I get different results.

I can't find the source of these differences and I would like some insight on that.

script (only the case of even vector length is presented here):

%% Calculate needed variables
DFT = zeros(length(data), 1); % Pre-allocate result vector
n = (1:length(data))/length(data); % Calculate the n/N factor
n = -1j * n * 2 * pi; % Calculate the (-j * n * 2 * pi)/N factor

dftLength = (length(data)/2) + 1; % Calculate the length of unique Fourier coefficients

% Calculate the coefficients
for i = 1:dftLength
    exponent = exp((i - 1) * n); % Calculate the exponent
    DFT(i) = exponent * data; % Multiply the signal with the exponent

    % Add the conjugate symmetric part of the spectrum
    if i ~= 1 
        % Skip first bin (it's DC)
        DFT(2 * dftLength - i) = conj(DFT(i)); % Add conjugate symmetric bin
    end
end

You can see the resulting plot here

Phase of the DFT and FFT

ZaellixA
  • 1,286
  • 2
  • 9
  • 18
  • I would vote to close this question on the grounds that is it purely a programming question and therefore of-topic for this site, but there are enough people who think that DSP is a subset of MATLAB and so will vehemently dispute the close vote as inappropriate. So, have fun, folks! – Dilip Sarwate Apr 07 '19 at 13:34
  • You do have a point here @DilipSarwate. Didn't think of it that way. Point taken and thanks for the tolerance. – ZaellixA Apr 07 '19 at 14:42
  • This isn't correct: "exponent = exp((i - 1) * n);" in any way. A pure DFT calculation should be a set of nested loops. The outer one for the bin, the inner one for the dot product of the basis vector times your signal. Yours looks nothing like that. Check my answer here for a single DFT bin calculation in Python.

    https://dsp.stackexchange.com/questions/54848/measure-sine-wave-amplitude-from-adc-signal/54879#54879

    Here is an article I wrote discussing an interpretation of the meaning of a DFT: https://www.dsprelated.com/showarticle/768.php

     – Cedron Dawg Apr 07 '19 at 14:49
  • Well, @CedronDawg thanks for the answer but as far as I understand the code, I have the "outer" loop implemented for the bins (like you suggest) and I take advantage of the vectorisation that MatLab does to avoid the inner loop that reads through the time-domain samples (in practice performing vector multiplication to implement the inner loop). Please correct me if I am wrong but after Fat32's answer the results are the same as MatLab's fft built-in routine. – ZaellixA Apr 07 '19 at 14:57
  • 1
    Sorry, I shouldn't have said anything as I don't use Matlab. There is the one based arrays vs zero based indexing issue, and all the implicit loops that I tend not to see. If it matches, then you got it right. I didn't look closely enough. – Cedron Dawg Apr 07 '19 at 15:54
  • No worries @CedronDawg. It's good that you spent your time trying to help. No worries at all :). Edit: Plus, you provided some good links over there (thumbs-up). – ZaellixA Apr 07 '19 at 17:01

1 Answers1

1

you should replace the line

n = (1:length(data))/length(data); % Calculate the n/N factor

with the line

n = (0:length(data)-1)/length(data); % Calculate the n/N factor

Note that in the standard DFT (and FFT), the indices for both $x[n]$, $X[k]$ range in $0 \leq n,k \leq N-1$.

Fat32
  • 28,152
  • 3
  • 24
  • 50