【修订版】本文已整理到我的个人博客,订正了原文的几处错误:巴特沃斯带阻滤波器的阶数与代码不一致(代码是 3 阶),代码节的标题有误,陷波滤波器连乘式的符号与图像尺寸的记号冲突。CSDN 这一版保留原样,之后的勘误与更新只在博客进行:
Frequency Domain Filter : 高通滤波器 | Finn's blogFrequency Domain Filter : 带阻滤波器与陷波滤波器 | Finn's blog
1.高通滤波器
1.1理想高通滤波器
1.2巴特沃斯高通滤波器
1.3高斯高通滤波器
1.4 上述三种滤波器的实现Matlab代码
close all; clear all; %% ---------Ideal Highpass Filters (Fre. Domain)------------ f = imread('characters_test_pattern.tif'); f = mat2gray(f,[0 255]); [M,N] = size(f); P = 2*M; Q = 2*N; fc = zeros(M,N); for x = 1:1:M for y = 1:1:N fc(x,y) = f(x,y) * (-1)^(x+y); end end F = fft2(fc,P,Q); H_1 = ones(P,Q); H_2 = ones(P,Q); for x = (-P/2):1:(P/2)-1 for y = (-Q/2):1:(Q/2)-1 D = (x^2 + y^2)^(0.5); if(D <= 60) H_1(x+(P/2)+1,y+(Q/2)+1) = 0; end if(D <= 160) H_2(x+(P/2)+1,y+(Q/2)+1) = 0; end end end G_1 = H_1 .* F; G_2 = H_2 .* F; g_1 = real(ifft2(G_1)); g_1 = g_1(1:1:M,1:1:N); g_2 = real(ifft2(G_2)); g_2 = g_2(1:1:M,1:1:N); for x = 1:1:M for y = 1:1:N g_1(x,y) = g_1(x,y) * (-1)^(x+y); g_2(x,y) = g_2(x,y) * (-1)^(x+y); end end %% -----show------- figure(); subplot(1,2,1); imshow(f,[0 1]); xlabel('a).Original Image'); subplot(1,2,2); imshow(log(1 + abs(F)),[ ]); xlabel('b).Fourier spectrum of a'); figure(); subplot(1,2,1); imshow(H_1,[0 1]); xlabel('c).Ideal Highpass filter(D=60)'); subplot(1,2,2); h = mesh(1:20:P,1:20:Q,H_1(1:20:P,1:20:Q)); set(h,'EdgeColor','k'); axis([0 P 0 Q 0 1]); xlabel('u');ylabel('v'); zlabel('|H(u,v)|'); figure(); subplot(1,2,1); imshow(log(1 + abs(G_1)),[ ]); xlabel('d).Result of filtering using c'); subplot(1,2,2); imshow(g_1,[0 1]); xlabel('e).Result image'); figure(); subplot(1,2,1); imshow(H_2,[0 1]); xlabel('f).Ideal Highpass filter(D=160)'); subplot(1,2,2); h = mesh(1:20:P,1:20:Q,H_2(1:20:P,1:20:Q)); set(h,'EdgeColor','k'); axis([0 P 0 Q 0 1]); xlabel('u');ylabel('v'); zlabel('|H(u,v)|'); figure(); subplot(1,2,1); imshow(log(1 + abs(G_2)),[ ]); xlabel('g).Result of filtering using e'); subplot(1,2,2); imshow(g_2,[0 1]); xlabel('h).Result image'); close all; clear all; %% ---------Butterworth Highpass Filters (Fre. Domain)------------ f = imread('characters_test_pattern.tif'); f = mat2gray(f,[0 255]); [M,N] = size(f); P = 2*M; Q = 2*N; fc = zeros(M,N); for x = 1:1:M for y = 1:1:N fc(x,y) = f(x,y) * (-1)^(x+y); end end F = fft2(fc,P,Q); H_1 = zeros(P,Q); H_2 = zeros(P,Q); for x = (-P/2):1:(P/2)-1 for y = (-Q/2):1:(Q/2)-1 D = (x^2 + y^2)^(0.5); D_0 = 100; H_1(x+(P/2)+1,y+(Q/2)+1) = 1/(1+(D_0/D)^2); H_2(x+(P/2)+1,y+(Q/2)+1) = 1/(1+(D_0/D)^6); end end G_1 = H_1 .* F; G_2 = H_2 .* F; g_1 = real(ifft2(G_1)); g_1 = g_1(1:1:M,1:1:N); g_2 = real(ifft2(G_2)); g_2 = g_2(1:1:M,1:1:N); for x = 1:1:M for y = 1:1:N g_1(x,y) = g_1(x,y) * (-1)^(x+y); g_2(x,y) = g_2(x,y) * (-1)^(x+y); end end %% -----show------- figure(); subplot(1,2,1); imshow(f,[0 1]); xlabel('a).Original Image'); subplot(1,2,2); imshow(log(1 + abs(F)),[ ]); xlabel('b).Fourier spectrum of a'); figure(); subplot(1,2,1); imshow(H_1,[0 1]); xlabel('c)Butterworth Lowpass (D_{0}=100,n=1)'); subplot(1,2,2); h = mesh(1:20:P,1:20:Q,H_1(1:20:P,1:20:Q)); set(h,'EdgeColor','k'); axis([0 P 0 Q 0 1]); xlabel('u');ylabel('v'); zlabel('|H(u,v)|'); figure(); subplot(1,2,1); imshow(log(1 + abs(G_1)),[ ]); xlabel('d).Result of filtering using c'); subplot(1,2,2); imshow(g_1,[0 1]); xlabel('e).Result image'); figure(); subplot(1,2,1); imshow(H_2,[0 1]); xlabel('f).Butterworth Lowpass (D_{0}=100,n=3)'); subplot(1,2,2); h = mesh(1:20:P,1:20:Q,H_2(1:20:P,1:20:Q)); set(h,'EdgeColor','k'); axis([0 P 0 Q 0 1]); xlabel('u');ylabel('v'); zlabel('|H(u,v)|'); figure(); subplot(1,2,1); imshow(log(1 + abs(G_2)),[ ]); xlabel('g).Result of filtering using e'); subplot(1,2,2); imshow(g_2,[0 1]); xlabel('h).Result image'); close all; clear all; %% ---------Gaussian Highpass Filters (Fre. Domain)------------ f = imread('characters_test_pattern.tif'); f = mat2gray(f,[0 255]); [M,N] = size(f); P = 2*M; Q = 2*N; fc = zeros(M,N); for x = 1:1:M for y = 1:1:N fc(x,y) = f(x,y) * (-1)^(x+y); end end F = fft2(fc,P,Q); H_1 = zeros(P,Q); H_2 = zeros(P,Q); for x = (-P/2):1:(P/2)-1 for y = (-Q/2):1:(Q/2)-1 D = (x^2 + y^2)^(0.5); D_0 = 60; H_1(x+(P/2)+1,y+(Q/2)+1) = 1 - exp(-(D*D)/(2*D_0*D_0)); D_0 = 160; H_2(x+(P/2)+1,y+(Q/2)+1) = 1 - exp(-(D*D)/(2*D_0*D_0)); end end G_1 = H_1 .* F; G_2 = H_2 .* F; g_1 = real(ifft2(G_1)); g_1 = g_1(1:1:M,1:1:N); g_2 = real(ifft2(G_2)); g_2 = g_2(1:1:M,1:1:N); for x = 1:1:M for y = 1:1:N g_1(x,y) = g_1(x,y) * (-1)^(x+y); g_2(x,y) = g_2(x,y) * (-1)^(x+y); end end %% -----show------- figure(); subplot(1,2,1); imshow(f,[0 1]); xlabel('a).Original Image'); subplot(1,2,2); imshow(log(1 + abs(F)),[ ]); xlabel('b).Fourier spectrum of a'); figure(); subplot(1,2,1); imshow(H_1,[0 1]); xlabel('c)Gaussian Highpass (D_{0}=60)'); subplot(1,2,2); h = mesh(1:20:P,1:20:Q,H_1(1:20:P,1:20:Q)); set(h,'EdgeColor','k'); axis([0 P 0 Q 0 1]); xlabel('u');ylabel('v'); zlabel('|H(u,v)|'); figure(); subplot(1,2,1); imshow(log(1 + abs(G_1)),[ ]); xlabel('d).Result of filtering using c'); subplot(1,2,2); imshow(g_1,[0 1]); xlabel('e).Result image'); figure(); subplot(1,2,1); imshow(H_2,[0 1]); xlabel('f).Gaussian Highpass (D_{0}=160)'); subplot(1,2,2); h = mesh(1:20:P,1:20:Q,H_2(1:20:P,1:20:Q)); set(h,'EdgeColor','k'); axis([0 P 0 Q 0 1]); xlabel('u');ylabel('v'); zlabel('|H(u,v)|'); figure(); subplot(1,2,1); imshow(log(1 + abs(G_2)),[ ]); xlabel('g).Result of filtering using e'); subplot(1,2,2); imshow(g_2,[0 1]); xlabel('h).Result image');1.5 锐化滤波器
close all; clear all; clc; %% ---------The High-Fre-Emphasis Filters (Fre. Domain)------------ f = imread('blurry_moon.tif'); f = mat2gray(f,[0 255]); [M,N] = size(f); P = 2*M; Q = 2*N; fc = zeros(M,N); for x = 1:1:M for y = 1:1:N fc(x,y) = f(x,y) * (-1)^(x+y); end end F = fft2(fc,P,Q); H_HP = zeros(P,Q); for x = (-P/2):1:(P/2)-1 for y = (-Q/2):1:(Q/2)-1 D = (x^2 + y^2)^(0.5); D_0 = 80; H_HP(x+(P/2)+1,y+(Q/2)+1) = 1 - exp(-(D*D)/(2*D_0*D_0)); end end G_HP = H_HP .* F; G_HFE = (1 + 1.1 * H_HP) .* F; g_1 = real(ifft2(G_HP)); g_1 = g_1(1:1:M,1:1:N); g_2 = real(ifft2(G_HFE)); g_2 = g_2(1:1:M,1:1:N); for x = 1:1:M for y = 1:1:N g_1(x,y) = g_1(x,y) * (-1)^(x+y); g_2(x,y) = g_2(x,y) * (-1)^(x+y); end end g = histeq(g_2); %g_1 = mat2gray(g_1); %% -----show------- figure(); subplot(1,2,1); imshow(f,[0 1]); xlabel('a).Original Image'); subplot(1,2,2); imshow(log(1 + abs(F)),[ ]); xlabel('b).Fourier spectrum of a'); figure(); subplot(1,2,1); imshow(g_1,[0 1]); xlabel('c).Result image of High-pass Filter'); subplot(1,2,2); imshow(log(1 + abs(G_HP)),[ ]); xlabel('d).Result of filtering using e'); figure(); subplot(1,2,1); imshow(g_2,[0 1]); xlabel('e).Result image of High-Fre-Emphasis Filter'); subplot(1,2,2); imshow(log(1 + abs(G_HFE)),[ ]); xlabel('f).Result of filtering using e');2.带阻滤波器
close all; clear all; clc; %% ---------------------Add Noise------------------------- f = imread('left.tif'); f = mat2gray(f,[0 255]); [M,N] = size(f); P = 2*M; Q = 2*N; fc = zeros(M,N); for x = 1:1:M for y = 1:1:N fc(x,y) = f(x,y) * (-1)^(x+y); end end F = fft2(fc,P,Q); H_NP = ones(P,Q); for x = (-P/2):1:(P/2)-1 for y = (-Q/2):1:(Q/2)-1 D = 2; v_k = -200; u_k = 150; D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5); H_NP(x+(P/2)+1,y+(Q/2)+1) = H_NP(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4); D_k = ((x-u_k)^2 + (y-v_k)^2)^(0.5); H_NP(x+(P/2)+1,y+(Q/2)+1) = H_NP(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4); v_k = 200; u_k = 150; D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5); H_NP(x+(P/2)+1,y+(Q/2)+1) = H_NP(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4); D_k = ((x-u_k)^2 + (y-v_k)^2)^(0.5); H_NP(x+(P/2)+1,y+(Q/2)+1) = H_NP(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4); v_k = 0; u_k = 250; D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5); H_NP(x+(P/2)+1,y+(Q/2)+1) = H_NP(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4); D_k = ((x-u_k)^2 + (y-v_k)^2)^(0.5); H_NP(x+(P/2)+1,y+(Q/2)+1) = H_NP(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4); v_k = 250; u_k = 0; D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5); H_NP(x+(P/2)+1,y+(Q/2)+1) = H_NP(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4); D_k = ((x-u_k)^2 + (y-v_k)^2)^(0.5); H_NP(x+(P/2)+1,y+(Q/2)+1) = H_NP(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4); H_NP(x+(P/2)+1,y+(Q/2)+1) = 1 + 700*(1 - H_NP(x+(P/2)+1,y+(Q/2)+1)); end end G_Noise = F .* H_NP; g_noise = real(ifft2(G_Noise)); g_noise = g_noise(1:1:M,1:1:N); for x = 1:1:M for y = 1:1:N g_noise(x,y) = g_noise(x,y) * (-1)^(x+y); end end %% ---------Bondpass Filters (Fre. Domain)------------ H_1 = ones(P,Q); for x = (-P/2):1:(P/2)-1 for y = (-Q/2):1:(Q/2)-1 D = (x^2 + y^2)^(0.5); D_0 = 250; W = 30; H_1(x+(P/2)+1,y+(Q/2)+1) = 1/(1+((D*W)/((D*D) - (D_0*D_0)))^6); end end G_1 = H_1 .* G_Noise; g_1 = real(ifft2(G_1)); g_1 = g_1(1:1:M,1:1:N); for x = 1:1:M for y = 1:1:N g_1(x,y) = g_1(x,y) * (-1)^(x+y); end end %% -----show------- close all; figure(); subplot(1,2,1); imshow(f,[0 1]); xlabel('a).Original Image'); subplot(1,2,2); imshow(log(1 + abs(F)),[ ]); xlabel('b).Fourier spectrum of a'); figure(); subplot(1,2,2); imshow(log(1 + abs(G_Noise)),[ ]); xlabel('c).Fourier spectrum of b'); subplot(1,2,1); imshow(g_noise,[0 1]); xlabel('b).Result of add noise'); figure(); subplot(1,2,1); imshow(H_1,[0 1]); xlabel('d).Butterworth Bandpass filter(D=355 W=40 n =2)'); subplot(1,2,2); h = mesh(1:20:Q,1:20:P,H_1(1:20:P,1:20:Q)); set(h,'EdgeColor','k'); axis([0 Q 0 P 0 1]); xlabel('u');ylabel('v'); zlabel('|H(u,v)|'); figure(); subplot(1,2,2); imshow(log(1 + abs(G_1)),[ ]); xlabel('e).Fourier spectrum of f'); subplot(1,2,1); imshow(g_1,[0 1]); xlabel('f).Result of denoise');3.陷波滤波器(Notch Filter)
close all; clear all; clc; %% ---------Butterworth Notch filter (Fre. Domain)------------ f = imread('car_75DPI_Moire.tif'); f = mat2gray(f,[0 255]); [M,N] = size(f); P = 2*M; Q = 2*N; fc = zeros(M,N); for x = 1:1:M for y = 1:1:N fc(x,y) = f(x,y) * (-1)^(x+y); end end F = fft2(fc,P,Q); H_NF = ones(P,Q); for x = (-P/2):1:(P/2)-1 for y = (-Q/2):1:(Q/2)-1 D = 30; v_k = 59; u_k = 77; D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5); H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4); D_k = ((x-u_k)^2 + (y-v_k)^2)^(0.5); H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4); v_k = 59; u_k = 159; D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5); H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4); D_k = ((x-u_k)^2 + (y-v_k)^2)^(0.5); H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4); v_k = -54; u_k = 84; D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5); H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4); D_k = ((x-u_k)^2 + (y-v_k)^2)^(0.5); H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4); v_k = -54; u_k = 167; D_k = ((x+u_k)^2 + (y+v_k)^2)^(0.5); H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4); D_k = ((x-u_k)^2 + (y-v_k)^2)^(0.5); H_NF(x+(P/2)+1,y+(Q/2)+1) = H_NF(x+(P/2)+1,y+(Q/2)+1) * 1/(1+(D/D_k)^4); end end G_1 = H_NF .* F; g_1 = real(ifft2(G_1)); g_1 = g_1(1:1:M,1:1:N); for x = 1:1:M for y = 1:1:N g_1(x,y) = g_1(x,y) * (-1)^(x+y); end end %% -----show------- close all; figure(); subplot(1,2,1); imshow(f,[0 1]); xlabel('a).Original Image'); subplot(1,2,2); imshow(log(1 + abs(F)),[ ]); xlabel('b).Fourier spectrum of a'); figure(); subplot(1,2,1); imshow(H_NF,[0 1]); xlabel('c).Butterworth Notch filter(D=30 n=2)'); subplot(1,2,2); h = mesh(1:10:Q,1:10:P,H_NF(1:10:P,1:10:Q)); set(h,'EdgeColor','k'); axis([0 Q 0 P 0 1]); xlabel('u');ylabel('v'); zlabel('|H(u,v)|'); figure(); subplot(1,2,2); imshow(log(1 + abs(G_1)),[ ]); xlabel('e).Fourier spectrum of d'); subplot(1,2,1); imshow(g_1,[0 1]); xlabel('d).Result image');