1. 問題描述
本人并非信號(hào)處理專業(yè)，僅在結(jié)構(gòu)監(jiān)測(cè)研究中遇到濾波問題，特總結(jié)常規(guī)的低通濾波技術(shù)，去除高頻噪音。
由于環(huán)境的干擾因素，監(jiān)測(cè)信號(hào)中總會(huì)包含噪音成分，影響信號(hào)處理過程，如下圖：
接收信號(hào)中出現(xiàn)很多“毛刺”，即為高頻噪音，預(yù)期通過低通濾波器過濾處理。
2. 技術(shù)背景
在MATLAB中有很多種濾波器可供選擇，本文僅介紹一筆者實(shí)現(xiàn)的濾波方式：切比雪夫?yàn)V波器。
低通濾波的技術(shù)要點(diǎn)有：
- 濾波器參數(shù)設(shè)置
[n,Wp]=cheb1ord(Wp,Ws,Rp,Rs); % Cheby1
[b,a]=cheby1(n,Rp,Wp);
freqz(b,a,2048,fs); % 查看設(shè)計(jì)濾波器的曲線
- 信號(hào)濾波運(yùn)算
y = filter(b,a,x);
此處僅說明代碼實(shí)現(xiàn)，理論問題不再說明。
3. 解決方案
濾波器參數(shù)的設(shè)置是有效濾波的關(guān)鍵，最重要的參數(shù)是確定濾波的范圍：
- 通過頻率$f_{pass}$
- 截止頻率$f_{stop}$
上圖可以看出，原信號(hào)的頻域范圍主要在100~300kHz。故可以設(shè)置：
- 通過頻率$f_{pass}= 300 kHz$
- 截止頻率$f_{stop}= 500 kHz$
即過濾掉500 kHz以上的高頻噪音。
4. 實(shí)施示例
4.1 數(shù)據(jù)讀入
%% 數(shù)據(jù)讀入
clc,clear,close all
[M,dt] = tools.getcsv(); % 讀入csv信號(hào)和采樣周期dt
fs = 1/dt; % 采樣頻率
t = M(:,1);
s = detrend(M(:,3)); % 去趨勢(shì)的信號(hào)
4.2 濾波參數(shù)設(shè)置
%% 參數(shù)設(shè)置
prompt0 = { % 對(duì)話框參數(shù)
'通過頻率 f-pass(kHz)', 300
'截止頻率 f-stop(kHz)', 500
'Passband ripple in decibels Rp',0.1
'衰減值Rs(Db)',30
};
dlg0.title = '濾波參數(shù)輸入-馬騁';
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為筆者自定義函數(shù)工具箱。
4.3 濾波器生成
%% cheby1低通濾波圖示
para.type = 1; % 濾波器類型：切比雪夫-1
s_lp = tools.lowp(s,para); % 濾波
可以看出，濾波器在頻域300-500 kHz范圍內(nèi)逐漸衰減。
4.4 濾波效果
%% 處理信號(hào)繪圖
figure
plot([t t],[s s_lp])
legend({'原始信號(hào)','低通濾波信號(hào)'})
title('cheby1低通濾波效果示例'),grid on
xlim([min(t) max(t)])
figure
subplot(211)
plot(t,s)
legend('原始信號(hào)'),grid on
xlim([min(t) max(t)])
subplot(212)
plot(t,s_lp)
legend('濾波信號(hào)'),grid on
xlim([min(t) max(t)])
顯然，濾波后的信號(hào)平滑很多。
5. 常見問題
濾波核心函數(shù)如下：
function y=lowp(x,para)
% 題目: 低通濾波器
% 輸入:
% x -- 原始信號(hào)序列
% para.
% f1 -- 通帶截止頻率
% f3 -- 阻帶截止頻率
% rp -- 邊帶區(qū)衰減DB數(shù)設(shè)置
% rs -- 截止區(qū)衰減DB數(shù)設(shè)置
% fs -- 序列x的采樣頻率
% type-- 濾波器類型
% 輸出：
% y -- 濾波后的信號(hào)
% 功能：
% 低通濾波，濾除高頻噪音
% Cheby1
% Butterworth
% 注意：
% 通帶或阻帶的截止頻率的選取范圍是不能超過采樣率的一半
% f1,f3的值都要小于fs/2
% rp=0.1;rs=30;%通帶邊衰減DB值和阻帶邊衰減DB值
% 作者: 未知
% 修改: 馬騁
% 2016.04.21 @HIT
%% 參數(shù)輸入
f1 = para.f1;
f3 = para.f3;
Rp = para.rp;
Rs = para.rs;
fs = para.fs;
%% 濾波器設(shè)計(jì)
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); % 查看設(shè)計(jì)濾波器的曲線
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'); % 計(jì)算濾波器系統(tǒng)函數(shù)分子分母多項(xiàng)式
[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); % 對(duì)序列x濾波后得到的序列y
end % lowp
注：此函數(shù)中，僅切比雪夫-1濾波器測(cè)試成功，2型濾波器測(cè)試失敗。
示例程序下載：
https://coding.net/u/frank0449/p/MATLAB_lowpassFilter/git
本文用時(shí) 25 m
