%二维DFT
clear
Filename='67079.bmp'
P=imread(Filename);
imshow(P)
P=double(P)/255;
PP=Colortogray(P);
figure
subplot(2, 3, 1)
imshow(PP)
%计算PP的DFT
[M,N]=size(PP);
FPP=FFT2(PP);
subplot(2, 3, 2)
imshow(abs(FPP)/max(max(abs(FPP))))
FPPs=FFTSHIFT(FPP);
subplot(2, 3, 3)
imshow(abs(FPPs)/max(max(abs(FPPs))))
IFPP=IFFT2(FPP);
subplot(2, 3, 4)
imshow(IFPP)
delta=1;
L=4;
w0=pi;
for gama=3:4
LL=floor(6.6*2^gama*delta+1);
m0=(2^gama*delta)^2/2;
MM=M+LL-1;
NN=N+LL-1;
PPP(MM,NN)=0;
for i=1:M
for j=1:N
PPP(i,j)=PP(i,j);
end
end
subplot(2, 3, 5)
imshow(PPP)
FPPP=FFT2(PPP);
for l=2:L-1
sita=l*pi/L;
for u=1:MM
for v=1:NN
FPPPP(u,v)=FPPP(u,v)*exp(-m0*((2*pi*u/MM-2^(-gama)*w0*cos(sita))^2+(2*pi*v/NN-2^(-gama)*w0*sin(sita))^2));
end
end
IFPPPP=IFFT2(FPPPP);
figure
imshow(real(IFPPPP))
end
end