巴特沃斯（Butterworth）滤波器的设计和幅频特性曲线绘制_matlab根据受扰音乐信号的频谱特性设计巴特沃斯低通滤波器,画出幅频特性图并得出

大佬的一些文章
设计滤波器的参考链接
绘制伯德图的源码参考链接

一、设计Butterworth滤波器

设计滤波器，其实就是设计传递函数，butterworth低通滤波器的传递函数如下：

对于butterworth高通滤波器，唯一的区别是分子项，从1变为s^N。
高低通Butterworth滤波器归一化后的系数如下：
（所谓归一化，就是只按照阶数来定系数，不考虑截止频率）

反归一化过程（考虑截止频率）：将s用以下公式代替，wc为截止频率，单位rad/s。

就得到考虑了阶数以及截至频率的Butterworth滤波器的传递函数。

1.1 利用matlab获取滤波器的归一化系数

在Matlab中，通过改变下面代码的N值（阶数），可以得到任意阶的归一化系数，（上面只给出了10阶的，你可以对比看看对不对）。

clc;clear
N = 10;%阶数
wc = 1;%截至频率设置为1rad/s，这样相当于归一化

[B,A]= butter(N,wc,'s')  %代入N和wn设计低通巴特沃斯模拟滤波器
%[B,A]= butter(N,wc,'high','s')  %代入N和wn设计高通巴特沃斯模拟滤波器
B=B./B(N);
A=A./B(N);
1
2
3
4
5
6
7
8

1.2 设计低通Butterworth滤波器

还记得上面说的吗？设计滤波器就是设计传递函数的分子项和分母项的系数。

 %% butterworth滤波器的归一化系数，共7阶
    coff=[1 1 0 0 0 0 0 0;
          1 sqrt(2) 1 0 0 0 0 0;
          1 2 2 1 0 0 0 0;
          1 2.6131 3.4142 2.6131 1 0 0 0;
          1 3.2361 5.2361 5.2361 3.2361 1 0 0;
          1 3.8637 7.4641 9.1416 7.4641 3.8637 1 0;
          1 4.4940 10.0978 14.5918 14.5918 10.0978 4.4940 1];
    N=4;%滤波器的阶数
    wc=90;%单位Hz
    wc=2*pi*wc;%转换为rad/s
    
    num1=[0 0 0 0 1];	%分子系数
    for i=1:(N+1)
        den1(i)=coff(N,i)*(1/wc)^(N-i+1);%分母系数
    end
    
    %整理系数
    num1=num1/den1(1);
    den1=den1/den1(1);
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

1.3 设计高通Butterworth滤波器

高低通的区别，仅仅只是分子项的区别。
对于低通，分子项为1；——以5阶为例为[0 0 0 0 1]
对于高通，分子项为s^n;——以5阶为例为[1 0 0 0 0]

 %% butterworth滤波器的归一化系数，共7阶
    coff=[1 1 0 0 0 0 0 0;
          1 sqrt(2) 1 0 0 0 0 0;
          1 2 2 1 0 0 0 0;
          1 2.6131 3.4142 2.6131 1 0 0 0;
          1 3.2361 5.2361 5.2361 3.2361 1 0 0;
          1 3.8637 7.4641 9.1416 7.4641 3.8637 1 0;
          1 4.4940 10.0978 14.5918 14.5918 10.0978 4.4940 1];
    N=4;%滤波器的阶数
    wc=90;%单位Hz
    wc=2*pi*wc;%转换为rad/s
    
    num1=[1 0 0 0 0];	%分子系数
    for i=1:(N+1)
        den1(i)=coff(N,i)*(1/wc)^(N-i+1);%分母系数
    end
    
    %整理系数
    num1=num1/den1(1);
    den1=den1/den1(1);
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

二、C++源码——绘制伯德图

头文件和cpp文件见最末尾，具体使用方法如下：

//在应用的头文件类中
#include "bode.h"

Bode HighPassCurveData;//高通滤波器数据

//在应用的cpp文件中
//butterworth滤波器的归一化系数，共8阶
double coff[8][9]={1,1,0,0,0,0,0,0,0,
               1,1.4142,1,0,0,0,0,0,0,
               1,2,2,1,0,0,0,0,0,
               1,2.6131,3.4142,2.6131,1,0,0,0,0,
               1,3.2361,5.2361,5.2361,3.2361,1,0,0,0,
               1,3.8637,7.4641,9.1416,7.4641,3.8637,1,0,0,
               1,4.4940,10.0978,14.5918,14.5918,10.0978,4.4940,1,0,
               1,5.1258309,13.1370712,21.846151,25.6883559,21.846151,13.1370712,5.1258309,1};

//更新高通滤波器的幅频曲线
HighPassCurveData._TF.n=OrderHighPass+1;
HighPassCurveData._TF.d=OrderHighPass+1;

//定义分子项系数
double tempvalue=1.0/(2*PI*CutOffFreqHighPass);
HighPassCurveData._TF.num[0]=pow(tempvalue,OrderHighPass);;
for(int i=1;i<HighPassCurveData._TF.n;++i)
{
     HighPassCurveData._TF.num[i]=0;
}

//定义分母项系数
for(int i=0;i<HighPassCurveData._TF.d;++i)
{
     HighPassCurveData._TF.den[i]=coff[OrderHighPass-1][i]*pow(tempvalue,OrderHighPass-i);
}

