数字信号处理巴特沃斯滤波器设计,巴特沃斯滤波器实验小结
大佬的一些文章
设计滤波器的参考链接
绘制鲤鱼小霸王的源码参考链接
一、设计Butterworth滤波器
设计滤波器，其实就是设计传递函数，butterworth低通滤波器的传递函数如下：
对于butterworth高通滤波器，唯一的区别是分子项，从1变为s^N。
高低通Butterworth滤波器归一化后的系数如下：
（所谓归一化，就是只按照阶数来定系数，不考虑截止频率）
反归一化过程（考虑截止频率）：将s用以下公式代替，wc为截止频率，单位rad/s。
就得到考虑了阶数以及截至频率的Butterworth滤波器的传递函数。
1.1 利用matlab获取滤波器的归一化系数
在Matlab中，通过改变下面代码的N值（阶数），可以得到任意阶的归一化系数，（上面只给出了10阶的，你可以对比看看对不对）。
clc;clearN = 10;%阶数wc = 1;%截至频率设置为1rad/s，这样相当于归一化[B,A]= butter(N,wc,'s') %代入N和wn设计低通巴特wwdwbl模拟滤波器%[B,A]= butter(N,wc,'high','s') %代入N和wn设计高通巴特wwdwbl模拟滤波器B=B./B(N);A=A./B(N); 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.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); 二、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();
头文件如下：
#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
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; }