巴特沃斯滤波器c语言.docx
《巴特沃斯滤波器c语言.docx》由会员分享，可在线阅读，更多相关《巴特沃斯滤波器c语言.docx（23页珍藏版）》请在冰豆网上搜索。
巴特沃斯滤波器c语言
巴特沃斯滤波器c语言
1.模拟滤波器的设计
1.1巴特沃斯滤波器的次数
根据给定的参数设计模拟滤波器，然后进行变数变换，求取数字滤波器的方法，称为滤波器的间接设计。
做为数字滤波器的设计基础的模拟滤波器，称之为原型滤波器。
这里，我们首先介绍的是最简单最基础的原型滤波器，巴特沃斯低通滤波器。
由于IIR滤波器不具有线性相位特性，因此不必考虑相位特性，直接考虑其振幅特性。
在这里，N是滤波器的次数，Ωc是截止频率。
从上式的振幅特性可以看出，这个是单调递减的函数，其振幅特性是不存在纹波的。
设计的时候，一般需要先计算跟所需要设计参数相符合的次数N。
首先，就需要先由阻带频率，计算出阻带衰减
将巴特沃斯低通滤波器的振幅特性，直接带入上式，则有
最后，可以解得次数N为
当然，这里的N只能为正数，因此，若结果为小数，则舍弃小数，向上取整。
1.2巴特沃斯滤波器的传递函数
巴特沃斯低通滤波器的传递函数，可由其振幅特性的分母多项式求得。
其分母多项式
根据S解开，可以得到极点。
这里，为了方便处理，我们分为两种情况去解这个方程。
当N为偶数的时候，
这里，使用了欧拉公式。
同样的，当N为奇数的时候，
同样的，这里也使用了欧拉公式。
归纳以上，极点的解为
上式所求得的极点，是在s平面内，在半径为Ωc的圆上等间距的点，其数量为2N个。
为了使得其IIR滤波器稳定，那么，只能选取极点在S平面左半平面的点。
选定了稳定的极点之后，其模拟滤波器的传递函数就可由下式求得。
1.3巴特沃斯滤波器的实现（C语言）
首先，是次数的计算。
次数的计算，我们可以由下式求得。
其对应的C语言程序为
[cpp] viewplaincopy
1.N = Ceil（0.5*（ log10 （ pow （10, Stopband_attenuation/10） - 1） /
2. log10 （Stopband/Cotoff） ））;
然后是极点的选择，这里由于涉及到复数的操作，我们就声明一个复数结构体就可以了。
最重要的是，极点的计算含有自然指数函数，这点对于计算机来讲，不是太方便，所以，我们将其替换为三角函数，
这样的话，实部与虚部就还可以分开来计算。
其代码实现为
[cpp] viewplaincopy
1.typedef struct
2.{
3. double Real_part;
4. double Imag_Part;
5.} COMPLEX;
6.
7.
8.COMPLEX poles[N];
9.
10.for（k = 0;k <= （（2*N）-1） ; k++）
11.{
12. if（Cotoff*cos（（k+dk）*（pi/N）） < 0）
13. {
14. poles[count].Real_part = -Cotoff*cos（（k+dk）*（pi/N））;
15. poles[count].Imag_Part= -Cotoff*sin（（k+dk）*（pi/N））;
16. count++;
17. if （count == N） break;
18. }
19.}
计算出稳定的极点之后，就可以进行传递函数的计算了。
传递的函数的计算，就像下式一样
这里，为了得到模拟滤波器的系数，需要将分母乘开。
很显然，这里的极点不一定是整数，或者来说，这里的乘开需要做复数运算。
其复数的乘法代码如下，
[cpp] viewplaincopy
1.int Complex_Multiple（COMPLEX a,COMPLEX b,
2. double *Res_Real,double *Res_Imag）
3.
4.{
5. *（Res_Real） = （a.Real_part）*（b.Real_part） - （a.Imag_Part）*（b.Imag_Part）;
6. *（Res_Imag）= （a.Imag_Part）*（b.Real_part） + （a.Real_part）*（b.Imag_Part）;
7. return （int）1;
8.}
有了乘法代码之后，我们现在简单的情况下，看看其如何计算其滤波器系数。
我们做如下假设
这个时候，其传递函数为
将其乘开，其大致的关系就像下图所示一样。
计算的关系一目了然，这样的话，实现就简单多了。
高阶的情况下也一样，重复这种计算就可以了。
其代码为
[cpp] viewplaincopy
1. Res[0].Real_part = poles[0].Real_part;
2. Res[0].Imag_Part= poles[0].Imag_Part;
3. Res[1].Real_part = 1;
4. Res[1].Imag_Part= 0;
5.
6.for（count_1 = 0;count_1 < N-1;count_1++）
7.{
8. for（count = 0;count <= count_1 + 2;count++）
9. {
10. if（0 == count）
11. {
12. Complex_Multiple（Res[count], poles[count_1+1],
13. &（Res_Save[count].Real_part）,
14. &（Res_Save[count].Imag_Part））;
15. }
16. else if（（count_1 + 2） == count）
17. {
18. Res_Save[count].Real_part += Res[count - 1].Real_part;
19. Res_Save[count].Imag_Part += Res[count - 1].Imag_Part;
20. }
21. else
22. {
23. Complex_Multiple（Res[count], poles[count_1+1],
24. &（Res_Save[count].Real_part）,
25. &（Res_Save[count].Imag_Part））;
26.1 Res_Save[count].Real_part += Res[count - 1].Real_part;
27. Res_Save[count].Imag_Part += Res[count - 1].Imag_Part;
28. }
29. }
30. *（b+N） = *（a+N）;
到此，我们就可以得到一个模拟滤波器巴特沃斯低通滤波器了。
2.双1次z变换
2.1双1次z变换的原理
我们为了将模拟滤波器转换为数字滤波器的，可以用的方法很多。
这里着重说说双1次z变换。
我们希望通过双1次z变换，建立一个s平面到z平面的映射关系，将模拟滤波器转换为数字滤波器。
和之前的例子一样，我们假设有如下模拟滤波器的传递函数。
将其做拉普拉斯逆变换，可得到其时间域内的连续微分方程式，
其中，x（t）表示输入，y（t）表示输出。
然后我们需要将其离散化，假设其采样周期是T，用差分方程去近似的替代微分方程，可以得到下面结果
然后使用z变换，再将其化简。
可得到如下结果
从而，我们可以得到了s平面到z平面的映射关系，即
由于所有的高阶系统都可以视为一阶系统的并联，所以，这个映射关系在高阶系统中，也是成立的。
然后，将关系式
带入上式，可得
到这里，我们可以就可以得到Ω与ω的对应关系了。
这里的Ω与ω的对应关系很重要。
我们最终的目的设计的是数字滤波器，所以，设计时候给的参数必定是数字滤波器的指标。
而我们通过间接设计设计IIR滤波器时候，首先是要设计模拟滤波器，再通过变换，得到数字滤波器。
那么，我们首先需要做的，就是将数字滤波器的指标，转换为模拟滤波器的指标，基于这个指标去设计模拟滤波器。
另外，这里的采样时间T的取值很随意，为了方便计算，一般取1s就可以。
2.2双1次z变换的实现（C语言）
我们设计好的巴特沃斯低通滤波器的传递函数如下所示。
我们将其进行双1次z变换，我们可以得到如下式子
可以看出，我们还是需要将式子乘开，进行合并同类项，这个跟之前说的算法相差不大。
其代码为。
[cpp] viewplaincopy
1.for（Count = 0;Count<=N;Count++）
2. {
3. for（Count_Z = 0;Count_Z <= N;Count_Z++）
4. {
5. Res[Count_Z] = 0;
6. Res_Save[Count_Z] = 0;
7. }
8. Res_Save [0] = 1;
9. for（Count_1 = 0; Count_1 < N-Count;Count_1++）
10. {
11. for（Count_2 = 0; Count_2 <= Count_1+1;Count_2++）
12. {
13. if（Count_2 == 0） Res[Count_2] += Res_Save[Count_2];
14. else if（（Count_2 == （Count_1+1））&&（Count_1 !
= 0））
15. Res[Count_2] += -Res_Save[Count_2 - 1];
16. else Res[Count_2] += Res_Save[Count_2] - Res_Save[Count_2 - 1];
17. for（Count_Z = 0;Count_Z<= N;Count_Z++）
18. {
19. Res_Save[Count_Z] = Res[C
升级会员 