引言
信号处理是现代工程和科学研究中不可或缺的一部分，而滤波器则是信号处理中最基本也最重要的工具之一！在各种噪声干扰的现实世界中，我们常常需要从"杂乱无章"的信号中提取有用信息，这正是滤波器大显身手的时刻。
MATLAB作为工程计算和数据分析的强大平台，提供了丰富而灵活的滤波器设计与实现工具。无论你是信号处理初学者，还是经验丰富的工程师，掌握MATLAB中滤波器的使用都能显著提升你处理实际问题的能力。
今天，我们将一起深入探讨MATLAB中三种最常用的滤波器类型：低通滤波器、带通滤波器和高通滤波器。我们不仅会了解它们的基本原理，还会通过实际的MATLAB代码示例来展示如何设计和应用这些滤波器。
滤波器基础知识
在开始具体实现之前，让我们先简单回顾一下滤波器的基本概念（非常重要）。
滤波器本质上是一种选择性地通过或抑制特定频率成分的系统。根据其功能，常见的滤波器可分为：
滤波器设计主要考虑两个域：时域和频域。在MATLAB中，我们可以使用多种方法设计这些滤波器，最常见的包括：
- 有限脉冲响应(FIR)滤波器设计
- 无限脉冲响应(IIR)滤波器设计
- 基于窗函数的滤波器设计
接下来，我们将使用MATLAB的Signal Processing Toolbox中的工具来实现这些滤波器，并通过实例来观察它们的性能。
准备工作
首先，确保你已安装MATLAB及其Signal Processing Toolbox。我们将使用以下函数：
- butter, cheby1, ellip – 设计IIR滤波器
- fir1, firls, firpm – 设计FIR滤波器
- freqz – 计算并显示滤波器的频率响应
- filter – 应用滤波器到信号
在开始之前，让我们创建一个测试信号，包含多个频率成分：
% 设置采样参数
fs = 1000; % 采样频率，Hz
t = 0:1/fs:1-1/fs; % 时间向量，1秒
% 创建包含多个频率成分的信号
f1 = 50; % 50Hz成分
f2 = 120; % 120Hz成分
f3 = 350; % 350Hz成分
% 合成信号
x = sin(2*pi*f1*t) + 0.5*sin(2*pi*f2*t) + 0.25*sin(2*pi*f3*t);
% 添加一些高斯噪声
noise = 0.2*randn(size(t));
x_noisy = x + noise;
% 绘制原始信号和含噪声信号
figure;
subplot(2,1,1);
plot(t, x);
title('原始信号');
xlabel('时间 (s)');
ylabel('幅度');
subplot(2,1,2);
plot(t, x_noisy);
title('含噪声信号');
xlabel('时间 (s)');
ylabel('幅度');
现在，我们有了一个包含50Hz、120Hz和350Hz三个频率成分，并且添加了一些高斯噪声的测试信号。接下来，我们将设计各种滤波器来处理这个信号。
低通滤波器设计与应用
低通滤波器允许低于截止频率的信号通过，而衰减高于截止频率的信号。这对于去除高频噪声特别有用。
使用Butterworth滤波器
Butterworth滤波器是一种经典的IIR滤波器，特点是通带内相当平坦，无波纹。
% 设计Butterworth低通滤波器
cutoff_freq = 100; % 截止频率100Hz
order = 6; % 滤波器阶数
Wn = cutoff_freq/(fs/2); % 归一化截止频率
[b, a] = butter(order, Wn, 'low');
% 应用滤波器
filtered_signal = filter(b, a, x_noisy);
% 绘制频率响应
figure;
freqz(b, a, 1024, fs);
title('Butterworth低通滤波器频率响应');
% 绘制时域信号对比
figure;
subplot(3,1,1);
plot(t, x);
title('原始无噪声信号');
xlabel('时间 (s)');
ylabel('幅度');
subplot(3,1,2);
plot(t, x_noisy);
title('含噪声信号');
xlabel('时间 (s)');
ylabel('幅度');
subplot(3,1,3);
plot(t, filtered_signal);
title('低通滤波后信号');
xlabel('时间 (s)');
ylabel('幅度');
在这个例子中，我们设计了一个截止频率为100Hz的低通滤波器。理论上，这应该保留50Hz的成分，大部分保留120Hz的成分，而显著衰减350Hz的成分和高频噪声。
使用FIR滤波器
FIR滤波器具有线性相位特性，对于某些应用非常重要。
% 设计FIR低通滤波器
order = 50; % 滤波器阶数（高一些以获得更好的性能）
Wn = cutoff_freq/(fs/2); % 归一化截止频率
b = fir1(order, Wn, 'low', hamming(order+1));
a = 1; % FIR滤波器的分母系数为1
% 应用滤波器
filtered_signal_fir = filter(b, a, x_noisy);
% 绘制频率响应
figure;
freqz(b, a, 1024, fs);
title('FIR低通滤波器频率响应');
% 绘制结果
figure;
plot(t, x, 'b-', t, filtered_signal_fir, 'r-');
legend('原始信号', 'FIR低通滤波后');
title('FIR低通滤波结果');
xlabel('时间 (s)');
ylabel('幅度');
高通滤波器设计与应用
高通滤波器允许高于截止频率的信号通过，衰减低于截止频率的信号。这对于去除直流偏置或低频干扰很有用。
% 设计Butterworth高通滤波器
cutoff_freq = 200; % 截止频率200Hz
order = 6; % 滤波器阶数
Wn = cutoff_freq/(fs/2); % 归一化截止频率
[b, a] = butter(order, Wn, 'high');
% 应用滤波器
filtered_signal = filter(b, a, x_noisy);
% 绘制频率响应
figure;
freqz(b, a, 1024, fs);
title('Butterworth高通滤波器频率响应');
% 绘制时域结果
figure;
subplot(2,1,1);
plot(t, x_noisy);
title('含噪声信号');
xlabel('时间 (s)');
ylabel('幅度');
subplot(2,1,2);
plot(t, filtered_signal);
title('高通滤波后信号');
xlabel('时间 (s)');
ylabel('幅度');
在这个高通滤波器例子中，我们设置截止频率为200Hz，这意味着它应该显著衰减50Hz和120Hz的成分，而主要保留350Hz的成分。
带通滤波器设计与应用
带通滤波器只允许特定频带内的信号通过，这对于提取特定频率范围内的信息很有用。
% 设计Butterworth带通滤波器
low_cutoff = 80; % 低截止频率80Hz
high_cutoff = 150; % 高截止频率150Hz
order = 6; % 滤波器阶数
Wn = [low_cutoff high_cutoff]/(fs/2); % 归一化截止频率范围
[b, a] = butter(order, Wn, 'bandpass');
% 应用滤波器
filtered_signal = filter(b, a, x_noisy);
% 绘制频率响应
figure;
freqz(b, a, 1024, fs);
title('Butterworth带通滤波器频率响应');
% 绘制时域结果
figure;
subplot(2,1,1);
plot(t, x_noisy);
title('含噪声信号');
xlabel('时间 (s)');
ylabel('幅度');
subplot(2,1,2);
plot(t, filtered_signal);
title('带通滤波后信号');
xlabel('时间 (s)');
ylabel('幅度');
这个带通滤波器应该主要保留120Hz附近的频率成分，而衰减50Hz和350Hz的成分。
更复杂的滤波器设计
在实际应用中，我们可能需要更精确地控制滤波器的特性。MATLAB提供了多种高级滤波器设计方法。
使用切比雪夫滤波器
切比雪夫滤波器在通带或阻带有波纹，但可以提供更陡峭的滚降特性。
% 设计切比雪夫I型带通滤波器
passband_ripple = 1; % 通带波纹，单位dB
[b, a] = cheby1(order, passband_ripple, Wn, 'bandpass');
% 应用滤波器
filtered_signal_cheby = filter(b, a, x_noisy);
% 绘制频率响应
figure;
freqz(b, a, 1024, fs);
title('切比雪夫I型带通滤波器频率响应');
% 绘制时域结果
figure;
plot(t, filtered_signal, 'b-', t, filtered_signal_cheby, 'r-');
legend('Butterworth', '切比雪夫I型');
title('不同带通滤波器结果比较');
xlabel('时间 (s)');
ylabel('幅度');
使用椭圆滤波器
椭圆滤波器在通带和阻带都有波纹，但能提供最陡峭的滚降。
% 设计椭圆滤波器
passband_ripple = 1; % 通带波纹，dB
stopband_atten = 60; % 阻带衰减，dB
[b, a] = ellip(order, passband_ripple, stopband_atten, Wn, 'low');
% 应用滤波器
filtered_signal_ellip = filter(b, a, x_noisy);
% 绘制频率响应
figure;
freqz(b, a, 1024, fs);
title('椭圆低通滤波器频率响应');
% 绘制时域结果对比
figure;
subplot(2,1,1);
plot(t, filtered_signal); % Butterworth结果
title('Butterworth低通滤波结果');
xlabel('时间 (s)');
ylabel('幅度');
subplot(2,1,2);
plot(t, filtered_signal_ellip);
title('椭圆低通滤波结果');
xlabel('时间 (s)');
ylabel('幅度');
滤波器性能评估
设计滤波器后，评估其性能是很重要的。我们可以从时域和频域两个角度来评估。
频域评估
% 计算原始信号和滤波后信号的频谱
N = length(x);
X = fft(x_noisy, N);
X_filtered = fft(filtered_signal, N);
% 计算频率向量
f = (0:N-1)*(fs/N);
f = f(1:N/2+1); % 只取正频率部分
% 绘制频谱
figure;
subplot(2,1,1);
plot(f, abs(X(1:N/2+1)));
title('含噪声信号频谱');
xlabel('频率 (Hz)');
ylabel('幅度');
subplot(2,1,2);
plot(f, abs(X_filtered(1:N/2+1)));
title('滤波后信号频谱');
xlabel('频率 (Hz)');
ylabel('幅度');
时域评估
% 计算滤波前后的信噪比
signal_power = sum(x.^2)/length(x);
noise_power_before = sum((x_noisy – x).^2)/length(x);
noise_power_after = sum((filtered_signal – x).^2)/length(x);
SNR_before = 10*log10(signal_power/noise_power_before);
SNR_after = 10*log10(signal_power/noise_power_after);
fprintf('滤波前信噪比: %.2f dB\\n', SNR_before);
fprintf('滤波后信噪比: %.2f dB\\n', SNR_after);
fprintf('信噪比提升: %.2f dB\\n', SNR_after – SNR_before);
零相位滤波
在某些应用中，滤波引入的相位延迟可能是不可接受的。MATLAB的filtfilt函数可以实现零相位滤波。
% 使用filtfilt进行零相位滤波
zero_phase_filtered = filtfilt(b, a, x_noisy);
% 比较常规滤波和零相位滤波
figure;
subplot(3,1,1);
plot(t, x);
title('原始无噪声信号');
xlabel('时间 (s)');
ylabel('幅度');
subplot(3,1,2);
plot(t, filtered_signal);
title('常规滤波 (有相位延迟)');
xlabel('时间 (s)');
ylabel('幅度');
subplot(3,1,3);
plot(t, zero_phase_filtered);
title('零相位滤波 (无相位延迟)');
xlabel('时间 (s)');
ylabel('幅度');
实际应用案例：心电图信号处理
让我们看一个更实际的例子：处理心电图(ECG)信号。ECG信号通常受到多种噪声的干扰，包括电源干扰(50/60Hz)、肌电干扰和基线漂移。
% 模拟ECG信号（简化版）
fs_ecg = 500; % 采样频率500Hz
t_ecg = 0:1/fs_ecg:10-1/fs_ecg; % 10秒记录
ecg_freq = 1.2; % 心跳频率，约72次/分钟
% 创建QRS复合波（简化为高斯脉冲）
qrs_width = 0.08; % QRS宽度，秒
qrs_amp = 1; % QRS振幅
% 生成基本ECG
ecg_clean = zeros(size(t_ecg));
beats = round(t_ecg(end) * ecg_freq);
for i = 1:beats
beat_center = i/ecg_freq;
beat_idx = abs(t_ecg – beat_center) < qrs_width;
ecg_clean(beat_idx) = qrs_amp * exp(-(t_ecg(beat_idx) – beat_center).^2/(qrs_width/4)^2);
end
% 添加噪声
powerline_noise = 0.2 * sin(2*pi*50*t_ecg); % 50Hz电源噪声
baseline_drift = 0.3 * sin(2*pi*0.3*t_ecg); % 基线漂移
random_noise = 0.1 * randn(size(t_ecg)); % 随机噪声
ecg_noisy = ecg_clean + powerline_noise + baseline_drift + random_noise;
% 绘制原始和含噪声信号
figure;
subplot(2,1,1);
plot(t_ecg, ecg_clean);
title('原始ECG信号');
xlabel('时间 (s)');
ylabel('振幅');
xlim([0 3]); % 只显示前3秒
subplot(2,1,2);
plot(t_ecg, ecg_noisy);
title('含噪声ECG信号');
xlabel('时间 (s)');
ylabel('振幅');
xlim([0 3]); % 只显示前3秒
% 设计带通滤波器处理ECG (常用范围：0.5-40Hz)
low_cut = 0.5;
high_cut = 40;
order = 4;
Wn_ecg = [low_cut high_cut]/(fs_ecg/2);
[b_ecg, a_ecg] = butter(order, Wn_ecg, 'bandpass');
% 使用零相位滤波
ecg_filtered = filtfilt(b_ecg, a_ecg, ecg_noisy);
% 绘制结果
figure;
subplot(3,1,1);
plot(t_ecg, ecg_clean);
title('原始ECG信号');
xlabel('时间 (s)');
ylabel('振幅');
xlim([0 3]);
subplot(3,1,2);
plot(t_ecg, ecg_noisy);
title('含噪声ECG信号');
xlabel('时间 (s)');
ylabel('振幅');
xlim([0 3]);
subplot(3,1,3);
plot(t_ecg, ecg_filtered);
title('滤波后ECG信号');
xlabel('时间 (s)');
ylabel('振幅');
xlim([0 3]);
总结与最佳实践
在MATLAB中使用滤波器时，可以记住以下几点最佳实践：
MATLAB的Signal Processing Toolbox提供了强大的工具集，使得滤波器设计和应用变得相对简单。通过理解不同滤波器的特性和适用场景，你可以更有效地处理各种实际信号处理问题。
希望这个教程能帮助你掌握MATLAB中滤波器的基本使用方法！通过实践和进一步探索，你将能够应对更复杂的信号处理挑战。
