% FFT 1D

% f(x)=delta(x) 
f=zeros(1,8);
f(1)=1
F=fft(f)
nf=ifft(F)


% f(x)=1
f=ones(1,8)
F=fft(f)
nf=ifft(F)


%f(x)=cos(x)
N=128;
k0=5;
x=0:N-1;
x=2*pi.*x*k0/N;
f=cos(x);
plot(f);
F=fft(f);
plot(abs(F));
SF=fftshift(F);
plot(abs(SF));
plot(real(SF));
plot(imag(SF)); axis([0 128 -1 1]);
plot(angle(SF));
plot(angle(SF).*abs(SF)>1);


%f(x)=sin(x)
N=128;
k0=5;
x=0:N-1;
x=2*pi.*x*k0/N;
f=sin(x);
plot(f);
F=fft(f);
plot(abs(F));
SF=fftshift(F);
plot(abs(SF));
plot(real(SF)); axis([0 128 -1 1]);
plot(imag(SF)); 
plot(angle(SF).*(abs(SF)>0.5));

% Gaussians with various std's
N=128;
for i=1:15,
   x=-64:63;
   f=zeros(1,N);
   f=exp(-(x.^2)./(i^2));
   subplot(2,1,1); plot(f); 
   F=fft(f);
	SF=fftshift(F);
   subplot(2,1,2); plot(abs(SF)); 
   pause;
end;



%f(x)=rect(x)
f=zeros(1,N);
M=10;
f(1:M)=1;  f(N:-1:N-M)=1;
subplot(2,1,1); plot(fftshift(f));
axis([0 140 0 2]);

F=fft(f);
SF=fftshift(F);
subplot(2,1,2); plot(abs(SF));


%Building rect wave
f=ones(1,32);
f(1:5)=-1; f(11:15)=-1; f(21:25)=-1;
figure;
subplot(3,1,1); plot(f); axis ([1 36 -2 2]);

F=fft(f);
F=fftshift(F);
T=zeros(1,32);
T(17)=F(17);
subplot(3,1,3); plot(real(ifft(fftshift(T)))); axis ([1 36 -2 2]);
for i=1:15,
   A=zeros(1,32);
   A(17+i)=F(17+i);
   A(17-i)=F(17-i);
   T=T+A;
   subplot(3,1,3); plot(real(ifft(fftshift(T)))); axis ([1 36 -2 2]);
   subplot(3,1,2); plot(real(ifft(fftshift(A)))); 
   pause;
end;
A=zeros(1,32);
A(1)=F(1);
T=T+A;
subplot(3,1,3); plot(real(ifft(fftshift(T)))); axis ([1 36 -2 2]);
subplot(3,1,2); plot(real(ifft(fftshift(A)))); 
% figure;
% plot(abs(F));
   



% the symmetric property of real functions f(x)=rect(x)
close;
f=zeros(1,N);
M=10;
f(1:M)=1;  f(N:-1:N-M)=1;
subplot(2,1,1); plot(fftshift(f));
axis([0 140 0 1.5]);
F=fft(f);
SF=fftshift(F);
subplot(2,1,2); plot(abs(SF));

figure;
subplot(2,1,1); plot(1:floor(N/2)+1,real(SF(1:floor(N/2)+1)),'r');
hold on; plot(floor(N/2)+1:N,real(SF(floor(N/2)+1:end)),'b');
title('real part');
subplot(2,1,2); plot(1:floor(N/2)+1,imag(SF(1:floor(N/2)+1)),'r');
hold on; plot(floor(N/2)+1:N,imag(SF(floor(N/2)+1:end)),'b');
title('imaginary part');



%%%% 2D


% 2d rect:  f(x,y)=rect(x,y)
close;
N=128;
[x,y]=meshgrid(-64:63);
im=abs(x)<10 & abs(y)<10;
figure; imagesc(im); colormap(gray);

IM=fft2(im);
imagesc(abs(IM));
imagesc(log(abs(IM)+1));
imagesc(fftshift(log(abs(IM)+1)));

% fft of Lena
close;
im=readImage('stroller.tif');
putImage(im);
IM=fft2(im);
imagesc(log(abs(fftshift(IM))+1));


%Separability of fft
close all;
N=128;
[x,y]=meshgrid(-64:63);
im=abs(y)<10 ;
imagesc(im); colormap(gray);

Im=fft(im); %apply fft in the y direction
figure;
imagesc(fftshift(log(abs(Im)+1)));colormap(gray);

