# different solutions in matlab / octave using dft and fft

i am new here in dsp.stackexchange and I am trying to do my first basic steps with fourier-transformation. Some years ago I learned the basic theory in university and also developed a fft implementation in matlab. Now I try to get back into the topic.

I started by reading the mathematic theory again and tried a dft example implementation with octave which i found here: http://www.bogotobogo.com/Matlab/Matlab_Tutorial_DFT_Discrete_Fourier_Transform.php

This example worked fine for the data array x = [2 3 -1 4] (see "example 1" in code). If you look into the below picture you see the same results for the DFT implementation (left plots) and for the build-in FFT implementation of octave (right plots).

But if I now use a longer data array like for example x = [2 3 -1 4 2 3 -1 4] (see "example 2" within the code), I get different solutions for DFT and FFT.

By debugging the code i found out, that the solutions are not that much different. If you compare the result of the DFT and FFT you see, that there are the same numbers with the only difference that in the FFT solution there is always a zero in between the values.

I don't really understand what is happening here. Can anyone maybe help me understanding this behavior or maybe give me a hint of where I have to search for the solution? I searched google a lot and found some topics about "zero padding" before using the fft function. But to be honest I do not really understand the advantage of zero padding and if this might have something to do with my problem.

Thank you very much in advance. anon1234

% example 1
%x = [2 3 -1 4];

% example 2
x = [2 3 -1 4 2 3 -1 4];

% example 3
%Fs = 1000;            % Sampling frequency
%T = 1/Fs;             % Sampling period
%L = 1500;             % Length of signal
%t = ((0:L-1)*T);     % Time vector
%x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);

