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

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

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