Unexpected result with DFT in MATLAB

Question

I have a problem when calculate discrete Fourier transform in MATLAB, apparently get the right result but when plot the amplitude of the frequencies obtained you can see values very close to zero which should be exactly zero. I use my own implementation:

function [y] = Discrete_Fourier_Transform(x)
    N=length(x);
    y=zeros(1,N);
    for k = 1:N
        for n = 1:N
            y(k) = y(k) + x(n)*exp( -1j*2*pi*(n-1)*(k-1)/N );
        end;
    end;
end

I know it's better to use fft of MATLAB, but I need to use my own implementation as it is for college.

The code I used to generate the square wave:

x = [ones(1,8), -ones(1,8)];
for i=1:63
    x = [x, ones(1,8), -ones(1,8)];
end

MATLAB version: R2013a(8.1.0.604) 64 bits

I have tried everything that has happened to me but I do not have much experience using MATLAB and I have not found information relevant to this issue in forums. I hope someone can help me.

Thanks in advance.

I assume this is a numerical problem. As these values are in the range of 1e-15 while the peaks have an amplitude of about 656 I wouldn't bother much. The sum of the squared error between the MATLAB FFT and your DFT routine is in the range of 1e-20, so basically zero. But maybe somebody else has a more detailed explenation? — hbaderts
I agree with @hbaderts. The values are very small, almost as small as eps, and can be considered 0 for all purposes here. Your FFT solution looks great. — siliconwafer
Thanks for your answers, i know this values are negligible and do not show in the plot but when i want to plot the phase with angle function obtain meaningless results under. I am using y((y < 0.00001) & (y > -0.00001)) = 0 to hide this values but it is not an elegant solution. — Adrián Álvarez Sáez

hbaderts hbaderts · Accepted Answer · 2015-04-30T14:39:54

This will be a numerical problem. The values are in the range of 1e-15, while the DFT of your signal has values in the range of 1e+02. Most likely this won't lead to any errors when doing further processing. You can calculate the total squared error between your DFT and the MATLAB fft function by

y = fft(x);
yh = Discrete_Fourier_Transform(x);
sum(abs(yh - y).^2)
ans = 
    3.1327e-20

which is basically zero. I would therefore conclude: your DFT function works just fine.

Just one small remark: You can easily vectorize the DFT.

n = 0:1:N-1;
k = 0:1:N-1;
y = exp(-1j*2*pi/N * n'*k) * x(:);

With n'*k you create a matrix with all combinations of n and k. You then take the exp(...) of each of those matrix elements. With x(:) you make sure x is a column vector, so you can do the matrix multiplication (...)*x which automatically sums over all k's. Actually, I just notice, this is exactly the well-known matrix form of the DFT.

Unexpected result with DFT in MATLAB

1 Answers