clc
%~~~读入图像~~~
path = input('D:\MATLAB7\work，please：','s');
imrgb = imread(path);
imrgb = imread('4.jpg');
imshow(imrgb)
title('原始图像');
[x,y,z] = size(imrgb); %rgb图像是三维的！若写成[x1,y1] = size(im)，则y1 = y*z~
%~~~转为灰度图像~~~
imgray = rgb2gray(imrgb);
figure,imshow(imgray)
title('转换成灰度图像');
[x2,y2,z2] = size(imgray); %超维时赋1，此时z2 = 1~
%~~~double处理~~~
imgray_dble = double(imgray); % 须转为double类型~
%~~~~低频居中~~~
for i = 1:x2
for j = 1:y2
imgray_dble(i,j) = imgray_dble(i,j)*(-1)^(i+j);
end
end
%~~~二维快速傅里叶变换~~~
im_fft2 = fft2(imgray_dble);
%~~~滤波前，log处理显示~~~
im_fft2_log =log(1+im_fft2);
im_fft2_log = histeq( uint8(real(im_fft2_log)) );
figure,imshow( im_fft2_log ),title('低通滤波前频域图像');
%~~~低通滤波，可以改变滤波半径~~~
for i=1:x2
for j=1:y2
if( ((i-x2/2)^2 + (j-y2/2)^2 )>10000 )
im_fft2(i,j)=0;
end
end
end
%~~~滤波后，log处理显示~~~
im_fft2_log = log(1+im_fft2);
im_fft2_log = histeq( uint8(real(im_fft2_log)) );
figure,imshow( im_fft2_log ),title('低通滤波后频域图像');
%~~~fouier反变换~~~
im_ifft2 = ifft2(im_fft2);
%取实部
im_ifft2_real=real(im_ifft2);
%~~~反低频居中~~~
for i = 1:x2
for j = 1:y2
im_ifft2_real(i,j) = im_ifft2_real(i,j)*(-1)^(i+j);
end
end
%~~~显示处理结果~~~
figure,imshow( uint8(im_ifft2_real) ),title('滤波后对应图像');
