别再只会用公式了！手把手教你用MATLAB实现一阶数字低通滤波器（附完整代码）

别再只会用公式了！手把手教你用MATLAB实现一阶数字低通滤波器（附完整代码）

在信号处理领域，滤波技术就像是一位隐形的"清洁工"，它能帮我们去除信号中的噪声和干扰，保留真正有价值的信息。对于工程师和科研人员来说，掌握滤波器的实际应用远比理解公式推导更为重要。本文将带你从零开始，用MATLAB实现一个实用的一阶数字低通滤波器，让你在处理传感器数据时游刃有余。

1. 理解一阶数字低通滤波器的核心

一阶数字低通滤波器之所以广受欢迎，是因为它在简单性和有效性之间取得了完美平衡。想象一下，当你用陀螺仪测量物体姿态时，原始信号往往包含高频噪声；或者用麦克风录音时，环境中的杂音干扰了主要声音。这时，一个设计得当的低通滤波器就能成为你的得力助手。

滤波器的核心差分方程看似简单：

y(n) = a*x(n) + (1-a)*y(n-1)

但这个方程中蕴含着几个关键点：

• x(n)：当前时刻的原始输入信号
• y(n)：当前时刻的滤波后输出
• y(n-1)：上一时刻的滤波结果
• 参数a：决定滤波器特性的关键因子

参数a的计算公式为：

a = (2πfcTs)/(1+2πfcTs)

其中：

• fc：截止频率(Hz)
• Ts：采样周期(1/fs)
• fs：采样频率(Hz)

注意：截止频率fc应该远小于采样频率fs，通常建议fs至少是fc的5-10倍，否则滤波效果会大打折扣。

2. MATLAB实现基础版本

让我们从最基本的实现开始。假设我们有一个被噪声污染的传感器信号，采样频率为1000Hz，我们希望滤除50Hz以上的噪声。

% 基本参数设置 fs = 1000; % 采样频率(Hz) fc = 50; % 截止频率(Hz) Ts = 1/fs; % 采样周期(s) wc = 2*pi*fc; % 截止角频率(rad/s) % 计算滤波系数a a = (wc*Ts)/(1+wc*Ts); % 生成测试信号：1Hz正弦波+高频噪声 t = 0:Ts:1; % 时间向量(1秒时长) signal = sin(2*pi*1*t) + 0.5*randn(size(t)); % 信号+噪声 % 初始化滤波后的信号 filtered_signal = zeros(size(signal)); filtered_signal(1) = signal(1); % 初始值设为第一个采样点 % 执行滤波(使用for循环实现) for n = 2:length(signal) filtered_signal(n) = a*signal(n) + (1-a)*filtered_signal(n-1); end % 绘制结果 figure; subplot(2,1,1); plot(t, signal); title('原始信号(含噪声)'); xlabel('时间(s)'); ylabel('幅值'); subplot(2,1,2); plot(t, filtered_signal); title('滤波后信号'); xlabel('时间(s)'); ylabel('幅值');

这段代码展示了最基本的实现方式，但实际应用中我们还需要考虑更多因素。

3. 高级实现与优化技巧

基础版本虽然简单，但在处理大数据量时效率不高。MATLAB的优势在于向量化运算，我们可以利用这一特性来优化性能。

3.1 向量化实现

function filtered_signal = vectorized_lowpass(signal, fc, fs) % 参数计算 wc = 2*pi*fc; Ts = 1/fs; a = (wc*Ts)/(1+wc*Ts); % 初始化输出 filtered_signal = zeros(size(signal)); filtered_signal(1) = signal(1); % 向量化滤波 for n = 2:length(signal) filtered_signal(n) = a*signal(n) + (1-a)*filtered_signal(n-1); end end

虽然这个版本看起来和之前类似，但我们可以进一步利用MATLAB的filter函数实现真正的向量化：

function filtered_signal = builtin_lowpass(signal, fc, fs) % 计算滤波器系数 wc = 2*pi*fc; Ts = 1/fs; a = (wc*Ts)/(1+wc*Ts); % 使用MATLAB内置filter函数 b = [a]; % 分子系数 a_filter = [1, -(1-a)]; % 分母系数 filtered_signal = filter(b, a_filter, signal); end

