different solutions in matlab / octave using dft and fft

Question

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'); 

Answer 1

The solution is simple, and it would have been sufficient to check the code against the DFT formula:

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