MATLAB中的IIR滤波器设计

news2025/2/3 11:34:53

在数字信号处理中，滤波器是消除噪声、提取特征或调整信号频率的核心工具。其中，无限脉冲响应（IIR）滤波器因其低阶数实现陡峭滚降的特性，被广泛应用于音频处理、通信系统和生物医学工程等领域。借助MATLAB强大的工具箱，设计IIR滤波器变得直观且高效。本文将结合MATLAB代码，系统介绍IIR滤波器的设计原理与实现方法。

一、IIR滤波器设计原理

IIR滤波器的核心思想是将模拟滤波器转换为数字滤波器。设计流程通常分为两步：

选择模拟原型滤波器（如巴特沃斯、切比雪夫等）。
通过双线性变换或脉冲响应不变法，将模拟滤波器转换为数字滤波器。

相较于FIR滤波器，IIR滤波器的相位非线性是其缺点，但其计算效率高（低阶数实现高性能），尤其适合实时处理场景。

二、MATLAB中的IIR滤波器类型与设计

1. 巴特沃斯滤波器（Butterworth）

特点：通带和阻带均无波纹，过渡带较宽，适合对相位失真不敏感的场景（如音频滤波）。
MATLAB函数：butter

% 设计一个5阶低通巴特沃斯滤波器，截止频率100Hz（采样率1000Hz）
fs = 1000;          % 采样率
fc = 100;           % 截止频率
n = 5;              % 阶数
[b, a] = butter(n, fc/(fs/2), 'low'); 
freqz(b, a);        % 绘制频率响应

2. 切比雪夫I型滤波器（Chebyshev Type I）

特点：通带等波纹，阻带单调衰减，适合需要快速过渡但允许通带波动的场景（如雷达信号处理）。
MATLAB函数：cheby1

% 设计4阶高通滤波器，通带波纹1dB，截止频率200Hz
rp = 1;             % 通带波纹（dB）
[b, a] = cheby1(4, rp, 200/(fs/2), 'high');

3. 切比雪夫II型滤波器（Chebyshev Type II）

特点：阻带等波纹，通带单调衰减，适合抑制特定频段干扰（如工频噪声消除）。
MATLAB函数：cheby2

% 设计6阶带阻滤波器，阻带衰减40dB，阻带范围150-250Hz
rs = 40;            % 阻带衰减（dB）
[b, a] = cheby2(6, rs, [150 250]/(fs/2), 'stop');

4. 椭圆滤波器（Elliptic）

特点：通带和阻带均为等波纹，过渡带最窄，适合对性能要求苛刻的场合（如通信系统中的信道滤波）。
MATLAB函数：ellip

% 设计3阶带通椭圆滤波器，通带波纹0.5dB，阻带衰减50dB
rp = 0.5; rs = 50;
[b, a] = ellip(3, rp, rs, [50 150]/(fs/2), 'bandpass');

5. 贝塞尔滤波器（Bessel）

特点：最大程度保持相位线性，但过渡带较宽，适合需要信号波形保真的场景（如心电图信号处理）。
设计步骤：需先设计模拟滤波器，再转换为数字滤波器。

% 设计4阶模拟低通贝塞尔滤波器，转换为数字滤波器
[b_analog, a_analog] = besself(4, 100);  % 模拟设计
[b_digital, a_digital] = bilinear(b_analog, a_analog, fs); % 双线性变换

三、如何选择IIR滤波器？

巴特沃斯：通用场景，平衡性能与复杂度。
切比雪夫I型：需快速过渡且允许通带波纹。
切比雪夫II型：需抑制特定阻带频率。
椭圆滤波器：对过渡带宽度要求极高。
贝塞尔滤波器：重视相位线性，牺牲过渡带性能。

四、Matlab编程实现

% MATLAB程序：模拟信号加入单频噪声，并使用5个低通滤波器进行滤波

% 参数设置
fs = 1000;                % 采样频率 (Hz)
t = 0:1/fs:1-1/fs;        % 时间向量 (1秒)
f_signal = 40;            % 信号频率 (Hz)
f_noise = 250;            % 噪声频率 (Hz)

% 原始信号：一个40Hz的正弦波
signal = sin(2*pi*f_signal*t);

