% dft1.m

dat = imread('test_image.tif');
f = double(dat(50,:,1)).'; % example data: 50th row of red channel

N=length(dat);
for k=1:N,
for j=1:N,
   A(k,j) = exp(-2i*pi*(k-1)*(j-1)/N);
end
end

F = A*f;

subplot(3,1,1);
plot(f,'k-'); xlabel('j'); ylabel('Intensity'); title('original data');
subplot(3,1,2);
plot(real(F),'k-'); axis([0 800 -1e4 1e4]); hold on;
plot(imag(F),'r-'); xlabel('k'); ylabel('F'); title('Fourier transformed data');


max(inv(A)*F - f) % indeed, it is the inverse...
max(F - fft(f))   % indeed, the conventions equal MATLAB's

F(100:700)=0;
subplot(3,1,3);
fnew = real(inv(A)*F);
plot(fnew,'k-'); xlabel('j'); ylabel('Intensity'); title('reconstructed from 200 low frequency modes');