# FFT a aplikace Fourierovy analyzy
function cviceni9

graphics_toolkit("gnuplot")

#DFT
L=10;
N=64;
delta=L/N;
x=(-N/2:N/2-1)*delta;
x0=5*delta;
f=(-N/2:N/2-1)/(N*delta);
g=exp(-pi*(x-x0).^2); #vzorkovana fce
G=ifftshift(fft(fftshift(g)))*delta; #DFT

#spojite verze
xs=linspace(x(1),x(end),500);
fs=linspace(f(1),f(end),500);
gs=exp(-pi*(xs-x0).^2);
Gs=exp(-i*2*pi*x0*fs).*exp(-pi*fs.^2);

#vzorkovani
figure(1)
plot(x,g,"x",xs,gs)
#velikost FT
figure(1)
plot(f,abs(G),"x",fs,abs(Gs))
#faze FT
figure(2)
plot(f,arg(G),"x",fs,arg(Gs))


#sinc funkce
N=64
L=32
delta=L/N;
x=(-N/2:N/2-1)*delta;
g=sinc(x);
G=ifftshift(fft(fftshift(g)));
figure(3)
plot(abs(G))


#konvoluce

N=128;
g(1:N)=0;
g(N/2-30:N/2+30)=1;
figure(1);
plot(g,"o")
axis('tight');
x=linspace(-5,5,N);
h=exp(-2*x.*x);
figure(2);
plot(h,"o");
axis('tight');
G=fft(g);
H=fft(h);
konv=ifftshift(ifft(G.*H));
figure(3);
plot(real(konv),"o");
axis('tight');


# konvoluce + nastaveni

L=6
N=256
del=L/N;
x=(-N/2:N/2-1)*del;
delf=1/del/N;
fx=(-N/2:N/2-1)*delf;

g=4+floor(x);
h=exp(-20*x.^2); 
h/=sum(h)*del;
sum(h)*del

figure(1)
plot(x,g,"linewidth",2)

figure(2)
plot(x,h,"linewidth",2)

G=del*fft(fftshift(g));
H=del*fft(fftshift(h));
S=G.*H;
s=real(ifftshift(ifft(S)))/del;

figure(3)
plot(x,s,"linewidth",2)


# nastaveni nulami

Npad=100; pad=zeros(1,Npad/2);
xx=x=(-N/2-Npad/2:N/2-1+Npad/2)*del;
gg=[pad g pad];
hh=[pad h pad];

figure(4)
plot(xx,gg,"linewidth",2)

figure(5)
plot(xx,hh,"linewidth",2)

GG=del*fft(fftshift(gg));
HH=del*fft(fftshift(hh));
SS=GG.*HH;
ss=real(ifftshift(ifft(SS)))/del;

figure(6)
plot(xx,ss,"linewidth",2)




# obraz

close all

g=zeros(N);
g(N/2-30:N/2+30,N/2)=1;
g(N/2,N/2-30:N/2+30)=1;
figure(1)
imshow(g)
[x,y]=meshgrid(linspace(-10,10,N));
h=exp(-x.^2-y.^2-x.*y);
figure(2)
imshow(h)
G=fft2(g);
H=fft2(h);
rozm=ifftshift(ifft2(G.*H));
rozm=rozm./max(max(rozm));
figure(3)
imshow(rozm)


endfunction
