大佬的一些文章
设计滤波器的参考链接
绘制伯德图的源码参考链接
一、设计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 设计低通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;
}
扫一扫
4054
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
