0
$\begingroup$

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

Improve this question
$\endgroup$
6
  • $\begingroup$ 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! $\endgroup$ Commented Apr 7, 2019 at 13:34
  • $\begingroup$ You do have a point here @DilipSarwate. Didn't think of it that way. Point taken and thanks for the tolerance. $\endgroup$ Commented Apr 7, 2019 at 14:42
  • $\begingroup$ 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. dsp.stackexchange.com/questions/54848/… Here is an article I wrote discussing an interpretation of the meaning of a DFT: dsprelated.com/showarticle/768.php $\endgroup$ Commented Apr 7, 2019 at 14:49
  • $\begingroup$ 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. $\endgroup$ Commented Apr 7, 2019 at 14:57
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Apr 7, 2019 at 15:54

1 Answer 1

Reset to default
1
$\begingroup$

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$.

Improve this answer
$\endgroup$
2
  • $\begingroup$ That was it! It slipped right past me. Thanks a lot for the quick and clean answer. $\endgroup$ Commented Apr 7, 2019 at 14:48
  • $\begingroup$ @ZaellixA your welcome... $\endgroup$ Commented Apr 7, 2019 at 20:14