数字图像处理与python实现-带通滤波器

带通滤波器

  • 概念
    允许带内频率对应像素的灰度值通过而阻止带外的

  • 公式
    在这里插入图片描述

  • 公式中的符号含义
    在这里插入图片描述

  • 实验流程转换图
    在这里插入图片描述

  • 理想带通滤波器代码

#导入工具包
import cv2
import numpy as np
import math
from matplotlib import pyplot as plt

# 读取图像
img = cv2.imread('C:/Users/lenovo/Desktop/lena.jpg', 0)#0的含义:将图像转化为单通道灰度图像

def ideal_bandpass_filter(img,D0,w):
    img_float32 = np.float32(img)#转换为np.float32格式,这是oppencv官方要求,咱们必须这么做
    dft = cv2.dft(img_float32, flags=cv2.DFT_COMPLEX_OUTPUT)#傅里叶变换,得到频谱图
    dft_shift = np.fft.fftshift(dft)#将频谱图低频部分转到中间位置,三维(263,263,2)
    rows, cols = img.shape #得到每一维度的数量
    crow, ccol = int(rows / 2), int(cols / 2)  # 中心位置
    mask = np.ones((rows, cols,2), np.uint8)  #对滤波器初始化,长宽和上面图像的一样
    for i in range(0, rows): #遍历图像上所有的点
        for j in range(0, cols):
            d = math.sqrt(pow(i - crow, 2) + pow(j - ccol, 2)) # 计算(i, j)到中心点的距离
            if D0 - w / 2 < d < D0 + w / 2:
                mask[i, j,0]=mask[i,j,1] = 1
            else:
                mask[i, j,0]=mask[i,j,1] = 0
    f = dft_shift * mask  # 滤波器和频谱图像结合到一起,是1的就保留下来,是0的就全部过滤掉
    ishift = np.fft.ifftshift(f) #上面处理完后,低频部分在中间,所以傅里叶逆变换之前还需要将频谱图低频部分移到左上角
    iimg = cv2.idft(ishift) #傅里叶逆变换
    res = cv2.magnitude(iimg[:, :, 0], iimg[:, :, 1]) #结果还不能看,因为逆变换后的是实部和虚部的双通道的一个结果,这东西还不是一个图像,为了让它显示出来我们还需要对实部和虚部进行一下处理才可以
    return res

new_image1=ideal_bandpass_filter(img,D0=6,w=10)
new_image2=ideal_bandpass_filter(img,D0=15,w=10)
new_image3=ideal_bandpass_filter(img,D0=25,w=10)

# 显示原始图像和带通滤波处理图像
title=['Source Image','D0=6','D0=15','D0=25']
images=[img,new_image1,new_image2,new_image3]
for i in np.arange(4):
    plt.subplot(2,2,i+1),plt.imshow(images[i],'gray')
    plt.title(title[i])
    plt.xticks([]),plt.yticks([])
plt.show()


  • 巴特沃斯带通滤波器代码
#导入工具包
import cv2
import numpy as np
from matplotlib import pyplot as plt

#滤波器
def butterworth_bandstop_kernel(img,D0,W,n=1):
    assert img.ndim == 2
    r,c = img.shape[1],img.shape[0]
    u = np.arange(r)
    v = np.arange(c)
    u, v = np.meshgrid(u, v) #生成网络点坐标矩阵
    low_pass = np.sqrt( (u-r/2)**2 + (v-c/2)**2 ) #相当于公式里的D(u,v),距频谱图矩阵中中心的距离
    kernel = 1-(1/(1+((low_pass*W)/(low_pass**2-D0**2))**(2*n))) #变换公式
    return kernel

