4.10 选择性滤波器

本节所用原始图片素材共享于原始图片，提取码 523o
数字图像处理matlab版给出了低通滤波器和高通滤波器的代码，但带阻滤波器和陷波器的代码没有给出来，本节基于高低通滤波器的代码，结合带阻滤波器、陷波滤波器的表达式给出这两种滤波器的matlab实现

带阻滤波器

带阻滤波器的的表达式如下，其中D(u,v)D(u,v)D(u,v)是距离频率中心点的距离，D0D_0D0是阻带中心频率， WWW是阻带宽度：

理想	巴特沃斯	高斯
H(u,v)={0若D0−W2≤D≤D0+W21其他H(u,v) = \begin{cases} 0 & 若D_0-\frac{W}{2} \le D \le D_0+\frac{W}{2} \cr 1 & 其他 \end{cases}H(u,v)={01若D0−2W≤D≤D0+2W其他	H(u,v)=11+[DWD2−D02]2nH(u,v)=\frac{1}{1+\left[ \frac{DW}{D^2 -D_0^2} \right]^{2n}}H(u,v)=1+[D2−D02DW]2n1	H(u,v)=1−e−[D2−D02DW]2H(u, v) = 1 - e^{-\left[ \frac{D^2-D_0^2}{DW} \right]^{2}}H(u,v)=1−e−[DWD2−D02]2

基于上面表格所列出的表达式，我们可以得到带阻滤波器的代码如下

function H = brfilter(type, M, N, D0, W, n)%BRFILTER Computes frequency domain band reject filters.
%   H = BRFILTER(TYPE, M, N, D0, W, n) creates the transfer function of
%   a band reject filter, H, of the specified TYPE and size (M-by-N). To
%   view the filter as an image or mesh plot, it should be centered
%   using H = fftshift(H). %
%   Valid values for TYPE, D0, W and n are:
%
%   'ideal'    Ideal band reject filter with cutoff frequency D0 and
%              reject band width W. n need not be supplied.
%              D0 and w must be positive.
%
%   'btw'      Butterworth band reject filter of order n, cutoff
%              D0 and reject band width W.  The default value for
%              n is 1.0.  D0 and W must be positive.
%
%   'gaussian' Gaussian band reject filter with cutoff (standard
%              deviation) D0 and reject band width W.  n need not
%              be supplied.  D0 and W must be positive. %   Author: caesar
%   Digital Image Processing
%   $Revision: 0.1 $  $Date: 2020/07/10 09:54:00 $% Use function dftuv to set up the meshgrid arrays needed for
% computing the required distances. [U, V] = dftuv(M, N);% Compute the distances D(U, V).
D = sqrt(U.^2 + V.^2);% Begin filter computations.
switch typecase 'ideal'H = double(D < (D0 - W/2)) + double(D > (D0 + W/2));case 'btw'if nargin == 5n = 1;   endH = 1./(1 + (D .* W ./ (D.^2-D0^2)).^(2*n));case 'gaussian'H = 1 - exp(-((D.^2 - D0^2)./(D.*W)).^2);otherwiseerror('Unknown filter type.')end

下图通过上面代码生成了一幅大小为512x512，阻带中心频率为150，阻带宽度为20的理想、巴特沃斯（4阶）、高斯带阻滤波器。

陷波滤波器

由傅里叶变换的对称性可知，零相移滤波器必须关于原点对称，一个中心位于(u0,v0)(u_0, v_0)(u0,v0)的陷波在位置(−u0,−v0)(-u_0,-v_0)(−u0,−v0)必须有一个对应的陷波，我们将其称为陷波对。
陷波带阻滤波器可以有若干个以(uk,vk)(u_k,v_k)(uk,vk)为中心的高通滤波器对的乘积来构造，这一点很容易解释：1、一个以(uk,vk)(u_k,v_k)(uk,vk)为中心的高通滤波器完成对频率(uk,vk)(u_k,v_k)(uk,vk)和(−u0,−v0)(-u_0,-v_0)(−u0,−v0)的筛选。2、频域将多个陷波对相乘，意味着空间域与这些陷波对逐个卷积，依次完成对各个陷波对中心频率的筛选。
陷波滤波器表达式如下
HNR(u,v)=∏k=1QHk(u,v)H−k(u,v)H_{NR}(u, v) = \prod_{k=1}^QH_k(u,v)H_{-k}(u,v)HNR(u,v)=k=1∏QHk(u,v)H−k(u,v)
其中Hk(u,v)H_k(u,v)Hk(u,v)和H−k(u,v)H_{-k}(u,v)H−k(u,v)是一对分别以(u0,v0)(u_0, v_0)(u0,v0)和(−u0,−v0)(-u_0,-v_0)(−u0,−v0)为中心的高通滤波器。以(u,v)(u,v)(u,v)为中心的理想、巴特沃斯、高斯高通滤波器的表达式分别如下，其中D0D_0D0是截止频率，n是巴特沃斯滤波器的阶数：

理想	巴特沃斯	高斯
H(u,v)={0若D(u,v)≤D01D(u,v)>D0H(u,v) = \begin{cases} 0 & 若D(u,v) \le D_0 \cr 1 & D(u,v) \gt D_0 \end{cases}H(u,v)={01若D(u,v)≤D0D(u,v)>D0	H(u,v)=11+[D0D(u,v)]2nH(u,v)=\frac{1}{1+\left[ \frac{D_0}{D(u, v)} \right]^{2n}}H(u,v)=1+[D(u,v)D0]2n1	H(u,v)=1−e−[D(u,v)2D0]2H(u, v) = 1 - e^{-\left[ \frac{D(u,v)}{2D_0} \right]^{2}}H(u,v)=1−e−[2D0D(u,v)]2

这样，我们就可以得到包含多个陷波对的陷波滤波器。例如
HNR(u,v)=∏k=13[11+[D0k/Dk(u,v)]2n][11+[D0k/D−k(u,v)]2n]H_{NR}(u,v)=\prod_{k=1}^3 \left[ \frac{1}{1+[D_{0k}/D_k(u,v)]^{2n}} \right]\left[ \frac{1}{1+[D_{0k}/D_{-k}(u,v)]^{2n}} \right]HNR(u,v)=k=1∏3[1+[D0k/Dk(u,v)]2n1][1+[D0k/D−k(u,v)]2n1]
就是一个包含三个陷波对的巴特波斯陷波滤波器，其中D0kD_{0k}D0k是每个陷波对的陷波带宽，不同的陷波对此参数可以不同。
基于以上讨论，可以得到陷波滤波器的matlab实现如下：

function NH = notchfilter(type, M, N, D, C, n)%NOTCHFILTER Computes frequency domain notch filters.%   H = NOTCHFILTER(TYPE, M, N, C, n) creates the transfer function of
%   a notch filter, H, of the specified TYPE and size (M-by-N).
%   Valid values for TYPE, D, C and n are:
%
%   'ideal'    Ideal notch filter with cutoff frequency Dk at C(Uk, Vk) and
%              C(-Uk, -Vk). n need not be supplied. D must be positive.
%
%   'btw'      Butterworth notch filter of order n, and cutoff
%              Dk at C(Uk, Vk) and (-Uk, -Vk).  The default value
%              for n is 1.0. D must be positive.
%
%   'gaussian' Gaussian notch filter with cutoff (standard
%              deviation) Dk at C(Uk, Vk) and C(-Uk, -Vk). n need
%              not be supplied. D must be positive.%   Author: caesar
%   Digital Image Processing
%   $Revision: 0.1 $  $Date: 2020/07/10 10:28:22 $% The transfer function NH at (U,V) of a notch filter is Hhp(U, V),
% where Hhp(U, V) is the transfer function of the corresponding highpass
% filter centered at (U, V).if nargin == 5n = 1; % Default value of n.end% K is the number of notch point
[K, ~] = size(C);% all the notch points have same cutoff
if length(D) == 1Dk = D * ones(1, K);
% each notch points have specified cutoff
elseif length(D) == KDk = D;
elseerror('Invalie length of D, which may be either a single number or a vector of same length as C(:,1)')
end[V, U] = dftuv(M, N);% Generate notch filter.
NH = 1;
switch typecase 'ideal'for k = 1:K% centered at(U, V)Dik = sqrt((U - C(k,1)).^2 + (V - C(k,2)).^2);        Hik = double(Dik >= Dk(k));% centered at(-U, -V)Di_k = sqrt((U + C(k,1)).^2 + (V + C(k,2)).^2);        Hi_k = double(Di_k >= Dk(k));NH = NH .* Hik .* Hi_k;endcase 'btw'if nargin == 5n = 1;    endfor k = 1:K% centered at(U, V)Dik = sqrt((U - C(k,1)).^2 + (V - C(k,2)).^2);      Hik = 1./(1 + (Dk(k)./Dik).^(2*n));% centered at(-U, -V)Di_k = sqrt((U + C(k,1)).^2 + (V + C(k,2)).^2);       Hi_k = 1./(1 + (Dk(k)./Di_k).^(2*n));NH = NH .* Hik .* Hi_k;endcase 'gaussian'for k = 1:K% centered at(U, V)Dik = sqrt((U - C(k,1)).^2 + (V - C(k,2)).^2);Hik = 1 - exp(-(Dik.^2)./(2*(Dk(k)^2)));     % centered at(-U, -V)Di_k = sqrt((U + C(k,1)).^2 + (V + C(k,2)).^2);        Hi_k = 1 - exp(-(Di_k.^2)./(2*(Dk(k)^2)));     NH = NH .* Hik .* Hi_k;endotherwiseerror('Unknown filter type.')end

下图显示了一幅大小为512x512，陷波点为[128, 64]、 [128, 192]、[-128, 64]、[-128, 192]，陷波带宽为20的理想、巴特沃斯(4阶)、高斯陷波滤波器：

上面示例所用代码如下：

C = [128, 64; 128, 192; -128, 64; -128, 192];H = notchfilter('ideal', 512, 512, 20, C);
figure
subplot(1, 2, 1)
mesh(fftshift(abs(H)))
title('理想陷波滤波器(透视图)')
subplot(1, 2, 2)
imshow(fftshift(abs(H)), [ ]);
title('理想陷波滤波器')H = notchfilter('btw', 512, 512, 20, C, 4);
figure
subplot(1, 2, 1)
mesh(fftshift(abs(H)))
title('巴特沃斯陷波滤波器(透视图)')
subplot(1, 2, 2)
imshow(fftshift(abs(H)), [ ]);
title('巴特沃斯陷波滤波器')H = notchfilter('gaussian', 512, 512, 20, C);
figure
subplot(1, 2, 1)
mesh(fftshift(abs(H)))
title('高斯陷波滤波器(透视图)')
subplot(1, 2, 2)
imshow(fftshift(abs(H)), [ ]);
title('高斯陷波滤波器')

下图(a)显示了一幅含有莫尔条纹的原始图像，图(b)显示了使用对数标定的图(a)的幅谱，从幅谱中可以看到导致这些莫尔条纹的频率点，为了去掉摩尔条纹，我们生成了一个四阶巴特沃斯陷波滤波器，如图(c)，通过频域处理再反变换回空间域，我们得到了如图(d)的不再含有莫尔条纹的图像。
上面关于去除莫尔条纹示例的代码如下：

D0 = 10 ;
n = 4;f = imread('Fig0464(a)(car_75DPI_Moire).tif');
[M, N] =  size(f);figure
subplot(1, 2, 1);
imshow(f);
title('(a)带莫尔条纹的原始图像')subplot(1, 2, 2);
F = fft2(f);
imshow(intrans(abs(F), 'log'), [ ])
title('(b)幅谱(对数标定)')C = [30, 38; 30, 81; -30, 38; -30, 81];
H = notchfilter('btw', M, N, 10, C, 4);
figure
subplot(1, 2, 1)
imshow(fftshift(abs(H)), [ ]);
title('(c)巴特沃斯陷波器')subplot(1, 2, 2)
G = H .* F;
g = real(ifft2(G));
imshow(g, [])
title('(d)使用4阶巴特沃斯陷波滤波器过滤结果')