MATLAB Calculating Discrete Fourier Transform Coefficient with Matlab

[drdizzard]

I've been asked to write a function (.m file) in Matlab to calculate the discrete Fourier transform coefficient for an arbitrary function x. So far this is what I've done:

function a = mydft(x,N)
%MYDFT Calculates the discrete Fourier transform
%usage:
%[a]=mydft(x)
%x=[ x[0] x[1] ... x[N-1] ] - vector containing discrete time sequence
%a=[ a_0 a_1 ... a_N-1 ] - vector of discrete time Fourier Series
%coefficients, k=0,...,N-1
for k=0:N-1
    for m=0:N-1
    w=exp(-j*2*pi/N.*m.*k);
    end
    a=sum(x(m).*w);
end
end

when I use it on a signal given to me and then plot kw_hat vs. |a| I get a plot that looks nothing like what I do for the same signal by hand. I don't think I'm getting the correct DFT coefficients out of this function

This is the first time I've used Matlab so I'm not very familiar with it and would appreciate any advice.

 

[osilmag]

I think you might have an index error. You should allocate vector space. So, for "a", you should write a(1), a(2), a(3), etc. Otherwise, your program might overwrite previous values that need to be there.

 

[olivermsun]

The code seems to mix and match the traditional C/FORTRAN programming style (inner and outer loops with accumulation) with Matlab (vector products), so you may be confusing yourself there.If you want to do it the traditional way then you want to calculate the product a = x(m)*w at every step of the loop and accumulate. You can do this by putting something like a(k) = a(k) + x(m)*w inside the innermost loop.

You can also compute the basis function w(k) = exp(...) once for each k and then integrate the product by
using point wise multiplication (.*) and summing. This allows you to eliminate the innermost loop completely. Just set up a vector

m = 0:N-1;   % only need to define this once

then the loop can look like

for k = 0:N-1
    w = exp(-j*2*pi/N.*m*k);
    a(k+1) = sum(x.*w);  % watch the index: a(1) is the first array element
end

(note that you have to index on k as osilmag pointed out so that you retain the individual a's for each k)

Finally, you can eliminate both loops if you pass to exp() an NxN array which varies in m in one dimension and k in the other...