3.2 实时处理实现

在实际应用中，我们经常需要实时处理数据流。这时，我们需要维护滤波器的状态：

classdef RealTimeLowPassFilter < handle properties a % 滤波系数 last_output % 上一次的输出值 is_initialized % 是否已初始化 end methods function obj = RealTimeLowPassFilter(fc, fs) % 构造函数 wc = 2*pi*fc; Ts = 1/fs; obj.a = (wc*Ts)/(1+wc*Ts); obj.is_initialized = false; end function output = process(obj, input) % 处理单个输入样本 if ~obj.is_initialized obj.last_output = input; obj.is_initialized = true; end output = obj.a*input + (1-obj.a)*obj.last_output; obj.last_output = output; end function reset(obj) % 重置滤波器状态 obj.is_initialized = false; end end end

使用示例：

% 创建滤波器实例 filter = RealTimeLowPassFilter(50, 1000); % 模拟实时处理 for i = 1:length(raw_data) filtered_value = filter.process(raw_data(i)); % 使用filtered_value... end

4. 参数选择与性能分析

选择合适的滤波器参数至关重要。让我们通过实验来分析不同参数的影响。

4.1 截止频率的影响

我们固定采样频率为1000Hz，观察不同截止频率的效果：

% 测试不同截止频率 fc_values = [10, 50, 100, 200]; fs = 1000; t = 0:1/fs:1; signal = sin(2*pi*5*t) + 0.3*randn(size(t)); figure; for i = 1:length(fc_values) fc = fc_values(i); filtered = builtin_lowpass(signal, fc, fs); subplot(length(fc_values),1,i); plot(t, signal, 'b', t, filtered, 'r', 'LineWidth', 1.5); title(sprintf('截止频率 = %d Hz', fc)); legend('原始信号', '滤波后'); end

4.2 采样频率的影响

截止频率固定为50Hz，改变采样频率：

fc = 50; fs_values = [100, 200, 500, 1000]; figure; for i = 1:length(fs_values) fs = fs_values(i); Ts = 1/fs; t = 0:Ts:1; signal = sin(2*pi*5*t) + 0.3*randn(size(t)); % 确保有足够的数据点 if fs < 2*fc warning('采样频率低于奈奎斯特频率，可能出现混叠'); end filtered = builtin_lowpass(signal, fc, fs); subplot(length(fs_values),1,i); plot(t, signal, 'b', t, filtered, 'r', 'LineWidth', 1.5); title(sprintf('采样频率 = %d Hz (fc=%dHz)', fs, fc)); legend('原始信号', '滤波后'); end

重要提示：采样频率fs必须至少是截止频率fc的2倍（奈奎斯特准则），实际应用中建议保持fs ≥ 5fc以获得良好效果。

5. 实际应用案例分析

让我们看几个实际工程中的应用场景。

5.1 陀螺仪数据平滑

假设我们从IMU获取的陀螺仪数据存在高频噪声：

% 模拟陀螺仪数据 fs = 200; % 采样频率200Hz t = 0:1/fs:10; % 10秒数据 true_velocity = 50 + 10*sin(2*pi*0.2*t); % 真实角速度(°/s) noisy_velocity = true_velocity + 5*randn(size(t)); % 添加噪声 % 设计滤波器(fc=5Hz) fc = 5; filtered_velocity = builtin_lowpass(noisy_velocity, fc, fs); % 绘制结果 figure; plot(t, noisy_velocity, 'b', t, filtered_velocity, 'r', t, true_velocity, 'g--', 'LineWidth', 1.5); legend('含噪声数据', '滤波后', '真实值'); xlabel('时间(s)'); ylabel('角速度(°/s)'); title('陀螺仪数据滤波效果'); grid on;

5.2 音频信号处理

处理含有高频噪声的音频信号：

% 读取音频文件 [audio, fs] = audioread('noisy_audio.wav'); % 假设已有含噪声的音频文件 % 设计滤波器(fc=4kHz用于语音) fc = 4000; filtered_audio = builtin_lowpass(audio, fc, fs); % 保存结果 audiowrite('filtered_audio.wav', filtered_audio, fs); % 绘制频谱对比 n = length(audio); f = (0:n-1)*(fs/n); audio_fft = abs(fft(audio)); filtered_fft = abs(fft(filtered_audio)); figure; subplot(2,1,1); plot(f(1:n/2), audio_fft(1:n/2)); title('原始音频频谱'); xlabel('频率(Hz)'); subplot(2,1,2); plot(f(1:n/2), filtered_fft(1:n/2)); title('滤波后音频频谱'); xlabel('频率(Hz)');

