clc,clear,close all % 清理命令区、清理工作区、关闭显示图形
warning off % 消除警告
feature jit off % 加速代码运行
D0 = 20; % 阻止的频率点与频域中心的距离
n = 2; % 阶次
im = imread(‘coloredChips.png’); % 原图像
R = imnoise(im(:,:,1),‘gaussian’,0,0.01); % R + 白噪声
G = imnoise(im(:,:,2),‘gaussian’,0,0.01); % G + 白噪声
B = imnoise(im(:,:,3),‘gaussian’,0,0.01); % B + 白噪声
im = cat(3,R,G,B); % 原图像 + 白噪声
R1 = freqfilter_btw_lp(R,D0,n); % 巴特沃斯低通滤波器
G1 = freqfilter_btw_lp(G,D0,n); % 巴特沃斯低通滤波器
B1 = freqfilter_btw_lp(B,D0,n); % 巴特沃斯低通滤波器
im1 = cat(3,R1,G1,B1);
figure(‘color’,[1,1,1])
subplot(121),imshow(im,[]); title(‘原始图像’)
subplot(122),imshow(im1,[]); title(‘巴特沃斯低通滤波图像’);
function im5 = freqfilter_btw_lp(im,D0,n)
if ~isa(im,‘double’)
im1 = double(im)/255;
end
im2 = fft2(im1); % 傅里叶变换
im3 = fftshift(im2); % 中心化
[N1, N2] = size(im3);
n1 = fix(N1 / 2);
n2 = fix(N2 / 2);
for i = 1:N1
for j = 2:N2
d = sqrt((i-n1)2+(j-n2)2);
h = 1/(1 + 0.414 * (d / D0)^(2*n)); % 巴特沃斯低通滤波器
result(i,j) = h * im3(i,j);
end
end
result = ifftshift(result); % 反中心化
im4 = ifft2(result); % 反变换
im5 = im2uint8(real(im4)); % 滤波图像
end