function [df  M omegas W ff] = RealDFT(f)

N=length(f);
n=N/2;

if n==ceil(n)
    
 % M is the matrix whose columns span the real record space. 
 % We list those corresponding to cosine values first, then sines.
 % When you start with an even number of data points (as we do here) 
 % there will be two more cosine coefficients than sine coefficients.

for m=1:(n+1)

for k=1:N

M(k,m)=cos(2*pi*(m-1)*(k-1)/N);

end

end

for m=2:n

for k=1:N

M(k,m+n)=sin(2*pi*(m-1)*(k-1)/N);

end

end

df=M\f;

% omegas is the basis vector whose ``powers'' span the complex record space

for k=1:N
omegas(k,1)=exp(2*pi*i*(k-1)/N);
end

% W has columns that are a basis for the complex record space

for k=1:N
W(:,k)=omegas.^(k-1);
end

% ff is just the discrete Fourier transform calculated with MATLAB's built-in algorithm

ff=fft(f);

else 

display('the length of the input vector must be even')

end