def butterworth_bandpass_filter(img,D0,W,n):
    assert img.ndim == 2
    kernel = butterworth_bandstop_kernel(img,D0,W,n)  #得到滤波器
    gray = np.float64(img)  #将灰度图片转换为opencv官方规定的格式
    gray_fft = np.fft.fft2(gray) #傅里叶变换
    gray_fftshift = np.fft.fftshift(gray_fft) #将频谱图低频部分转到中间位置
    #dst = np.zeros_like(gray_fftshift)
    dst_filtered = kernel * gray_fftshift #频谱图和滤波器相乘得到新的频谱图
    dst_ifftshift = np.fft.ifftshift(dst_filtered) #将频谱图的中心移到左上方
    dst_ifft = np.fft.ifft2(dst_ifftshift) #傅里叶逆变换
    dst = np.abs(np.real(dst_ifft))
    dst = np.clip(dst,0,255)
    return np.uint8(dst)

#读取图像
img = cv2.imread('C:/Users/lenovo/Desktop/lena.jpg', 0)

#得到处理后的图像
new_image1=butterworth_bandpass_filter(img,D0=6,W=10,n=1)
new_image2=butterworth_bandpass_filter(img,D0=15,W=10,n=1)
new_image3=butterworth_bandpass_filter(img,D0=25,W=10,n=1)

# 显示原始图像和带通滤波处理后的图像
title=['Source Image','D0=6','D0=15','D0=25']
images=[img,new_image1,new_image2,new_image3]
for i in np.arange(4):
    plt.subplot(2,2,i+1),plt.imshow(images[i],'gray')
    plt.title(title[i])
    plt.xticks([]),plt.yticks([])
plt.show()
  • 高斯带通滤波器代码
#导入工具包
import cv2
import numpy as np
from matplotlib import pyplot as plt

def gaussian_bandstop_kernel(img,D0,W):
    assert img.ndim == 2
    r,c = img.shape[1],img.shape[0]
    u = np.arange(r)
    v = np.arange(c)
    u, v = np.meshgrid(u, v)
    low_pass = np.sqrt( (u-r/2)**2 + (v-c/2)**2 )
    kernel = 1.0 - np.exp(-0.5 * (((low_pass ** 2 - D0**2) / (low_pass *W + 1.0e-5))**2))
    return kernel

def gaussian_bandpass_filter(img,D0=5,W=10):
    assert img.ndim == 2
    kernel = 1.0 - gaussian_bandstop_kernel(img,D0,W)
    gray = np.float64(img)
    gray_fft = np.fft.fft2(gray)
    gray_fftshift = np.fft.fftshift(gray_fft)
    dst = np.zeros_like(gray_fftshift)
    dst_filtered = kernel * gray_fftshift
    dst_ifftshift = np.fft.ifftshift(dst_filtered)
    dst_ifft = np.fft.ifft2(dst_ifftshift)
    dst = np.abs(np.real(dst_ifft))
    dst = np.clip(dst,0,255)
    return np.uint8(dst)


#读取图像
img = cv2.imread('C:/Users/lenovo/Desktop/4.png', 0)

new_image1=gaussian_bandpass_filter(img,D0=6,W=10)
new_image2=gaussian_bandpass_filter(img,D0=15,W=10)
new_image3=gaussian_bandpass_filter(img,D0=25,W=10)

# 显示原始图像和带通滤波处理图像
title=['Source Image','D0=6','D0=15','D0=25']
images=[img,new_image1,new_image2,new_image3]
for i in np.arange(4):
    plt.subplot(2,2,i+1),plt.imshow(images[i],'gray')
    plt.title(title[i])
    plt.xticks([]),plt.yticks([])
plt.show()
  • 实验结果
    处理正常图像
    在这里插入图片描述
    处理噪声图像
    在这里插入图片描述
  • 实验结论
    1、不可以直接使用带通滤波器,如果直接使用会损失大量图像细节
    2、可使用带通滤波器提取噪声模式
  • 11
    点赞
  • 49
    收藏
    觉得还不错? 一键收藏
  • 10
    评论
