function imt = fftcenter(ima) % 图像中心化转换函数

[xs,ys] = size(ima);

for i = 1:xs

for j = 1:ys

imt(i,j) = ima(i,j)*power((-1),(i+j)) ;%输入图像每个像素点乘以(-1)^x+y因子

end

end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 2-D FFT 程序包

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function proj04();

ima = imread('Fig4.41(a).jpg');

imshow(ima);

title('原始图像');

ima_f = fft2(ima); % 对图像ima求二维快速傅里叶变换

mod_ima = abs(ima_f) % 要求幅度谱，需对傅里叶变换求模

mod_ima_log = uint8(log2((mod_ima))) % 要求幅度谱需对FFT后得到的模取对数

phase_ima = atan2(imag(ima_f),real(ima_f)); % 求相位谱，arctan必须使用一个四象限反正切来计算，调用MATLAB的atan2函数

figure,subplot(1,2,1),imshow(mod_ima_log,[]); % []的作用是将mod_ima_log的最大最小值分别作为纯白255和纯黑0，中间值映射到0到255之间

title('输入图像的二维傅里叶幅度谱');

subplot(1,2,2),imshow(phase_ima,[]);

title('输入图像的二维傅里叶相位谱');

threeD_mod_ima_log = log2(mod_ima)+1;

[x1,y1] = size(threeD_mod_ima_log);

[a1,b1]=meshgrid(1:x1,1:y1); % 使用函数mesh绘制三维傅立叶图需要先进行meshgrid

figure,subplot(1,2,1),mesh(a1,b1,threeD_mod_ima_log);

shading interp; % 在网格片内采用颜色插值处理，使得三维表面图显得光滑

title('输入图像的三维傅里叶幅度谱');

[x2,y2] = size(phase_ima);

[a2,b2] = meshgrid(1:x2,1:y2); % 使用函数mesh绘制三维傅立叶图需要先进行meshgrid

subplot(1,2,2),mesh(a2,b2,phase_ima);

shading interp;

title('输入图像的三维傅里叶相位谱');

%

s = fftshift(ima_f); % 调用fftcenter函数实现图像的中心化变换(为了对比fftcenter)

mod_s = abs(s) % 要求幅度谱，需对傅里叶变换求模

mod_s_log = uint8(log2(mod_s)) % 要求幅度谱需对FFT后得到的模取对数

phase_s = atan2(imag(s),real(s)); % 求相位谱，arctan必须使用一个四象限反正切来计算，调用MATLAB的atan2函数

figure,imshow(mod_s_log,[]); % []的作用是将mod_ima_log的最大最小值分别作为纯白255和纯黑0，中间值映射到0到255之间

title('fftshift函数的中心变换图的二维傅里叶幅度谱');

imt = fftcenter(ima); % 调用自编函数对图像进行中心化变换

imt_f = fft2(imt); % 对中心化变换后的图像求二维傅里叶变换

figure,imshow(imt_f,[]);

title('使用自编函数fftcenter()中心变换图的二维傅里叶图谱');

mod_imt = abs(imt_f); % 要求幅度谱，需对傅里叶变换求模

mod_imt_log = uint8(log2(mod_imt)); % 要求幅度谱需对FFT后得到的模取对数

phase_imt = atan2(imag(imt_f),real(imt_f)); % 求相位谱，arctan必须使用一个四象限反正切来计算，调用MATLAB的atan2函数

figure,subplot(1,2,1),imshow(mod_imt_log,[ ]);

title('中心变换图的二维傅里叶幅度谱');

subplot(1,2,2),imshow(phase_imt,[]),

title('中心变换图的二维傅里叶相位谱');

%

threeD_mod_imt_log = log10(mod_imt)+1;

[x3,y3] = size(threeD_mod_imt_log);

[a3,b3]=meshgrid(1:x3,1:y3); % 使用函数mesh绘制三维傅立叶图需要先进行meshgrid

figure,subplot(1,2,1),mesh(a3,b3,threeD_mod_imt_log);

shading interp;

%figure,plot3((1:x),(1:y),ak); shading interp;

title('中心变换图的三维傅里叶幅度谱');

[x4,y4] = size(phase_imt);

[a4,b4] = meshgrid(1:x4,1:y4); % 使用函数mesh绘制三维傅立叶图需要先进行meshgrid

subplot(1,2,2),mesh(a4,b4,phase_imt);

shading interp;

title('中心变换图的三维傅里叶相位谱');

%

imt_if = ifft2(imt_f);

imt_if_abs = abs(imt_if);

figure,subplot(1,2,1),imshow(imt_if_abs,[]);

title('中心化变化图像');

ima = imread('Fig4.41(a).jpg');

subplot(1,2,2),imshow(ima);

title('原始图像');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 设计一个高斯滤波器

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function im=Lowpassfiltering()

st = imread('Fig4.41(a).jpg');

[M,N] = size(st);

D0 = 20; % 设定截止频率D0

x = 2*M;

y = 2*N;

% 进行高斯滤波之前先对图像进行填充，填充后的图像尺寸为(2M*2N)

for i=1:x

for j=1:y

if((i<=M)&&(j<=N))

st(i,j)=st(i,j);

else

st(i,j)=0;

end

end

end

figure,imshow(st);

title('对输入图像进行填充');

imwrite(st,'K1.jpg');

figure(2);

st = fftshift(fft2(st)); % 自编函数中心化效果有问题，因此调用fftshift函数实现中心化

for i = 1:x

for j = 1:y

D(i,j) = sqrt((i-M)^2+(j-N)^2); % 中心点为(M,N)

h(i,j) = 1*exp(-1/2*(D(i,j)^2/D0^2)); % 系统函数H(u,v)

st(i,j) = h(i,j)*st(i,j); % 滤波过程

end

end

st = fftshift(ifft2(st)); % 进行傅里叶反变换时再进行一次中心化

st_abs = abs(st);

for i=1:M

for j=1:N

st_k(i,j)=st_abs(i,j); % 对图像进行裁剪，去除填充部分

end

end

imshow(st_k,[])

title('GLPF滤波后的输出图像(已裁剪)');

st_k1 = imread('Fig4.41(a).jpg');

for i=1:M

for j=1:N

st_k2(i,j) = st_k1(i,j) - ( st_k(i,j)); % 原始图像减去GLPF图像，得到锐化图像

end

end

figure,imshow(st_k2,[])

title('原始图像减去GLPF滤波图像后的锐化图像');