低通、⾼通数字滤波器——C语⾔单⽚机实现
博主刚好进⼊研⼆，研究的⽅向刚好涉及到数字滤波这⼀块，因此花了⼀周时间钻研了下数字滤波的实现。
由于本科是电⽓专业，所以没有数字信号处理相关知识，在⼀开始看数字信号处理相关理论的时候就显得⽐较⼒不从⼼，尤其是难懂的数学公式。相⽐看到这⾥的读者多多少少也有类似的体会。 好在功夫不负有⼼⼈，本博主从繁琐的公式中，加上其他博主的博客讲解，领悟了如何使⽤C代码实现⼏种经典数字滤波器，可以使其在VS 或者单⽚机上运⾏⽽不受限于在matlab上跑仿真。
当然在本篇博客中，主要是以巴特沃斯滤波器做讲解，其他经典滤波器由于数学公式和讲解都⽐较少，因此没研究出来。若有读者发现较好的资源，欢迎留⾔⼀起分享。
好的，读到这⾥，先讲解⼀个同学进场会遇到却不是太理解的概念——滤波器的阶数代表什么意思？其实博主⼀开始也不理解，后来⽤C代码实现的时候突然顿悟的。答案是：阶数N，代表着在该滤波器中，需要前N个时刻的采样点数据参与运算。
读到这⾥可能有的读者还不是太理解。没关系，当你读完整篇博客之后，再回过头都这句话试⼀下。（其实是博主整理了半天语⾔也没想到⼀句通俗易懂的话来描述，如果想到了后期再回来更新吧）
⼀阶滤波器
⼀阶滤波器呢，是使⽤数字滤波器来模拟硬件RC滤波器，以抑制⼲扰信号。⼀阶低通数字滤波器公式如下：
其中：
好，看到这⾥，我们可以看到整个公式中，唯⼀的需要读者⾃⼰计算的就是常数q，有了q值，整个数学公式，相⽐读者就可以使⽤C语⾔实现。那么常数q是如何根据我们所需的截⽌频率和采样间隔求出的呢？请看如下公式：
其中：
到这⾥，我们可以得到时间常数q的计算公式为：
下⾯举个例⼦教⼤家如何计算。
例如需要⼀个20Hz的低通滤波器，对信号的采样率为1000Hz（即0.001秒采样⼀次），因此q的值为206.280.001=0.125。将q值带⼊到数学公式中，以C语⾔实现即可。
⾼阶滤波器
讲解完了⼀阶滤波器，下⾯我们来看下⾼阶滤波器的实现。
看到这⾥，有读者可能会问，这个阶数怎么选，是选择使⽤⼀阶滤波器还是⾼阶的，⾼阶滤波器究竟选择⼏阶？这其实是有⼀个数学公式计算你所需的滤波器阶数的，但是由于太过繁琐，不好记，因此在⼯程领域中不推荐记忆。
这⾥给⼤家⼀个经验，就是滤波器阶数越⾼，滤波效果越好。
因此，若读者的处理器性能⾜够强，直接上⾼阶滤波器即可（阶数越⾼，计算量越⼤），⼀般四到六阶即可解决⼤部分问题。（有的场景需要⾼阶滤波器才会完全滤除杂波，但是受到处理器性能影响滤除不⼲净，因此可以适当降低阶数，因此使⽤公式计算所需阶数基本没有意义，处理器性能达不到，只能降阶层，因此⾃⼰⽤代码测⼀下波形⽐通过公式计算简单得多）
下⾯是IIR滤波器的差分⽅程（直接型）。
其中：
看到这⾥，可能有的读者已经觉得公式有点繁琐。不过没有关系，我们继续往下看，代码实现后在会过来品味⼀下这个公式。
下⾯是通过matlab⼯具计算两个系数。
⾸先在matlab命令⾏中输⼊ fdatool或者filterDesigner，回车执⾏后会弹出⼀个框，如下图所⽰。
根据业务需求，选择低通或者⾼通滤波器（注：低通、⾼通、带通滤波器均为同⼀个公式，系数不同⽽已），并设计阶数，输⼊采样频率和截⾄频率，点击设计滤波器。
观察幅频响应曲线可以看到，从20Hz开始衰减，50Hz是衰减达到50dB，此处需读者根据⾃⼰使⽤情况调节截⽌频率的参数。
如下图所⽰，点击两个按钮，将滤波器系数转换为单节结构，得到滤波器系数。分母为y(n)的系数ak，分⼦为x(n)系数bk。
将两组参数带⼊，得到下⾯公式。
那本⽂中，列举了低通滤波器的实现⽅法，⾼通滤波器和带通滤波器的实现⽅法同上，更改系数即可实现。
下⾯是代码展⽰：
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.Threading.Tasks;
using 肌电可视化界⾯.ViewModel;
using NWaves;
using NWaves.Signals;
using NWaves.Filters;
using NWaves.Filters.Bessel;
namespace 肌电可视化界⾯
{
class SingleProcess
{
int FilterDataNum_VAR = 10000; //全局变量--采集信号的长度
int FilterDataFrequence_VAR = 10000; //全局变量--采集信号的频率
int LowPassCutoffFre_VAR = 150; //全局变量--信号预处理低通截⾄频率
public Tuple<double[],double[]> Pretreatment()
{
double FilterSamplingInterval = 1.0 / Convert.ToDouble( FilterDataFrequence_VAR); //全局变量--采集信号的间隔
double[] Sin200Hz=new double[FilterDataNum_VAR],Sin300Hz = new double[FilterDataNum_VAR], Sin1Hz = new double[FilterDataNum_VAR], Sin 50Hz = new double[FilterDataNum_VAR], Sin100Hz = new double[FilterDataNum_VAR];
double[] InputData = new double[FilterDataNum_VAR];
int x = 0;
for (double i = 0;x< FilterDataNum_VAR; i=i+ FilterSamplingInterval,x++)
{
//Sin200Hz[x] = Math.Sin(Math.PI * 400 * i);
Sin300Hz[x] = Math.Sin(Math.PI * 600 * i);
/
/Sin1Hz[x] = Math.Sin(Math.PI * 2.0 * i);
Sin50Hz[x] = Math.Sin(Math.PI * 100 * i);
//Sin100Hz[x] = Math.Sin(Math.PI * 200 * i);
InputData[x] = Sin200Hz[x] + Sin300Hz[x] + Sin1Hz[x]+ Sin50Hz[x]+ Sin100Hz[x];
}
//var OutputData=FirstOrderLowPass(InputData);
//var OutputData = IIRLowPassFilter(InputData);
//var OutputData = IIRNotchFilter(InputData);
//var OutputData = IIRLowPassFilter_5(InputData);
var OutputData = IIRLowPassFilter_6(InputData);
//var OutputData = IIRBandPassFilter_6(InputData);
/
/OutputData = IIRHighPassFilter_3(OutputData);
//OutputData = IIRNotchFilter(OutputData);
//OutputData = IIRLowPassFilter(OutputData);
//var OutputData = FourthOrderRCLowpassFilter(InputData);
Tuple<double[], double[]> value = new Tuple<double[], double[]>(InputData, OutputData);
return value;
}
public double[] FirstOrderLowPass(double[] ProSingle) //⼀阶低通滤波
{
double[] FilteredData =new double[FilterDataNum_VAR]; //滤波后数据
double LowPassCoefficient = (1.0 / FilterDataFrequence_VAR) * 2 * 3.1415926 * LowPassCutoffFre_VAR; //低通滤波器系数
FilteredData[0] = ProSingle[0];
for(int i=1;i< FilterDataNum_VAR;i++)
{
FilteredData[i] =LowPassCoefficient * ProSingle[i] + ((1.0 - LowPassCoefficient) * FilteredData[i - 1]);
}
return FilteredData;
}
public double[] IIRLowPassFilter(double[] ProSingle) //IIR3阶低通滤波
{
double[] FilteredData = new double[FilterDataNum_VAR]; //滤波后数据
FilteredData[0] = ProSingle[0];
FilteredData[1] = ProSingle[1];
FilteredData[2] = ProSingle[2];
FilteredData[3] = ProSingle[3];
for (int i = 3; i < FilterDataNum_VAR; i++)
{
FilteredData[i] = 0.00009544*ProSingle[i] + 0.0002863 * ProSingle[i - 1] + 0.0002863*ProSingle[i - 2] + 0.00009844 * ProSingle[i - 3] - (-2.81157) * FilteredData[i - 1] - 2.64048 * FilteredData[i - 2] - (-0.828146) * FilteredData[i - 3];
}
return FilteredData;
}
public double[] IIRNotchFilter(double[] ProSingle) //IIR陷波
{
double[] FilteredData = new double[FilterDataNum_VAR]; //滤波后数据
double[] W = new double[3];
double[] IIR_B = new[] { 0.9984039, -1.995822607, 0.998403955 };
double[] IIR_A =new []{1, -1.995822607, 0.9968079114};
W[0] = 0; W[1] = 0; W[2] = 0;
for (int i = 0;i< FilterDataNum_VAR;i++)
{
W[0] = ProSingle[i] - IIR_A [1]* W[1] - IIR_A [2]* W[2];
FilteredData[i] = IIR_B[0] * W[0] + IIR_B[1] * W[1] + IIR_B[2] * W[2];
W[2] = W[1];
W[1] = W[0];
}
return FilteredData;
}
public double[] IIRLowPassFilter_5(double[] ProSingle) //IIR5阶低通滤波
{
double[] FilteredData = new double[FilterDataNum_VAR]; //滤波后数据
FilteredData[0] = ProSingle[0];
FilteredData[1] = ProSingle[1];
FilteredData[2] = ProSingle[2];
FilteredData[3] = ProSingle[3];
FilteredData[4] = ProSingle[4];
FilteredData[5] = ProSingle[5];
for (int i = 5; i < FilterDataNum_VAR; i++)
{
FilteredData[i] = 0.000000200242067359682832116980072401091 * ProSingle[i]
+ 0.000001001210336798414160584900362005456 * ProSingle[i - 1]
+ 0.000002002420673596828321169800724010912 * ProSingle[i - 2]
+ 0.000002002420673596828321169800724010912 * ProSingle[i - 3]
+ 0.000001001210336798414160584900362005456 * ProSingle[i - 4]
+ 0.000000200242067359682832116980072401091 * ProSingle[i - 5]
- (-4.695040626100188241309751902008429169655) * FilteredData[i - 1]
-
8.82614592256382479718013200908899307251 * FilteredData[i - 2]
- (-8.30396669308456125691009219735860824585) * FilteredData[i - 3]
- 3.909893994115305027037265972467139363289 * FilteredData[i - 4]
- (-0.737026189748224092035400190070504322648) * FilteredData[i - 5];
}
return FilteredData;
}
public double[] IIRLowPassFilter_6(double[] ProSingle) //IIR6阶低通滤波
{
double[] FilteredData = new double[FilterDataNum_VAR]; //滤波后数据
FilteredData[0] = ProSingle[0];
FilteredData[1] = ProSingle[1];
FilteredData[2] = ProSingle[2];
FilteredData[3] = ProSingle[3];