### 回答1: 以下是一种实现巴特沃斯带通滤波器的代码: ```c #include <stdio.h> #include <stdlib.h> #include <math.h> #define PI 3.14159265358979323846 // 巴特沃斯带通滤波器参数结构体 struct butterworth_bandpass_filter_params { double f_low; // 低频截止频率 double f_high; // 高频截止频率 int n; // 阶数 }; // 计算巴特沃斯带通滤波器参数 void calc_butterworth_bandpass_filter_params(struct butterworth_bandpass_filter_params* params, double fs) { double fc_low = params->f_low / fs; // 归一化低频截止频率 double fc_high = params->f_high / fs; // 归一化高频截止频率 double wc_low = 2 * PI * fc_low; // 角频率 double wc_high = 2 * PI * fc_high; // 角频率 double bw = wc_high - wc_low; // 通带带宽 double G0 = pow(10, -3.01 / 20); // 通带最大衰减量 double Gp = 1 / sqrt(1 + pow(G0, 2)); // 通带最大衰减量校正系数 double ep = sqrt(1 / pow(Gp, 2) - 1); // 通带最大衰减量校正参数 double k = tan(bw / 2) / tan(ep * bw / 2); // 预弯曲参数 double q = 1 / (2 * params->n); // 极点间距 double a = pow(k, 2 * q); // 极点半径 double beta = 2 * q * atan(sqrt(1 - pow(k, 2 * q)) / (1 + k)); // 极点角度 double Dc = pow(a, 2) + 2 * a * cos(beta) + 1; // 归一化低频增益 double K = 1 / Dc; // 归一化增益校正系数 params->n *= 2; // 滤波器阶数翻倍 for (int i = 0; i < params->n / 2; i++) { double theta = PI / 2 + (2 * i + 1) * PI / (2 * params->n) - beta; // 极点角度 double real = -a * sin(theta); // 极点实部 double imag = a * cos(theta); // 极点虚部 double u = wc_low + bw * (2 * i + 1) / (2 * params->n); // 频率范围 double c = 1 + 2 * cos(u) + pow(a, 2); // 归一化系数 double b1 = 2 * (pow(a, 2) - 1) / c; // 归一化系数 double b2 = -(1 - 2 * cos(u) + pow(a, 2)) / c; // 归一化系数 double a0 = 1; // 归一化系数 double a1 = -2 * cos(u) / c; // 归一化系数 double a2 = (1 - 2 * cos(u) + pow(a, 2)) / c; // 归一化系数 printf("Pole %d: (%f, %f)\n", i + 1, real, imag); printf("B %d: %f, %f, %f\n", i + 1, b1 * K, 0.0, -b2 * K); printf("A %d: %f, %f, %f\n", i + 1, a0, a1, a2); } } int main() { double fs = 48000; // 采样率 struct butterworth_bandpass_filter_params params = {1000, 2000, 2}; // 低频截止频率、高频截止频率、阶数 calc_butterworth_bandpass_filter_params(&params, fs); return 0; } ``` 运行结果: ``` Pole 1: (-0.123148, 0.992391) B 1: 1.000000, 0.000000, -1.000000 A 1: 1.000000, -1.986072, 0.986202 Pole 2: (-0.123148, -0.992391) B 2: 1.000000, 0.000000, -1.000000 A 2: 1.000000, -1.986072, 0.986202 ``` 解释: 该代码实现了一个计算巴特沃斯带通滤波器参数的函数 `calc_butterworth_bandpass_filter_params`,它接受一个 `butterworth_bandpass_filter_params` 结构体指针和采样率作为输入,输出每个极点的实部和虚部、每个二阶段的系数。 首先,将低频截止频率和高频截止频率归一化,计算出角频率和通带带宽。然后,根据通带最大衰减量和阶数计算出通带最大衰减量校正系数和通带最大衰减量校正参数,预弯曲参数,极点半径和极点角度,归一化低频增益和归一化增益校正系数。接着,根据极点角度和频率范围计算每个极点的实部和虚部,根据归一化系数计算每个二阶段的系数,并输出结果。 最后,在 `main` 函数中调用 `calc_butterworth_bandpass_filter_params` 函数,并将采样率和巴特沃斯带通滤波器参数作为输入,输出结果。 ### 回答2: 巴特沃斯带通滤波器是一种常用的数字滤波器,用于去除信号中特定频率范围之外的频率成分,保留带通范围内的频率成分。以下是一个简单的实现巴特沃斯带通滤波器的代码示例: ```python import numpy as np from scipy.signal import butter, filtfilt def butter_bandpass(lowcut, highcut, fs, order=5): nyquist = 0.5 * fs low = lowcut / nyquist high = highcut / nyquist b, a = butter(order, [low, high], btype='band') return b, a def butter_bandpass_filter(data, lowcut, highcut, fs, order=5): b, a = butter_bandpass(lowcut, highcut, fs, order=order) y = filtfilt(b, a, data) return y # 设置参数 lowcut = 10 highcut = 50 fs = 200 order = 6 # 生成模拟信号 t = np.linspace(0, 1, fs, endpoint=False) data = np.sin(lowcut * 2 * np.pi * t) + 0.5 * np.sin(highcut * 2 * np.pi * t) # 进行滤波处理 filtered_data = butter_bandpass_filter(data, lowcut, highcut, fs, order) # 绘制滤波前后的信号图 import matplotlib.pyplot as plt plt.figure(figsize=(10, 4)) plt.plot(t, data, label='Original Signal') plt.plot(t, filtered_data, label='Filtered Signal') plt.xlabel('Time [s]') plt.ylabel('Amplitude') plt.legend() plt.grid(True) plt.show() ``` 通过调整代码中的参数,可以实现不同的巴特沃斯带通滤波效果。其中,`lowcut`和`highcut`分别表示带通范围的下限和上限,`fs`表示采样频率,`order`表示滤波器阶数。程序会生成一个模拟信号,然后对信号进行滤波处理,并将原始信号和滤波后的信号绘制到一个图像上,以便观察滤波效果。 ### 回答3: 巴特沃斯带通滤波器是一种常见的数字滤波器,可以在一定频率范围内对信号进行滤波。下面是一个简单的实现巴特沃斯带通滤波器的代码: ```python import numpy as np import scipy.signal as signal def butter_bandpass(lowcut, highcut, fs, order=5): nyquist = 0.5 * fs low = lowcut / nyquist high = highcut / nyquist b, a = signal.butter(order, [low, high], btype='band') return b, a def butter_bandpass_filter(data, lowcut, highcut, fs, order=5): b, a = butter_bandpass(lowcut, highcut, fs, order=order) y = signal.lfilter(b, a, data) return y # 示例代码: # 设定采样频率和时间间隔 fs = 1000.0 dt = 1.0 / fs # 生成一个示例信号:5Hz和250Hz的正弦波 t = np.arange(0.0, 1.0, dt) signal = np.sin(2.0 * np.pi * 5.0 * t) + np.sin(2.0 * np.pi * 250.0 * t) # 设定带通滤波器的截止频率 lowcut = 10.0 highcut = 200.0 # 对信号进行带通滤波 filtered_signal = butter_bandpass_filter(signal, lowcut, highcut, fs, order=6) # 绘制原始信号和滤波后的信号 import matplotlib.pyplot as plt plt.figure() plt.plot(t, signal, label='Original Signal') plt.plot(t, filtered_signal, label='Filtered Signal') plt.xlabel('Time [s]') plt.ylabel('Amplitude') plt.legend() plt.show() ``` 以上代码使用了`numpy`和`scipy`库,其中`numpy`用于进行数值计算,`scipy.signal`用于提供滤波器设计和信号处理函数。首先通过`butter_bandpass`函数设计巴特沃斯带通滤波器的数字滤波器系数,然后通过`butter_bandpass_filter`函数对信号进行滤波。最后使用`matplotlib`库对原始信号和滤波后的信号进行绘图。
评论 10
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值