1. 问题描述
本人并非信号处理专业，仅在结构监测研究中遇到滤波问题，特总结常规的低通滤波技术，去除高频噪音。
由于环境的干扰因素，监测信号中总会包含噪音成分，影响信号处理过程，如下图：
接收信号中出现很多“毛刺”，即为高频噪音，预期通过低通滤波器过滤处理。
2. 技术背景
在MATLAB中有很多种滤波器可供选择，本文仅介绍一笔者实现的滤波方式：切比雪夫滤波器。
低通滤波的技术要点有：
- 滤波器参数设置
[n,Wp]=cheb1ord(Wp,Ws,Rp,Rs); % Cheby1
[b,a]=cheby1(n,Rp,Wp);
freqz(b,a,2048,fs); % 查看设计滤波器的曲线
- 信号滤波运算
y = filter(b,a,x);
此处仅说明代码实现，理论问题不再说明。
3. 解决方案
滤波器参数的设置是有效滤波的关键，最重要的参数是确定滤波的范围：
- 通过频率$f_{pass}$
- 截止频率$f_{stop}$
上图可以看出，原信号的频域范围主要在100~300kHz。故可以设置：
- 通过频率$f_{pass}= 300 kHz$
- 截止频率$f_{stop}= 500 kHz$
即过滤掉500 kHz以上的高频噪音。
4. 实施示例
4.1 数据读入
%% 数据读入
clc,clear,close all
[M,dt] = tools.getcsv(); % 读入csv信号和采样周期dt
fs = 1/dt; % 采样频率
t = M(:,1);
s = detrend(M(:,3)); % 去趋势的信号
4.2 滤波参数设置
%% 参数设置
prompt0 = { % 对话框参数
'通过频率 f-pass(kHz)', 300
'截止频率 f-stop(kHz)', 500
'Passband ripple in decibels Rp',0.1
'衰减值Rs(Db)',30
};
dlg0.title = '滤波参数输入-马骋';
dlg0.save = 'lp';
para_input = tools.paradlg(prompt0,dlg0);
para.f1 = para_input{1}*1e3;
para.f3 = para_input{2}*1e3;
para.rp = para_input{3};
para.rs = para_input{4};
para.fs = fs;
注：以上tools为笔者自定义函数工具箱。
4.3 滤波器生成
%% cheby1低通滤波图示
para.type = 1; % 滤波器类型：切比雪夫-1
s_lp = tools.lowp(s,para); % 滤波
可以看出，滤波器在频域300-500 kHz范围内逐渐衰减。
4.4 滤波效果
%% 处理信号绘图
figure
plot([t t],[s s_lp])
legend({'原始信号','低通滤波信号'})
title('cheby1低通滤波效果示例'),grid on
xlim([min(t) max(t)])
figure
subplot(211)
plot(t,s)
legend('原始信号'),grid on
xlim([min(t) max(t)])
subplot(212)
plot(t,s_lp)
legend('滤波信号'),grid on
xlim([min(t) max(t)])
显然，滤波后的信号平滑很多。
5. 常见问题
滤波核心函数如下：
function y=lowp(x,para)
% 题目: 低通滤波器
% 输入:
% x -- 原始信号序列
% para.
% f1 -- 通带截止频率
% f3 -- 阻带截止频率
% rp -- 边带区衰减DB数设置
% rs -- 截止区衰减DB数设置
% fs -- 序列x的采样频率
% type-- 滤波器类型
% 输出：
% y -- 滤波后的信号
% 功能：
% 低通滤波，滤除高频噪音
% Cheby1
% Butterworth
% 注意：
% 通带或阻带的截止频率的选取范围是不能超过采样率的一半
% f1,f3的值都要小于fs/2
% rp=0.1;rs=30;%通带边衰减DB值和阻带边衰减DB值
% 作者: 未知
% 修改: 马骋
% 2016.04.21 @HIT
%% 参数输入
f1 = para.f1;
f3 = para.f3;
Rp = para.rp;
Rs = para.rs;
fs = para.fs;
%% 滤波器设计
Wp = f1/(fs/2); % 采用fs/2归一化,Nyquist frequency.
Ws = f3/(fs/2);
if para.type==1
[n,Wp]=cheb1ord(Wp,Ws,Rp,Rs); % Cheby1
[b,a]=cheby1(n,Rp,Wp);
freqz(b,a,2048,fs); % 查看设计滤波器的曲线
title(sprintf('n = %d Cheby1 Lowpass Filter',n))
xlim([0 f3])
else
[n,Wn] = buttord(Wp,Ws,Rp,Rs,'s'); % Butterworth
[b,a] = butter(n,Wn,'s'); % 计算滤波器系统函数分子分母多项式
[z,p,k] = butter(n,Wn);
sos = zp2sos(z,p,k);
freqz(sos,2048,fs)
title(sprintf('n = %d Butterworth Lowpass Filter',n))
xlim([0 f3])
end
%% 滤波
y = filter(b,a,x); % 对序列x滤波后得到的序列y
end % lowp
注：此函数中，仅切比雪夫-1滤波器测试成功，2型滤波器测试失败。
示例程序下载：
https://coding.net/u/frank0449/p/MATLAB_lowpassFilter/git
本文用时 25 m