5.3 传感器数据实时处理

使用前面创建的RealTimeLowPassFilter类进行实时处理：

% 模拟实时数据流 fs = 100; % 采样频率100Hz duration = 10; % 10秒 num_samples = fs * duration; % 创建滤波器(fc=2Hz) filter = RealTimeLowPassFilter(2, fs); % 初始化存储 filtered_data = zeros(1, num_samples); raw_data = sin(2*pi*0.5*(0:num_samples-1)/fs) + 0.2*randn(1, num_samples); % 实时处理循环 for i = 1:num_samples filtered_data(i) = filter.process(raw_data(i)); % 这里可以添加其他实时处理逻辑 % 例如: if mod(i, fs) == 0 % 每秒执行一次操作 end % 绘制结果 t = (0:num_samples-1)/fs; figure; plot(t, raw_data, 'b', t, filtered_data, 'r', 'LineWidth', 1.5); legend('原始数据', '滤波后'); xlabel('时间(s)'); title('实时滤波效果');

6. 常见问题与调试技巧

在实际应用中，你可能会遇到以下问题：

6.1 滤波器响应太慢

症状：滤波后的信号明显滞后于原始信号，变化迟缓。

可能原因：

• 截止频率设置过低
• 采样频率与截止频率比例不当

解决方案：

1. 适当提高截止频率fc
2. 检查fs/fc比值，确保至少为5:1
3. 使用以下代码测试不同参数：

fc_test_values = linspace(1, 20, 5); % 测试1Hz到20Hz for fc_test = fc_test_values % 测试代码... end

6.2 滤波效果不明显

症状：滤波前后信号几乎没有变化。

可能原因：

• 截止频率设置过高
• 噪声频率与信号频率重叠
• 采样频率不足

解决方案：

1. 先分析信号频谱，确定噪声频率：

[pxx, f] = pwelch(signal, [], [], [], fs); figure; plot(f, 10*log10(pxx)); xlabel('频率(Hz)'); ylabel('功率谱密度(dB/Hz)');

2. 根据频谱分析结果调整截止频率
3. 考虑使用更高阶滤波器或其它滤波技术

6.3 数值不稳定

症状：滤波后的信号出现NaN或异常值。

可能原因：

• 采样频率过低导致a值接近1
• 数值累积误差

解决方案：

1. 检查a值计算：

a = (wc*Ts)/(1+wc*Ts); disp(['a值为: ', num2str(a)]); % a值应在0.01到0.3之间较为理想

2. 使用双精度计算
3. 定期重置滤波器状态（对实时滤波器）

6.4 相位失真问题

一阶滤波器会导致相位延迟，这在某些应用中可能不可接受。如果需要零相位延迟，可以考虑使用双向滤波：

function zero_phase_signal = zero_phase_lowpass(signal, fc, fs) % 正向滤波 forward_filtered = builtin_lowpass(signal, fc, fs); % 反向滤波 reversed_signal = forward_filtered(end:-1:1); backward_filtered = builtin_lowpass(reversed_signal, fc, fs); % 再次反向 zero_phase_signal = backward_filtered(end:-1:1); end

注意：双向滤波会引入处理延迟，只适用于非实时应用。

7. 性能优化与高级话题

当处理大规模数据或对性能有严格要求时，可以考虑以下优化策略：

7.1 使用C-MEX加速

对于MATLAB，可以编写C-MEX函数来加速滤波计算：

// lowpass_filter.c - MEX函数实现 #include "mex.h" void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { double *signal, *filtered_signal, a; mwSize n, i; // 检查输入输出参数 if (nrhs != 2 || nlhs != 1) mexErrMsgTxt("用法: filtered = lowpass_filter(signal, a)"); // 获取输入 signal = mxGetPr(prhs[0]); a = mxGetScalar(prhs[1]); n = mxGetNumberOfElements(prhs[0]); // 创建输出 plhs[0] = mxCreateDoubleMatrix(1, n, mxREAL); filtered_signal = mxGetPr(plhs[0]); // 执行滤波 filtered_signal[0] = signal[0]; for (i = 1; i < n; i++) filtered_signal[i] = a*signal[i] + (1-a)*filtered_signal[i-1]; }

