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

$\endgroup$
  • $\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$ – Dilip Sarwate Apr 7 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$ – ZaellixA Apr 7 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$ – Cedron Dawg Apr 7 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$ – ZaellixA Apr 7 at 14:57
  • $\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$ – Cedron Dawg Apr 7 at 15:54
0
$\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$.

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

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.