% 添加噪声：一个250Hz的单频噪声
noise = 0.5 * sin(2*pi*f_noise*t);

% 带噪声信号
noisy_signal = signal + noise;

% 设计滤波器：使用巴特沃斯、切比雪夫I型、切比雪夫II型、椭圆和贝塞尔滤波器
% 1. 巴特沃斯滤波器
[b_butter, a_butter] = butter(4, 0.2);  % 低通滤波器，截止频率0.2*fs

% 2. 切比雪夫I型滤波器
[b_cheby1, a_cheby1] = cheby1(4, 0.5, 0.2);  % 通带波纹0.5dB

% 3. 切比雪夫II型滤波器
[b_cheby2, a_cheby2] = cheby2(4, 40, 0.2);  % 阻带衰减40dB

% 4. 椭圆滤波器
[b_ellip, a_ellip] = ellip(4, 0.5, 40, 0.2);  % 通带波纹0.5dB，阻带衰减40dB

% 5. 贝塞尔滤波器
[b_bessel, a_bessel] = besself(4, 0.2);  % 贝塞尔滤波器，最大化相位线性

% 滤波处理
filtered_signal_butter = filter(b_butter, a_butter, noisy_signal);
filtered_signal_cheby1 = filter(b_cheby1, a_cheby1, noisy_signal);
filtered_signal_cheby2 = filter(b_cheby2, a_cheby2, noisy_signal);
filtered_signal_ellip = filter(b_ellip, a_ellip, noisy_signal);
filtered_signal_bessel = filter(b_bessel, a_bessel, noisy_signal);

% 频谱分析：使用FFT计算信号的频谱
n = length(t);                        % 信号长度
f = (0:n-1)*(fs/n);                    % 频率向量
signal_fft = abs(fft(noisy_signal));   % 噪声信号的频谱
butter_fft = abs(fft(filtered_signal_butter));
cheby1_fft = abs(fft(filtered_signal_cheby1));
cheby2_fft = abs(fft(filtered_signal_cheby2));
ellip_fft = abs(fft(filtered_signal_ellip));
bessel_fft = abs(fft(filtered_signal_bessel));

% 绘制结果
figure;

% 时域图：显示原始信号与不同滤波器处理后的信号
subplot(2, 2, 1);
plot(t, noisy_signal, 'k'); hold on;
plot(t, filtered_signal_butter);
title('时域 - 巴特沃斯滤波');
legend('带噪声信号', '巴特沃斯滤波结果');

subplot(2, 2, 2);
plot(t, noisy_signal, 'k'); hold on;
plot(t, filtered_signal_cheby1);
title('时域 - 切比雪夫I型滤波');
legend('带噪声信号', '切比雪夫I型滤波结果');

subplot(2, 2, 3);
plot(t, noisy_signal, 'k'); hold on;
plot(t, filtered_signal_cheby2);
title('时域 - 切比雪夫II型滤波');
legend('带噪声信号', '切比雪夫II型滤波结果');

subplot(2, 2, 4);
plot(t, noisy_signal, 'k'); hold on;
plot(t, filtered_signal_ellip);
title('时域 - 椭圆滤波');
legend('带噪声信号', '椭圆滤波结果');

% 频域图：显示原始信号与不同滤波器处理后的频谱
figure;

subplot(2, 2, 1);
plot(f, signal_fft, 'k'); hold on;
plot(f, butter_fft);
title('频域 - 巴特沃斯滤波');
legend('带噪声信号', '巴特沃斯滤波结果');

subplot(2, 2, 2);
plot(f, signal_fft, 'k'); hold on;
plot(f, cheby1_fft);
title('频域 - 切比雪夫I型滤波');
legend('带噪声信号', '切比雪夫I型滤波结果');

subplot(2, 2, 3);
plot(f, signal_fft, 'k'); hold on;
plot(f, cheby2_fft);
title('频域 - 切比雪夫II型滤波');
legend('带噪声信号', '切比雪夫II型滤波结果');

subplot(2, 2, 4);
plot(f, signal_fft, 'k'); hold on;
plot(f, ellip_fft);
title('频域 - 椭圆滤波');
legend('带噪声信号', '椭圆滤波结果');