0
$\begingroup$

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

enter image description here

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.

enter image description here

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.

enter image description here

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 1

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

Example 3 (the same here) example 3

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'); 
$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – MBaz Apr 1 '18 at 18:14
  • $\begingroup$ 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... :) $\endgroup$ – Cris Luengo Apr 2 '18 at 0:42
  • $\begingroup$ 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. $\endgroup$ – anon1234 Apr 2 '18 at 9:05
  • $\begingroup$ 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. $\endgroup$ – Peter K. Apr 2 '18 at 16:31
1
$\begingroup$

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

$\endgroup$
  • $\begingroup$ 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!! $\endgroup$ – anon1234 Apr 2 '18 at 9:03

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.