OpenCV 例程200篇102. 陷波带阻滤波器的传递函数

Posted 2022-02-23 小白YouCans

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

欢迎关注 『OpenCV 例程200篇 100 篇』 系列，持续更新中
欢迎关注 『Python小白的OpenCV学习课』 系列，持续更新中

4. 频率域滤波图像复原

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

4.1 陷波滤波器（Notch Filter）

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

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

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

其中，滤波器的距离计算公式为：
$$\begin{gathered} D_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} \end{gathered}$$
例如，具有 3个陷波对的 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 = radius
        halfW = w/2
        kernel = np.piecewise(D, [D<=D0+halfW, D<=D0-halfW], [1, 0])
        kernel = 1 - kernel  # 带阻
        return kernel

    def 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_+k
        Dkp = np.sqrt((u - M//2 + uk)**2 + (v - N//2 + vk)**2)  # D_-k
        D0 = radius
        n2 = n * 2
        epsilon = 1e-6
        kernel = (1 / (1 + (D0 / (Dkm+epsilon))**n2)) * (1 / (1 + (D0 / (Dkp+epsilon))**n2))
        return kernel

    def 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_+k
        Dkp = np.sqrt((u - M//2 + uk)**2 + (v - N//2 + vk)**2)  # D_-k
        D0 = radius
        kernel = (1 - np.exp(-(Dkm**2)/(D0**2))) * (1 - np.exp(-(Dkp**2)/(D0**2)))
        return kernel

    def 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_+k
        Dkp = np.sqrt((u - M//2 + uk)**2 + (v - N//2 + vk)**2)  # D_-k
        D0 = radius

        k1 = Dkm.copy()
        k1[Dkm>D0] = 1
        k1[Dkm<=D0] = 0
        k2 = Dkp.copy()
        k2[Dkp>D0] = 1
        k2[Dkp<=D0] = 0

        kernel = k1 * k2
        return kernel

    # 理想、高斯、巴特沃斯陷波带阻滤波器传递函数
    img = np.zeros([128, 128])
    shape = img.shape
    radius = 15
    INRF = 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()

（本节完）

---

版权声明：

youcans@xupt 原创作品，转载必须标注原文链接

Copyright 2021 youcans, XUPT
Crated：2022-2-1

欢迎关注 『OpenCV 例程200篇 100 篇』 系列，持续更新中
欢迎关注 『Python小白的OpenCV学习课』 系列，持续更新中

【OpenCV 例程200篇】01. 图像的读取（cv2.imread）
【OpenCV 例程200篇】02. 图像的保存（cv2.imwrite）
【OpenCV 例程200篇】03. 图像的显示（cv2.imshow）
【OpenCV 例程200篇】04. 用 matplotlib 显示图像（plt.imshow）
【OpenCV 例程200篇】05. 图像的属性（np.shape）
【OpenCV 例程200篇】06. 像素的编辑（img.itemset）
【OpenCV 例程200篇】07. 图像的创建（np.zeros）
【OpenCV 例程200篇】08. 图像的复制（np.copy）
【OpenCV 例程200篇】09. 图像的裁剪（cv2.selectROI）
【OpenCV 例程200篇】10. 图像的拼接（np.hstack）
【OpenCV 例程200篇】11. 图像通道的拆分（cv2.split）
【OpenCV 例程200篇】12. 图像通道的合并（cv2.merge）
【OpenCV 例程200篇】13. 图像的加法运算（cv2.add）
【OpenCV 例程200篇】14. 图像与标量相加（cv2.add）
【OpenCV 例程200篇】15. 图像的加权加法（cv2.addWeight）
【OpenCV 例程200篇】16. 不同尺寸的图像加法
【OpenCV 例程200篇】17. 两张图像的渐变切换
youcans 的 OpenCV 例程 200 篇103. 陷波带阻滤波器消除周期噪声干扰

OpenCV 完整例程89. 带阻滤波器的传递函数

OpenCV 完整例程89. 带阻滤波器的传递函数

OpenCV 例程200篇100. 自适应局部降噪滤波器

OpenCV 例程200篇101. 自适应中值滤波器

OpenCV 例程200篇101. 自适应中值滤波器