Calculating own dft via matlab?

Question

We are asked to code our own dft function from the formula :

If everything is done correctly it should give the same result with matlab's own dft function, in the end I'm comparing them but they are not the same... We are not allowed to use for loops this is why Im trying to do it with matrices. Below is my code:

x=linspace(0,pi,10);
y=cos(x);

ones_vec = ones(10,1);
ones_vec=(ones_vec*y);

k=linspace(1,10,10);
n=linspace(0,9,10);
ones_vec_n = ones(10,1);
ones_vec_n = (ones_vec_n*n).';

ones_vec2 = ones(10,1);
ones_vec2=(ones_vec2*k*-1i*2*pi/10);
%n = n*ones_vec2.'
ones_vec2= ones_vec2 .* ones_vec_n;

dft_my = sum(ones_vec2,1)

fft_matlab = fft(y)

I would be glad if you can give a hand!

Answer 1

I checked the following in Octave:

N = 128;             % FFT length
M = exp(-1i*(0:N-1)'*(0:N-1)*2*pi/N);    % DFT matrix
x = randn(128,1);    % input vector
X2 = M*x;            % DFT
X = fft(x);          % built-in FFT
error = max(abs(X-X2));    % max. magnitude of error (= 5.6474e-13)
                           % (due to numerical errors)

Answer 2

Let's see if this helps:

dk = 0.01;
N  = 10;

[k,n] = meshgrid(0:dk:10, 0:1:N-1);

x  = cos(n);
X  = sum(x.*exp(-1i*k.*n*2*pi/N));
plot(k,abs(X));

I'm not very familiar with discrete Fourier's transform and I'm not sure if $X(n)$, which is a complex value, is what really should be.

I hope this helps.

Cheers!

Answer 3

I don't have access to matlab at the moment, so I can't double check this, but it's really easy to construct the DFT matrix. All you have to do is take the outer product of 0:N-1 with itself in the phase of the complex exponential:

N = 256;
D = exp(1j * 2*pi/N * (0:N-1)' * (0:N-1));

Then the DFT of a length N vector x is then just the matrix-vector product:

D*x

To check, use something like this:

maxerr = max(abs(fft(x) - D*x))