编译和使用：

mex lowpass_filter.c % 编译 filtered = lowpass_filter(signal, a); % 使用

7.2 多通道并行处理

如果需要同时处理多个信号通道，可以利用MATLAB的矩阵运算：

function filtered_data = multi_channel_lowpass(data, fc, fs) % data: 每列代表一个通道的信号 [num_samples, num_channels] = size(data); % 计算系数 wc = 2*pi*fc; Ts = 1/fs; a = (wc*Ts)/(1+wc*Ts); % 初始化输出 filtered_data = zeros(size(data)); filtered_data(1,:) = data(1,:); % 滤波处理 for n = 2:num_samples filtered_data(n,:) = a*data(n,:) + (1-a)*filtered_data(n-1,:); end end

7.3 与其它滤波技术结合

一阶滤波器可以与其他滤波技术结合使用：

% 先使用移动平均滤波预处理 window_size = 5; smoothed = movmean(noisy_signal, window_size); % 再应用一阶低通滤波 fc = 10; fs = 1000; final_filtered = builtin_lowpass(smoothed, fc, fs);

这种组合方式可以在保持计算效率的同时获得更好的滤波效果。

8. 完整工具箱实现

为了方便日常使用，我们可以将上述功能封装成一个完整的工具箱：

classdef LowPassFilterToolbox methods (Static) function a = calculate_coefficient(fc, fs) % 计算滤波系数a wc = 2*pi*fc; Ts = 1/fs; a = (wc*Ts)/(1+wc*Ts); end function filtered = filter_signal(signal, fc, fs) % 基本滤波函数 a = LowPassFilterToolbox.calculate_coefficient(fc, fs); filtered = zeros(size(signal)); filtered(1) = signal(1); for n = 2:length(signal) filtered(n) = a*signal(n) + (1-a)*filtered(n-1); end end function filtered = zero_phase_filter(signal, fc, fs) % 零相位滤波 forward = LowPassFilterToolbox.filter_signal(signal, fc, fs); backward = LowPassFilterToolbox.filter_signal(forward(end:-1:1), fc, fs); filtered = backward(end:-1:1); end function analyze_response(fc, fs, varargin) % 分析滤波器频率响应 a = LowPassFilterToolbox.calculate_coefficient(fc, fs); % 计算频率响应 [h, w] = freqz([a], [1, -(1-a)], 1024, fs); % 绘制幅频响应 figure; subplot(2,1,1); plot(w, 20*log10(abs(h))); title('幅频响应'); xlabel('频率(Hz)'); ylabel('增益(dB)'); grid on; % 绘制相频响应 subplot(2,1,2); plot(w, unwrap(angle(h))*180/pi); title('相频响应'); xlabel('频率(Hz)'); ylabel('相位(度)'); grid on; end function compare_with_builtin(fc, fs) % 与MATLAB内置滤波器比较 a = LowPassFilterToolbox.calculate_coefficient(fc, fs); % 设计数字滤波器 [b, a_filter] = butter(1, fc/(fs/2)); % 比较频率响应 figure; [h1, w] = freqz([a], [1, -(1-a)], 1024, fs); [h2, ~] = freqz(b, a_filter, 1024, fs); semilogx(w, 20*log10(abs([h1, h2]))); legend('一阶滤波器', 'Butterworth一阶'); title('频率响应比较'); xlabel('频率(Hz)'); ylabel('增益(dB)'); grid on; end end end

使用示例：

% 使用工具箱 signal = randn(1, 1000); % 测试信号 fc = 50; fs = 1000; % 基本滤波 filtered = LowPassFilterToolbox.filter_signal(signal, fc, fs); % 零相位滤波 zero_phase = LowPassFilterToolbox.zero_phase_filter(signal, fc, fs); % 分析响应 LowPassFilterToolbox.analyze_response(fc, fs); % 与内置滤波器比较 LowPassFilterToolbox.compare_with_builtin(fc, fs);