1
votes

I want to understand the discrete fourier transform by implementing it by myself. While the result returned by my DFT is not correct the in matlab included version returns the correct frequencies of the original signal. So the question is, where went I wrong. Is it a math or a implementation problem?

%% Initialisation
samples=2000;
nfft = 1024;
K = nfft / 2 + 1;
c = 264;
e = 330;
t = -1:1/samples:1-1/samples;
[~, N] = size(t);
f = (sin(2*c*pi*t)+cos(2*e*pi*t)).*exp(-pi*(2*t-1).^2);
X = zeros(nfft, 1);

%% Discrete Fourier Transform
if true
    for k=1:nfft
        for n=1:nfft
            X(k) = X(k) + f(n)*exp(-j*2*pi*(k-1)*(n-1)/N);
        end
    end
else
    X=fft(f, nfft);
end
R = abs(X(1:K));
[V,I] = sort(R,'descend');
F1 = samples*(I(1)-1)/nfft;
F2 = samples*(I(2)-1)/nfft;
disp(F1)
disp(F2)
plot(1:K, R, 1:K, real(X(1:K)), 1:K, imag(X(1:K)))
1

1 Answers

3
votes

The issue lies in the number of samples for which the transform is done.

Xall = fft(f);
plot(abs(Xall(1:500)),'b');
hold on
plot(abs(X(1:500)),'r');

What you compute matches the result from the FFT done on all samples (i.e. with 4000 real samples in and 4000 complex values out).

Now, if you read the documentation of FFT with doc fft you will see that the signal is truncated if the output size is smaller than the input size. If you try:

Y = zeros(nfft, 1);
for k=1:nfft
    for n=1:nfft
        Y(k) = Y(k) + f(n)*exp(-1j*2*pi*(k-1)*(n-1)/nfft);
    end
end
Y2 = fft(f(:),nfft); %make it a column
abs(sum(Y-Y2)) %6.0380e-12 , result within precision of the double float format