【youcans 的 OpenCV 例程 300 篇】102. 陷波带阻滤波器的传递函数

通过频率域滤波可以有效分析并滤除周期噪声，其理论基础是傅里叶变换后周期噪声在对应周期干扰的频率显示为集中突发的能量，因此可以使用选择性滤波器来分离滤除噪声。

4.1 陷波滤波器（Notch Filter）

陷波滤波器阻止或通过预定的频率矩形邻域中的频率，可以很好地复原被周期性噪声干扰的图像。在《5.2 陷波滤波器（Notch Filter）》已进行介绍并给出了例程。

陷波滤波器可以在某一个频率点迅速衰减输入信号，以达到阻碍此频率信号通过的滤波效果的滤波器。

陷波带阻滤波器的传递函数是中心平移到陷波中心的各个高通滤波器的乘积：
$HNR(u,v)=∏k=1QHk(u,v)H−k(u,v)H_{NR}(u,v) = \prod_{k=1}^Q H_k(u,v) H_{-k}(u,v)$

其中，滤波器的距离计算公式为：
$Dk(u,v)=(u−M/2−uk)2+(v−N/2−vk)2D−k(u,v)=(u−M/2+uk)2+(v−N/2+vk)2D_k(u,v) = \sqrt{(u-M/2-u_k)^2 + (v-N/2-v_k)^2} \\ D_{-k}(u,v) = \sqrt{(u-M/2+u_k)^2 + (v-N/2+v_k)^2}$
例如，具有 3个陷波对的 n 阶巴特沃斯陷波带阻滤波器为：
$HNR(u,v)=∏k=13[11+[D0k/Dk(u,v)]n][11+[D−k/Dk(u,v)]n]H_{NR}(u,v) = \prod_{k=1}^3 [\frac{1}{1+[D_{0k}/D_k(u,v)]^n}] [\frac{1}{1+[D_{-k}/D_k(u,v)]^n}]$

例程 9.16：陷波带阻滤波器的传递函数

    # 9.16: 陷波带阻滤波器的传递函数def ideaBondResistFilter(shape, radius=10, w=5):  # 理想带阻滤波器u, v = np.meshgrid(np.arange(shape[1]), np.arange(shape[0]))D = np.sqrt((u - shape[1]//2)**2 + (v - shape[0]//2)**2)D0 = radiushalfW = w/2kernel = np.piecewise(D, [D<=D0+halfW, D<=D0-halfW], [1, 0])kernel = 1 - kernel  # 带阻return kerneldef butterworthNRFilter(img, radius=10, uk=60, vk=80, n=1):  # 巴特沃斯陷波带阻滤波器M, N = img.shape[1], img.shape[0]u, v = np.meshgrid(np.arange(M), np.arange(N))Dkm = np.sqrt((u - M//2 - uk)**2 + (v - N//2 - vk)**2)  # D_+kDkp = np.sqrt((u - M//2 + uk)**2 + (v - N//2 + vk)**2)  # D_-kD0 = radiusn2 = n * 2epsilon = 1e-6kernel = (1 / (1 + (D0 / (Dkm+epsilon))**n2)) * (1 / (1 + (D0 / (Dkp+epsilon))**n2))return kerneldef gaussNRFilter(img, radius=10, uk=60, vk=80):  # 高斯陷波带阻滤波器M, N = img.shape[1], img.shape[0]u, v = np.meshgrid(np.arange(M), np.arange(N))Dkm = np.sqrt((u - M//2 - uk)**2 + (v - N//2 - vk)**2)  # D_+kDkp = np.sqrt((u - M//2 + uk)**2 + (v - N//2 + vk)**2)  # D_-kD0 = radiuskernel = (1 - np.exp(-(Dkm**2)/(D0**2))) * (1 - np.exp(-(Dkp**2)/(D0**2)))return kerneldef ideaNRFilter(img, radius=10, uk=60, vk=80):  # 高斯陷波带阻滤波器M, N = img.shape[1], img.shape[0]u, v = np.meshgrid(np.arange(M), np.arange(N))Dkm = np.sqrt((u - M//2 - uk)**2 + (v - N//2 - vk)**2)  # D_+kDkp = np.sqrt((u - M//2 + uk)**2 + (v - N//2 + vk)**2)  # D_-kD0 = radiusk1 = Dkm.copy()k1[Dkm>D0] = 1k1[Dkm<=D0] = 0k2 = Dkp.copy()k2[Dkp>D0] = 1k2[Dkp<=D0] = 0kernel = k1 * k2return kernel# 理想、高斯、巴特沃斯陷波带阻滤波器传递函数img = np.zeros([128, 128])shape = img.shaperadius = 15INRF = ideaNRFilter(img, radius=radius, uk=20, vk=30)  # (uk,vk) 陷波中心GNRF = gaussNRFilter(img, radius=radius, uk=20, vk=30)BNRF = butterworthNRFilter(img, radius=radius, uk=20, vk=30, n=2)filters = ["INRF", "GNRF", "BNRF"]u, v = np.mgrid[-1:1:2.0/shape[0], -1:1:2.0/shape[1]]fig = plt.figure(figsize=(10, 8))for i in range(3):nrFilter = eval(filters[i]).copy()ax1 = fig.add_subplot(3, 3, 3*i+1)ax1.imshow(nrFilter, 'gray')ax1.set_title(filters[i]), ax1.set_xticks([]), ax1.set_yticks([])ax2 = plt.subplot(3,3,3*i+2, projection='3d')ax2.set_title("transfer function")ax2.plot_wireframe(u, v, nrFilter, rstride=2, linewidth=0.5, color='c')ax2.set_xticks([]), ax2.set_yticks([]), ax2.set_zticks([])ax3 = plt.subplot(3,3,3*i+3)profile = nrFilter[:, 30]ax3.plot(profile), ax3.set_title("profile"), ax3.set_xticks([]), ax3.set_yticks([])plt.show()