% Generate sinusoidal component image and plot the spectrum calculated and
% displayed in a separable manner.
% need to be in directory c:\jl\4F8_ImageProcessing\ImageProcessing1\ImageProc1Mat\ImageProc1Demo

% sinusoid_demo

clear,clc,close all,home

disp('Use following values in demo. N=64, f1=1/16, f2=1/32');
disp('Note that number of cycles =f*N');
N=input('Number of samples in horizontal and vertical direction, N = ');
freq1=input('Normalised horizontal frequency(<0.5) = ');
freq2=input('Normalised vertical frequency(<0.5) = ');
n=(0:N-1);
g=zeros(N,N);

vec1=sin(2*pi*freq1*n);

vec2=sin(2*pi*freq2*n);
g=vec2'*vec1;

figure(1);
mesh(g)
title(['2-D Sinusoid   f1=',num2str(freq1),'  f2=',num2str(freq2)]);
  Az=[-70,-30,-5];
  view(Az); rotate3d on;
  pause
%end
% Compute Spectrum in u1 direction
Gu1=zeros(N,N);
for i=1:N
	Gu1(i,:)=fftshift(fft(g(i,:))); 
%       note, fftshift shifts the DC component to the centre of the 
%       spectrum; useful for visualisation
end

figure(2);
f=(n-N/2)/N;
mesh(f,(0:N-1),abs(Gu1)); rotate3d on;
title(['Spectrum in x direction  f1=',num2str(freq1),'  f2=',num2str(freq2)]);
pause


figure(3)
% now compute spectrum in u2 direction
G1 = zeros(N,N);
for j=1:N
	G1(:,j)=fftshift(fft(Gu1(:,j))); 
%       note, fftshift shifts the DC component to the centre of the 
%       spectrum; useful for visualisation
end
mesh(f,f,abs(G1)); rotate3d on;
title(['2-D Spectrum   f1=',num2str(freq1),'  f2=',num2str(freq2)]);
pause


% note that we get precisely that same end result if we use the Matlab command fft2

figure(4)
G=fftshift(fft2(g));
%G = fft2(fftshift(g));
mesh(f,f,abs(G)); rotate3d on;
title(['2-D Spectrum   f1=',num2str(freq1),'  f2=',num2str(freq2)]);
pause

% see what happens if we do not use fftshift
figure(5)
G2=(fft2(g));
mesh(f,f,abs(G2)); rotate3d on;
title(['2-D Spectrum   f1=',num2str(freq1),'  f2=',num2str(freq2)]);
% should really rescale our axes in this case.

% look at this with 1d case

figure(6);
plot(vec1);
fvec=zeros(N);
fvec=fft(vec1);
figure(7);
plot(abs(fvec));
figure(8);
plot(abs(fftshift(fvec)));



break

