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