% solution by build-in fft function
XFFT = (fft(x)');

% solution by two different dft functions
N = length(x);
X = zeros(N,1);
for k = 0:N-1
% calc real and imaginary part of solution
for n = 0:N-1
X(k+1) = X(k+1) + x(n+1)*exp(-j*pi/2*n*k);      % solution by example
end

end

t = 0:N-1
%-----------------------------------
% plots of example solution
subplot(3,2,1)
stem(t,x);
xlabel('Time (s)');
ylabel('Amplitude');
title('1 Time domain - Input sequence')

subplot(3,2,3)
stem(t,abs(X))
xlabel('Frequency');
ylabel('|X(k)|');
title('3 Frequency domain - Magnitude response')

subplot(3,2,5)
stem(t,angle(X))
xlabel('Frequency');
ylabel('Phase');
title('5 Frequency domain - Phase response')

% ------
subplot(3,2,2)
stem(t,x);
xlabel('Time (s)');
ylabel('Amplitude');
title('2 Time domain - Input sequence')

subplot(3,2,4)
stem(t,abs(XFFT))
xlabel('Frequency');
ylabel('|X(k)|');
title('4 Frequency domain - Magnitude response')

subplot(3,2,6)
stem(t,angle(XFFT))
xlabel('Frequency');
ylabel('Phase');
title('6 Frequency domain - Phase response')


One thing in advance, I actually wanted to post this as a comment on Matt's answer. But i couldn't find a way to add images there. Is this possible?

Okay I now updated the example code and also tried my own code again. The calculation of the magnitude now works really fine. But there is still something wrong with the angle. I added my own implementation for angle calculation and if you compare the angles of the example code which are calculated by the build-in angle() function and the angle of my implementation it produces the same results. But if you use the build-in angle() function for the results of fft() there is still a difference.

Example 1 (the left diagrams are of the example code, the diagrams in the middle are from my own code and the right diagrams are from the build-in fft function) In this example it looks for me like there is an error in the sign

Example 2 (the same here) But here it doesn't look like a sign-error anymore

Example 3 (the same here)

does anyone know why there is a difference in the angle?

# example 1
#x = [2 3 -1 4];

# example 2
x = [2 3 -1 4 2 3 -1 4];

# example 3
#Fs = 1000;            % Sampling frequency
#T = 1/Fs;             % Sampling period
#L = 1500;             % Length of signal
#t = ((0:L-1)*T);     % Time vector
#x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);

# solution by build-in fft function
XFFT = (fft(x)');

# solution by two different dft functions
N = length(x);
X = zeros(N,1);
X2R = zeros(N,1);
X2I = zeros(N,1);
XMagnitude = zeros(N,1);
XAngle = zeros(N,1);
for k = 0:N-1
# calc real and imaginary part of solution
for n = 0:N-1
X(k+1) = X(k+1) + x(n+1)*exp(-j*2*pi*n*k/N);      # solution by example
X2R(k+1) = X2R(k+1) + x(n+1)* ( cos(2*pi*n*k/N)); # my solution
X2I(k+1) = X2I(k+1) + x(n+1)* (-sin(2*pi*n*k/N)); # my solution
end

# calc magnitude
XMagnitude(k+1) = sqrt(X2R(k+1)^2+X2I(k+1)^2);

# calc angle
if (X2R(k+1) > 0)
XAngle(k+1) = atan(X2I(k+1)/X2R(k+1));
elseif ( (X2R(k+1) < 0) & (X2I(k+1) >= 0) )
XAngle(k+1) = atan(X2I(k+1)/X2R(k+1)) + pi;
elseif ( (X2R(k+1) < 0) & (X2I(k+1) < 0) )
XAngle(k+1) = atan(X2I(k+1)/X2R(k+1)) - pi;
elseif ( (X2R(k+1) == 0) & (X2I(k+1) > 0) )
XAngle(k+1) = pi/2;
elseif ( (X2R(k+1) == 0) & (X2I(k+1) > 0) )
XAngle(k+1) = -pi/2;
elseif ( (X2R(k+1) == 0) & (X2I(k+1) = 0) )
XAngle(k+1) = 0;
endif

if (mod(k,50) == 0)
fprintf('k %s\n',num2str(k));
fprintf('X %s\n',num2str(X(k+1)));
fprintf('XR %s\n',num2str(X2R(k+1))); # my solution
fprintf('XI %s\n',num2str(X2I(k+1))); # my solution
fprintf('XFFT %s\n',num2str(XFFT(k+1)));
fprintf('Mag %s\n',num2str(abs(X(k+1))));
fprintf('Mag %s\n',num2str(XMagnitude(k+1))); # my solution
fprintf('Mag %s\n',num2str(abs(XFFT(k+1))));
fprintf('Angle %s\n',num2str(angle(X(k+1))));
fprintf('Angle %s\n',num2str(XAngle(k+1)));   # my solution
fprintf('Angle %s\n',num2str(angle(XFFT(k+1))));
fflush(stdout);
endif
end

t = 0:N-1
#-----------------------------------
# plots of example solution
subplot(3,3,1)
stem(t,x);
xlabel('Time (s)');
ylabel('Amplitude');
title('1 Time domain - Input sequence')

subplot(3,3,4)
stem(t,abs(X))
xlabel('Frequency');
ylabel('|X(k)|');
title('4 Frequency domain - Magnitude response')

subplot(3,3,7)
stem(t,angle(X))
xlabel('Frequency');
ylabel('Phase');
title('6 Frequency domain - Phase response')

#-----------------------------------
# plots of my solution
subplot(3,3,2)
stem(t,x);
xlabel('Time (s)');
ylabel('Amplitude');
title('2 Time domain - Input sequence')

subplot(3,3,5)
stem(t,XMagnitude)
xlabel('Frequency');
ylabel('|X(k)|');
title('5 Frequency domain - Magnitude response')

subplot(3,3,8)
stem(t,XAngle)
xlabel('Frequency');
ylabel('Phase');
title('8 Frequency domain - Phase response')

#-----------------------------------
# plots of FFT solution
subplot(3,3,3)
stem(t,x);
xlabel('Time (s)');
ylabel('Amplitude');
title('3 Input sequence')

subplot(3,3,6);
stem(0:N-1,abs(fft(x)));
xlabel('Frequency');
ylabel('|X(k)|');
title('6 Magnitude Response');

subplot(3,3,9);
stem(0:N-1,angle(fft(x)));
xlabel('Frequency');
ylabel('Phase');
title('9 Phase Response');

• Personally, I think that starting from random code found somewhere on the internet is generally a bad idea. My suggestion is to delete that code, and develop your own starting from the DFT equation.
– MBaz
Apr 1, 2018 at 18:14
• The phase is wrong in the N=4 example also. Your question is basically "why is this random code off the Internet wrong?" The answer is: a lot of stuff on the internet is wrong. Be glad you didn't use that function in production code... :) Apr 2, 2018 at 0:42
• Thank's a lot for your replies. Of course you are right, that the other way around would have made more sense. I actually even started by developing my own code. But then I noticed that I have not understood the topic enough. And so I decided to search for a working solution and debug this solution until i understand what happens here. Apr 2, 2018 at 9:05
• Converted your "Answer" to an edit. I suggest asking the question in a new question on this site, given you've already awarded the check mark to Matt.
– Peter K.
Apr 2, 2018 at 16:31

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}\tag{1}$$
The code does not correctly implement Eq. $(1)$. The argument of the exponential function should be -j*2*pi*n*k/N, where N is the DFT length. For N=4 (as in ex. 1), the code happens to be correct. However, for any other DFT length it will fail (as it does for ex. 2 where we have N=8).
• Hi @Matt L. thank's a lot for your fast reply!! I actually even wondered about that the term -j*pi/2*n*k differed from other solutions i found in the internet, but unfortunately i did not question it enough. But anyway, thank's a lot for your help. Now it works fine!! Apr 2, 2018 at 9:03