HighPassCurveData._wlen=LEN;//幅频特性曲线的数据长度
HighPassCurveData.freData  = new struct fre[LEN];
HighPassCurveData.BodeData =new struct BodeNum[LEN];
for(int i=0;i<LEN;i++)
{
     HighPassCurveData.freData[i].f=Xaxisdata[i];//幅频特性曲线的X轴坐标，单位Hz
     HighPassCurveData.freData[i].w=HighPassCurveData.freData[i].f*2*PI;//转换单位为rad/s
}
//计算幅频特性曲线各频点的幅度值
HighPassCurveData.compute();

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45

头文件如下：

#ifndef BODE_H
#define BODE_H
/******************************************************************************
* 文件       ：    bode.h
* 作者       ：     dhs 746769845@qq.com
* 版本       ：    V1.0
* 日期       ：    2020-7-31
* 描述       ：    给定传递函数绘制，计算对应的波特图坐标的数据（幅值、相角）
*
******************************************************************************/
#include <complex>
#include <cmath>
using namespace std;
 
#define PI 3.1415926535
/*传递函数结构体*/
struct TransferFunction
{
   double num[10]; //传递函数分子项
   double den[10]; //传递函数分母项
   char  n;       //分子个数
   char  d;       //分母个数
};
/*幅值、相角结构体*/
struct BodeNum
{
   double mag;      //幅值 db
   double phase;    //相角 度（角度）
};
/*频率、角频率结构体*/
struct fre
{
   double f;       //频率 Hz
   double w;       //角频率 rad/s
};
class Bode
{
public:
   Bode(struct TransferFunction TF);  //传入传递函数
   ~Bode();
   struct BodeNum *compute();         //完成计算
   struct fre *logspace(int start, int stop, int num=50, int base=10.0);//产生输入的频率数组
   int getWlen(){return _wlen;}      //获取数组长度
 
private:
   struct TransferFunction _TF;      //传递函数
   double *_w;                        //频率指针
   int _wlen;                        //频率长度
   struct fre *freData;              //频率、角频率
   struct BodeNum *BodeData;         //幅值、相角
};
 
#endif // BODE_H
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53

cpp文件如下：

 
/******************************************************************************
* 文件       ：    bode.h
* 作者       ：     dhs 746769845@qq.com
* 版本       ：    V1.0
* 日期       ：    2020-7-31
* 描述       ：    给定传递函数绘制，计算对应的波特图坐标的数据（幅值、相角）
*
******************************************************************************/
#include "bode.h"
 
 
/*******************************************************************************
* 函 数 名         : Bode
* 函数功能         : bode类，构造函数
* 输入参数         : TF          --> 传入传递函数
* 返 回 值         : 无
*******************************************************************************/
Bode::Bode(struct TransferFunction TF)
   :_TF(TF)
{
   freData  = 0;
   BodeData = 0;
}
/*******************************************************************************
* 函 数 名         : compute()
* 函数功能         : bode图数据计算
* 输入参数         : 无
* 返 回 值         : bode图计算后的数据（数组指针）
*******************************************************************************/
struct BodeNum *Bode::compute()
{
   int i=0;
   int j=0;
   complex<double> ds,ms,s,j1;
   if(freData ==0)
   {
       return 0;
   }
   ds=ms={0,0};
   j1={0,1};                //虚数单位
   for (i=0; i<_wlen; i++)
   {
       s=j1*freData[i].w; //传递函数中的 s用jw代入
       ms ={0,0};         //保存分子项
       ds ={0,0};         //保存分母项
       for (j=0; j<_TF.n; j++)
       {
           ms= ms * s + _TF.num[j];
       }
       for (j=0; j<_TF.d; j++)
       {
           ds= ds * s + _TF.den[j];
       }
       s = ms/ds;
       BodeData[i].mag = 20.0 * log10(abs(s));       //20倍 log10幅度值
       BodeData[i].phase = atan2(s.imag(),s.real()) * 180.0 / PI;  //相角值
   }
   return BodeData;
}
/*******************************************************************************
* 函 数 名         : logspace
* 函数功能         : 生成频率
* 输入参数         : start-->对数开始的冥
*                   stop -->对数结束的冥
*                   num  --> 生成频率的个数 默认50
*                   base --> 对数的底  默认为10
*                  例：logspace(-2, 5, 200, 10)  //表示从0.01Hz~0.1MHz按照10倍频的间距 生成200个频率
* 返 回 值         : 返回生成的频率（数组指针）
*******************************************************************************/
struct fre *Bode::logspace(int start, int stop, int num, int base)
{
   double step=0;
   int i =0;
   _wlen = num;
   freData  = new struct fre[num];
   BodeData =new struct BodeNum[num];
   step =(stop-start)/(num-1.0);
   for(i=0;i<num; i++)
   {
       freData[i].f=pow(base, start+step*i);
       freData[i].w=freData[i].f*2*PI;
   }
 
   return freData;
 
}
/*******************************************************************************
* 函 数 名         : ~Bode
* 函数功能         : 析构
* 输入参数         :
* 返 回 值         :
*******************************************************************************/
Bode::~Bode()
{
   delete BodeData;
   delete freData;
 
}
 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100