巴特沃斯滤波器

---

一、MATLAB 实现

1. 巴特沃斯滤波器函数（支持图像/信号）

function H = butterworth_filter(D0, size, n, mode)
% BUTTERWORTH_FILTER 生成巴特沃斯滤波器
%   - D0: 截止频率
%   - size: 滤波器尺寸（图像：[height, width]，信号：[1, length]）
%   - n: 阶数（默认1）
%   - mode: 'low'（低通）或 'high'（高通，默认低通）

if nargin < 3, n = 1; end
if nargin < 4, mode = 'low'; end

% 生成中心对称坐标网格
[x, y] = meshgrid(1:size(2), 1:size(1));
x_centered = x - floor(size(2)/2) - 1;
y_centered = y - floor(size(1)/2) - 1;
D = sqrt(x_centered.^2 + y_centered.^2);

% 计算滤波器响应
H = 1 ./ (1 + (D/D0).^(2*n));

% 高通模式取反
if strcmpi(mode, 'high')
    H = 1 - H;
end
end

2. 频域中心化函数

function data_centered = centralize(data)
% CENTRALIZE 频域中心化（图像/信号通用）
%   - data: 输入频域数据（复数矩阵）

mask = (-1).^(meshgrid(1:size(data,2), 1:size(data,1)) + ...
            meshgrid(1:size(data,1), 1:size(data,2))');
data_centered = data .* mask;
end

3. 使用示例（图像处理）

% 读取图像并转为灰度
img = im2double(rgb2gray(imread('lena.jpg')));

% 傅里叶变换并中心化
F = centralize(fft2(img));

% 生成低通滤波器（D0=50, n=2）
H = butterworth_filter(50, size(img), 2, 'low');

% 滤波并反变换
filtered_F = F .* H;
filtered_img = real(ifft2(centralize(filtered_F)));

% 显示结果
imshowpair(img, filtered_img, 'montage');

---

二、Python 实现（使用 NumPy 和 OpenCV）

1. 巴特沃斯滤波器函数

import numpy as np

def butterworth_filter(D0, size, n=1, mode='low'):
    """
    生成巴特沃斯滤波器（支持2D图像/1D信号）
    - D0: 截止频率
    - size: 滤波器尺寸（图像：(height, width)，信号：(length,)）
    - n: 阶数
    - mode: 'low'或'high'
    """
    # 生成坐标网格
    if len(size) == 2:
        rows, cols = size
        y, x = np.ogrid[-rows//2:rows//2, -cols//2:cols//2]
    else:  # 1D信号
        length = size[0]
        x = np.arange(-length//2, length//2)
        y = 0
    
    D = np.sqrt(x**2 + y**2)
    D = np.where(D == 0, 1e-6, D)  # 避免除以零
    
    # 计算滤波器响应
    H = 1 / (1 + (D/D0)**(2*n))
    
    # 高通模式
    if mode.lower() == 'high':
        H = 1 - H
    return H

2. 频域中心化函数

def centralize(data):
    """
    频域中心化（支持2D图像/1D信号）
    - data: 输入频域数据（复数矩阵）
    """
    if data.ndim == 2:
        m, n = data.shape
        mask = (-1)**(np.arange(m)[:, None] + np.arange(n))
    else:  # 1D信号
        length = data.shape[0]
        mask = (-1)**np.arange(length)
    return data * mask

3. 使用示例（信号处理）

import matplotlib.pyplot as plt

# 生成含噪声的1D信号
t = np.linspace(0, 1, 1000)
signal = np.sin(2*np.pi*5*t) + 0.5*np.random.randn(1000)

# 傅里叶变换并中心化
F = centralize(np.fft.fft(signal))

# 生成高通滤波器（D0=20, n=3）
H = butterworth_filter(20, (1000,), n=3, mode='high')

# 滤波并反变换
filtered_F = F * H
filtered_signal = np.real(np.fft.ifft(centralize(filtered_F)))

# 绘制结果
plt.figure()
plt.plot(t, signal, label='Original')
plt.plot(t, filtered_signal, label='Filtered')
plt.legend()

---