IM=fft(Im')';  %apply fft in the x direction
imagesc(fftshift(log(abs(IM)+1)));

% Reverse order
Im=fft(im')'; %apply fft in the x direction
figure;
imagesc(fftshift(log(abs(Im)+1)));colormap(gray);

IM=fft(Im);  %apply fft in the y direction
imagesc(fftshift(log(abs(IM)+1)));

% plot(fftshift(abs(IM(:,1))));

%convolution theorem (1D)
a=zeros(64,1);
a(1:10)=1;
subplot(3,1,1); plot(a); title('signal f(x)'); 
axis([0 70 0 1.5]);

b=conv(a,a);
subplot(3,1,2); plot(b(1:64)); title('signal f(x)*f(x)'); 

A=fft(a);
B=A.*A;
b=real(ifft(B));
subplot(3,1,3); plot(b); title('ifft(F(x).F(x))'); 


%%%%%%%%%%%%%%%% Filters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;


% image scaling by optimal interpolation

close all
img = readImage('stroller.tif');
showImage(img);
fimg = (fft2(img)); %the real fft shifted for symmetry
fimgRep = repmat(fimg,2,2);
% show the image of the replicated transform
nimg = ifft2(fimgRep);
nimg = round(real(nimg));
showImage(nimg);
% create rect
rect = zeros(size(fimgRep));
rect(size(fimgRep,1)./4+1 : 3.*size(fimgRep,1)./4,size(fimgRep,2)./4+1 : 3.*size(fimgRep,2)./4) = 1;
showImage(rect.*255);
%interpolate by mult FFT with rect
fimg_interp = fimgRep.*rect;
img_inter = ifft2(fftshift(fimg_interp));
nimg = round(real(img_inter));
showImage(nimg);
image(4.*nimg);



%filtering
close;

img = readImage('stroller.tif');
ffs = fftshift(fft2(img)); %the real fft shifted for symmetry
fimg = log(abs(ffs)+1); %the log of fft to be able to see it in good visible proportions

%this is a meshgrid to use with the calculation of the round filter
[x y] = meshgrid(-1:(2/size(img,1)):1-(2/size(img,1)), -1:(2/size(img,2)):1-(2/size(img,2)));
circledist = sqrt(x.^2+y.^2);

nimg_fig = figure;
image(img); axis off; axis image; colormap(gray(256));

ffilter_fig = figure;
image([-64:63],[-64:63],zeros(size(fimg))); colormap(gray); axis image;

fimg_fig = figure;
imagesc([-64:63],[-64:63],fimg); colormap(gray); axis image;

pause;

%Ideal Low pass
for i=0.5:-0.05:0.05,
  filter = circledist<i;
  figure(ffilter_fig);
  imagesc(filter); colormap(gray); axis image;
  figure(fimg_fig);
  imagesc(abs(fimg.*filter));  colormap(gray);   axis image;
  figure(nimg_fig);
  imagesc(real(ifft2(fftshift(ffs.*filter))));  colormap(gray); axis image;
  pause;
  %imagesc([-64:63],[-64:63],ffl); colormap(gray); 
end;


%Gaussian low pass
for i=0.5:-0.05:0.05,
  filter = exp(-circledist.^2/(i^2));
  figure(ffilter_fig);
  imagesc(filter); colormap(gray); axis image;
  figure(fimg_fig);
  imagesc(abs(fimg.*filter));  colormap(gray);   axis image;
  figure(nimg_fig);
  imagesc(real(ifft2(fftshift(ffs.*filter))));  colormap(gray); axis image;
  pause;
  %imagesc([-64:63],[-64:63],ffl); colormap(gray); 
end;

%Ideal high pass
for i=0.1:0.1:1,
  filter = circledist>(i) ; 
  figure(ffilter_fig);
  imagesc(filter); colormap(gray);  axis image;
  figure(fimg_fig);
  imagesc(abs(fimg.*filter));  colormap(gray);   axis image;
  figure(nimg_fig);
  imagesc(real(ifft2(fftshift(ffs.*filter))));  colormap(gray);  axis image;
  pause;
  %imagesc([-64:63],[-64:63],ffl); colormap(gray); 
end;

%band pass;
for i=0.0:0.1:0.9,  
  filter = circledist<(i+0.1) & circledist >(i) ; 
  figure(ffilter_fig);
  imagesc(filter); colormap(gray);  axis image;
  figure(fimg_fig);
  imagesc(abs(fimg.*filter));  colormap(gray);   axis image;
  figure(nimg_fig);
  imagesc(real(ifft2(fftshift(ffs.*filter))));  colormap(gray);  axis image;
  pause;
  %imagesc([-64:63],[-64:63],ffl); colormap(gray); 
end;



