在matlab里设计使用Bessel Filter（贝塞尔滤波器）进行数据降噪的简单例子

作为一个土木工程系的研究生，现在大家做数据处理也用得越来越多了，自己接触的东西也多。作为土木工程师，很大部分时间其实并不想弄清楚滤波器的原理，只是希望拿来用，拿来把自己采集到的数据降噪一下，拿来做个傅立叶分析看看结构的固有频率变没变。但是matlab自带的文档很多时候都默认大家知道滤波器的原理，很多专有名词搞得不太友好，所以在这里以问题解答的形式直接给出我在科研中遇到的一个简单的用贝塞尔滤波器的例子，拿来主义万岁。

背景知识：采集到的数据都是discrete（离散）的数据，采集数据的频率为sampling rate（采样频率）。贝塞尔滤波器为analog filter（模拟滤波器），需要转换为digital filter（数字滤波器）才能直接使用。总体

首先题设是：

已知以1000 Hz采样频率得到的数据x，里面会有若干噪声，求使用贝塞尔滤波器滤掉噪声。
这里令采集到的数据为

    Fs = 1000;
    t  = linspace(0,1,Fs);
    x  = cos(2*pi*(round(rand*100))*t)+0.5*randn(size(t));

可以看到这里结构固有频率应为某个特定值round(rand*100)，信号里有均匀分布噪声。

解题思路：

第一步 - 快速傅立叶分析决定 cut-off frequency（截止频率），也就是确定我们要滤掉多少Hz以上的数据

第二步 - 使用截止频率来建立贝塞尔滤波器的模拟形式（analog filter）

第三步 - 将模拟形式转换为数字形式（digital filter）

第四步 - 应用滤波器

第五步 - 检查结果

---

解：

第一步：
用FFT分析原始数据，得到原始数据的频率分布如下图

我们可以看到该段数据的主要频率集中在82 Hz附近，所以我们可以设cut-off frequency在128 Hz以上。

第二步：
利用matlab里的自带函数besself建立贝塞尔滤波器,截止频率128 Hz
建立的滤波器需要另一个参数n，表示滤波器的阶数，越高则滤波器切断越明显，一般5阶左右，可以自由尝试。

    [b, a] = besself(5, 128);

这样我们就得到了滤波器的参数b和a。

第三步：
利用自带函数bilinear转换模拟滤波器至数字滤波器

    [num, den] = bilinear(b, a, 128);

这样我们得到了滤波器参数num和den，分别是两个向量。里面的数字就是滤波器表达式里的分子与分母的行列式。

第四步：
使用滤波器完成滤波，除去高于128Hz部分

    x_filtered = filter(num, den, x);

得到的x_filtered就是滤波后结果。

第五步：
检查结果，对x_filtered使用FFT结果如下

可以从频域图中看出，经过滤波之后的信号，噪音明显变小，150Hz以上成分几乎消失。我们降噪的目的也达到了。

至于为什么截止频率设为128Hz但结果却是150Hz以上成分明显变小？答案是这个与贝塞尔滤波器的特性相关。贝塞尔滤波器对于截止频率附近的数据会渐近地减少，直到接近于零，所以是从128Hz开始渐渐减少后面频率成分的amplitude，直到将近200Hz可以完全将amplitude降到零。其滤波器的频率响应图如下（可以看成是将该图与原数据的频域图做点乘）

最终我们得到的降噪之后的数据为 x_filtered

重复本例时的完整代码：

    Fs = 1000;
    t = linspace(0,1,Fs);
    x = cos(2*pi*82*t)+0.5*randn(size(t)); % 做模拟时将82改为round(rand*100)
    [b,a] = besself(10,128);
    [num,den] = bilinear(b,a,128);
    x_filtered = filter(num,den,x);
    plot(x); 
    hold on;
    plot(x_filtered,'r');
    legend('original signal','denoised signal')