一、背景
这种滤波器最先由英国工程师斯蒂芬·巴特沃斯(Stephen Butterworth)在1930年发表在英国《无线电工程》期刊的一篇论文中提出的。
二、特点
巴特沃斯滤波器的特点是通频带内的频率响应曲线最大限度平坦，没有起伏，而在阻频带则逐渐下降为零。 在振幅的对数对角频率的波特图上，从某一边界角频率开始，振幅随着角频率的增加而逐步减少,趋向负无穷大。一阶巴特沃斯滤波器的衰减率为每倍频6分贝，每十倍频20分贝。二阶巴特沃斯滤波器的衰减率为每倍频12分贝、三阶巴特沃斯滤波器的衰减率为每倍频18分贝、如此类推。巴特沃斯滤波器的振幅对角频率单调下降，并且也是唯一的无论阶数，振幅对角频率曲线都保持同样的形状的滤波器。只不过滤波器阶数越高，在阻频带振幅衰减速度越快。其他滤波器高阶的振幅对角频率图和低阶数的振幅对角频率有不同的形状。
由图可见，当N趋于无穷时，将会得到一个理想的低通滤波器。
三、原理
1.幅值平方
2.拉普拉斯变换
令
其中，
所以有
重新带入幅值平方
其中
设数组sb的各个元素为上式分母中关于s的幂次方项的系数
3.Z域变换（离散化）
由
得到
3.1线性化
试图使用泰勒公式对s进行展开，但是lnz在z=0处没有定义，于是令
于是
对上式在z=0处泰勒展开，有
带入s的定义，整理可得
展开没看明白的可以看一下张宇的基础三十讲（我真不是托）
取一阶近似，得
3.2双线性变换
双线性变换就是使用这种一阶估计法，将连续时间传递函数H(s)中的s替换成离散域中的z变量
即
进行整理后，最终得到
最终为
上式中的序列num、den就是离散化的目标。
四、巴特沃斯滤波器的设计与实现
1.确定参数
需要确定的参数如下
double *passF, //通带截止频率（滤波允许通过的范围）
double *stopF, //阻带截止频率（衰减的信号，不会出现在输出端输出）
double rp, //通带最大衰减（dB）(通带范围内的最小值与通带截止频率的差值)
double rs, //阻带最大衰减（dB）
// double Fs,
double fs, //采样频率
int filterType); //返回值：返回巴特沃斯低通滤波器的阶数N和截止频率Ws结构体
2.计算预畸参数
双线性变换具有从根本上避免脉冲响应不变法中的频率混叠现象，缺点是引入了频率失真，因此，引入预畸变的概念，即在双线性变换之前，进行预畸来校正频率失真的问题，使得进行双线性变换之后的频率与预期一致。
2.1预畸原理
由
进行变量代换
得到
所以有
2.2计算第一步确定的参数预畸值
3.计算低通滤波器阶数
由
可得，阶数N
4.计算截止频率\Omega_c
由
得
5.确定分之分母系数数组长度length
6.查表，获取拉式变换下传递函数的分母系数
即计算
其中，关于a可查表获得
(Note: 参考 陈佩青《数字信号处理教程》第二版266页 表6-4)
最终返回数组sb数组sb的各个元素为上式分母中关于s的幂次方项的系数
代码如下
//选择预畸方法
switch (butterValue.filterType)
{
case FILTER_IIR_BUTTER_LOW:
Wc = fabs(tan(butterValue.fc * M_PI / butterValue.fs));
break;
case FILTER_IIR_BUTTER_HIGH:
Wc = fabs(1/tan(butterValue.fc * M_PI / butterValue.fs));
break;
case FILTER_IIR_BUTTER_PASS:
Wc = 2 * M_PI * butterValue.fc / butterValue.fs;
Wc = fabs((butterValue.cosW0 - cos(Wc))/sin(Wc));
break;
case FILTER_IIR_BUTTER_STOP:
Wc = 2 * M_PI * butterValue.fc / butterValue.fs;
Wc = fabs(sin(Wc)/(cos(Wc) - butterValue.cosW0));
break;
default:
break;
}
for(int i = 0; i < length; i++)
{
returnSb[i] = g_butterPb[length - 1][i] * pow(Wc, length-i); //计算系数
}
returnSb[length] = 1.0; //最高次幂的系数为1
}7.双线性变换
根据前面的式子可得
7.1利用杨辉三角计算多项式
代码如下
void pascalTriangle(
int N,
int symbol,
int *vector)
{
vector[0] = 1;
if(N == 0)
{
return;
}
else if (N == 1)
{
if(symbol == SYMBOL_ADD)
{
vector[1] = 1;
}
else
{
vector[0] = -1; //如果是减号，则第二项系数是-1
vector[1] = 1;
}
return;
}
int length = N + 1; //数组长度
int *temp = (int *)malloc(sizeof(int) * length);
/*
int temp[length]; //定义中间变量
*/
temp[0] = 1;
temp[1] = 1;
for(int i = 2; i <= N; i++)
{
vector[i] = 1;
for(int j = 1; j < i; j++)
{
vector[j] = temp[j - 1] + temp[j]; //x[m] = x[m-1][n-1] + x[m-1]
}
if(i == N) //最后一次不需要给中间变量赋值
{
if(symbol == SYMBOL_SUB) //运算符为减号
{
for(int k = 0; k < length; k++)
{
vector[k] = vector[k] * pow((float)-1, length - 1 - k);
}
}
return;
}
for(int j = 1; j <= i; j++)
{
temp[j] = vector[j];
}
}
if (temp != NULL)
free(temp);
}7.3计算H（z）分母系数
分母形式：
其中
................
整体源码如下：
#include <stdlib.h>
#include <malloc.h>
#include "BTWF.h"
#define M_PI 3.14159265
#define MAX(a, b) ((a) > (b) ? (a) : (b))
#define MIN(a, b) ((a) < (b) ? (a) : (b))
//定义巴特沃斯滤波器pb系数列表（b0,b1,...,bn）
static const double g_butterPb[10][10] =
{
/* N = 1 */ {1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
/* N = 2 */ {1.0, 1.4142136, 0, 0, 0, 0, 0, 0, 0, 0},
/* N = 3 */ {1.0, 2.0000000, 2.0000000, 0, 0, 0, 0, 0, 0, 0},
/* N = 4 */ {1.0, 2.6131259, 3.4142136, 2.6131259, 0, 0, 0, 0, 0, 0},
/* N = 5 */ {1.0, 3.236068, 5.236068, 5.236068, 3.236068, 0, 0, 0, 0, 0},
/* N = 6 */ {1.0, 3.8637033, 7.4641016, 9.1416202, 7.4641016, 3.8637033, 0, 0, 0, 0},
/* N = 7 */ {1.0, 4.4939592, 10.0978347, 14.5917939, 14.5917939, 10.0978347, 4.4939592, 0, 0, 0},
/* N = 8 */ {1.0, 5.1258309, 13.1370712, 21.8461510, 25.6883559, 21.8461510, 13.1370712, 5.1258309, 0, 0},
/* N = 9 */ {1.0, 5.7587705, 16.5817187, 31.1634375, 41.9863857, 41.9863857, 31.1634375, 16.5817187, 5.7587705, 0},
/* N = 10 */{1.0, 6.3924532, 20.4317291, 42.8020611, 64.8823963, 74.2334292, 64.8823963, 42.8020611, 20.4317291, 6.3924532}
};
void pascalTriangle(
int N,
int symbol,
int *vector);
void coefficientEquation(
int *originalCoef,
int N,
int *nextCoef,
int nextN);
void coefficientEquation2(
double *originalCoef,
int N,
double *nextCoef,
int nextN);
/*======================================================================
* 方法名： filterIIRButterLowpass
* 方法功能：设计巴特沃斯样本低通示波器
*
* 变量名称：
* fpass - 通带截止频率（模拟频率）
* fstop - 阻带截止频率（模拟频率）
* rp - 通带最大衰减（dB）
* rs - 阻带最小衰减（dB）
* Fs - 采样频率
*
* 返回值： 返回巴特沃斯低通滤波器的阶数N和截止频率Ws结构体
*=====================================================================*/
ButterFilterStruct filterIIRButterLowpass(
double *passF,
double *stopF,
double rp,
double rs,
// double Fs,
double fs,
int filterType)
{
ButterFilterStruct nAndFc; //返回滤波器的阶数N和截止频率fc
nAndFc.filterType = filterType; //滤波器类型
double nOfN = 0.0;
double passW = 0.0, stopW = 0.0, wa, wc; //wa = stopW/passW或其导数
double passF1 = 0.0, passF2 = 0.0, stopF1 = 0.0, stopF2 = 0.0, w0 = 0.0;//w0 - 中心频率
double passW1 = 0.0, passW2 = 0.0, stopW1 = 0.0, stopW2 = 0.0, fc = 0.0;
rs = fabs(rs);
rp = fabs(rp);
passF1 = passF[0];
stopF1 = stopF[0];
//根据滤波器类型，选择不同的预畸变换式
switch (filterType)
{
case FILTER_IIR_BUTTER_LOW:
if(passF1 >= stopF1)
{
nAndFc.isFOK = FALSE;
printf("错误！应满足：passF < stopF \n");
return nAndFc;
}
nAndFc.isFOK = TRUE;
passW = tan(passF1 * M_PI / fs); //数字低通，频率预畸，W = tan(w/2)
stopW = tan(stopF1 * M_PI / fs);
wa = fabs(stopW/passW);
break;
case FILTER_IIR_BUTTER_HIGH:
if(passF1 <= stopF1)
{
nAndFc.isFOK = FALSE;
printf("错误！应满足：passF > stopF \n");
return nAndFc;
}
nAndFc.isFOK = TRUE;
passW = 1/tan(passF1 * M_PI / fs); //数字高通，频率预畸，W = cot(w/2)
stopW = 1/tan(stopF1 * M_PI / fs);
wa = fabs(stopW/passW);
break;
case FILTER_IIR_BUTTER_PASS:
passF2 = passF[1];
stopF2 = stopF[1];
if(!(stopF1 < passF1 && passF1 < passF2 && passF2 < stopF2))
{
nAndFc.isFOK = FALSE;
printf("错误！应满足：stopF[1] < passF[1] < passF[2] < stopF[2]");
return nAndFc;
}
nAndFc.isFOK = TRUE;
//转换为数字频率（不进行预畸）
passW1 = 2 * M_PI * passF1 / fs;
passW2 = 2 * M_PI * passF2 / fs;
stopW1 = 2 * M_PI * stopF1 / fs;
stopW2 = 2 * M_PI * stopF2 / fs;
nAndFc.cosW0 = cos((passW1 + passW2)/2)/cos((passW1 - passW2)/2); //保存cos(w0)
w0 = acos(nAndFc.cosW0);//求带通滤波器的中心频率
passW1 = (cos(w0)-cos(passW1))/sin(passW1); //通带截止频率
passW2 = (cos(w0)-cos(passW2))/sin(passW2);
stopW1 = (cos(w0)-cos(stopW1))/sin(stopW1);
stopW2 = (cos(w0)-cos(stopW2))/sin(stopW2);
passW = MAX(passW1, passW2); //通带截止频率
stopW = MIN(stopW1, stopW2); //阻带截止频率
wa = fabs(stopW/passW);
break;
case FILTER_IIR_BUTTER_STOP:
passF2 = passF[1];
stopF2 = stopF[1];
if(!(passF1 < stopF1 && stopF1 < stopF2 && stopF2 < passF2))
{
nAndFc.isFOK = FALSE;
printf("错误！应满足：passF[1] < stopF[1] < stopF[2] < passF[2]");
return nAndFc;
}
nAndFc.isFOK = TRUE;
//转换为数字频率（不进行预畸）
passW1 = 2 * M_PI * passF1 / fs;
passW2 = 2 * M_PI * passF2 / fs;
stopW1 = 2 * M_PI * stopF1 / fs;
stopW2 = 2 * M_PI * stopF2 / fs;
nAndFc.cosW0 = cos((stopW1 + stopW2)/2)/cos((stopW1 - stopW2)/2); //保存cos(w0)
w0 = acos(nAndFc.cosW0);//求带通滤波器的中心频率
passW1 = sin(passW1)/(cos(passW1)-nAndFc.cosW0); //通带截止频率
passW2 = sin(passW2)/(cos(passW2)-nAndFc.cosW0);
stopW1 = sin(stopW1)/(cos(stopW1)-nAndFc.cosW0);
stopW2 = sin(stopW2)/(cos(stopW2)-nAndFc.cosW0);
passW = MAX(passW1, passW2); //通带截止频率
stopW = MIN(stopW1, stopW2); //阻带截止频率
wa = fabs(stopW/passW);
break;
default:
break;
}
nAndFc.fs = fs; //采样频率
nAndFc.N = ceil(0.5 * log10((pow(10, 0.1*rs)-1)/(pow(10, 0.1*rp)-1))/log10(wa)); //计算N
nOfN = (float)nAndFc.N; //将N转化为float型
//根据滤波器类型，选择不同的预畸变换式
switch (filterType)
{
case FILTER_IIR_BUTTER_LOW:
wc = stopW / pow((pow(10, 0.1*rs) - 1), 1/(2*nOfN));
nAndFc.fc = fs/M_PI*atan(wc); //计算截止频率(3dB)Hz
nAndFc.length = nAndFc.N + 1; //系数数组长度
break;
case FILTER_IIR_BUTTER_HIGH:
wc = stopW / pow((pow(10, 0.1*rs) - 1), 1/(2*nOfN));
//wc = passW / pow((pow(10, 0.1*rp) - 1), 1/(2*nOfN));
nAndFc.fc = fs/M_PI*atan(1/wc); //计算截止频率(3dB)Hz
nAndFc.length = nAndFc.N + 1; //系数数组长度
break;
case FILTER_IIR_BUTTER_PASS:
wc = stopW1 / pow((pow(10, 0.1*rs) - 1), 1/(2*nOfN));
fc =asin((2*cos(w0)*wc + sqrt(pow(2*cos(w0)*wc, 2)-4*(wc*wc+1)*(cos(w0)*cos(w0)-1)))/(2*wc*wc+2));
// wc = passW1 / pow((pow(10, 0.1*rp) - 1), 1/(2*nOfN));
// fc =asin((2*cos(w0)*wc + sqrt(pow(2*cos(w0)*wc, 2)-4*(wc*wc+1)*(cos(w0)*cos(w0)-1)))/(2*wc*wc+2));
nAndFc.fc = fs / (2*M_PI) * fc;
nAndFc.length = 2 * nAndFc.N + 1; //系数数组长度
break;
case FILTER_IIR_BUTTER_STOP:
wc = -1/(stopW1 / pow((pow(10, 0.1*rs) - 1), 1/(2*nOfN)));
fc =asin((2*cos(w0)*wc + sqrt(pow(2*cos(w0)*wc, 2)-4*(wc*wc+1)*(cos(w0)*cos(w0)-1)))/(2*wc*wc+2));
nAndFc.fc = fs / (2*M_PI) * fc;
nAndFc.length = 2 * nAndFc.N + 1; //系数数组长度
break;
default:
break;
}
return nAndFc;
}
/*======================================================================
* 方法名： butterSbValue
* 方法功能：计算巴特沃斯滤波器分母多项式H(s)的系数Sb，注意：分子为Wc^N
* 说明： Sb[k] = Wc^(N-k) * Pb，其中Pb是归一化的分母多项式的根，可查表得到
* 系数由低次向高次排列
*
* 变量名称：
* butterValue - 存放滤波器参数（阶数和截止频率）的结构体变量
* returnSb - 计算结果
*
* 返回值： void
*=====================================================================*/
void butterSbValue(
ButterFilterStruct butterValue,
double *returnSb)
{
int length = butterValue.N; //滤波器阶数
double Wc = 0.0; //滤波器的截止频率
//选择预畸方法
switch (butterValue.filterType)
{
case FILTER_IIR_BUTTER_LOW:
Wc = fabs(tan(butterValue.fc * M_PI / butterValue.fs));
break;
case FILTER_IIR_BUTTER_HIGH:
Wc = fabs(1/tan(butterValue.fc * M_PI / butterValue.fs));
break;
case FILTER_IIR_BUTTER_PASS:
Wc = 2 * M_PI * butterValue.fc / butterValue.fs;
Wc = fabs((butterValue.cosW0 - cos(Wc))/sin(Wc));
break;
case FILTER_IIR_BUTTER_STOP:
Wc = 2 * M_PI * butterValue.fc / butterValue.fs;
Wc = fabs(sin(Wc)/(cos(Wc) - butterValue.cosW0));
break;
default:
break;
}
for(int i = 0; i < length; i++)
{
returnSb[i] = g_butterPb[length - 1][i] * pow(Wc, length-i); //计算系数
}
returnSb[length] = 1.0; //最高次幂的系数为1
}
/*======================================================================
* 方法名： butterLowOrHigh
* 方法功能：计算巴特沃斯低通（高通）滤波器系统方法的系数，包括分子和分母系数
*
* 变量名称：
* butterValue - 存放滤波器参数（阶数和截止频率）的结构体变量
* sb - 传入的模拟滤波器的系数，即H(s)的分母系数
* numerator - 计算后的分子系数数组
* denominator - 计算后的分母系数数组
*
* 返回值： void
*=====================================================================*/
void butterLowOrHigh(
ButterFilterStruct butterValue,
double *sb,
double *numerator,
double *denominator)
{
int length = butterValue.N; //滤波器阶数
int *tempCoef1 = (int *)malloc(sizeof(int) * (length+1));
int *tempCoef2 = (int *)malloc(sizeof(int) * (length+1));
/*
int tempCoef1[length + 1]; //定义系数数组，用于存放1 - z^(-1)、1 + z^(-1)每项次幂（0-N）系数，最低次在第一个
int tempCoef2[length + 1];
*/
int otherN; //1+z^(-1)的次数
double Fsx2 = 1; //butterValue.fs * 2; //计算2/T
for(int i = 0; i<= length; i++)
{
numerator[i] = 0.0; //初始化numerator和denominator
denominator[i] = 0.0;
}
for(int i = 0; i <= length; i++)
{
for(int j = 0; j<= length; j++)
{
tempCoef1[j] = 0; //tempCoef1和tempCoef2进行初始化
tempCoef2[j] = 0;
}
otherN = length - i;
if(butterValue.filterType == FILTER_IIR_BUTTER_LOW)
{
pascalTriangle(i, SYMBOL_SUB, tempCoef1); //利用杨辉三角计算1 - z^(-1)幂的系数
pascalTriangle(otherN, SYMBOL_ADD, tempCoef2); //利用杨辉三角计算1 + z^(-1)幂的系数
}
else
{
pascalTriangle(i, SYMBOL_ADD, tempCoef1); //利用杨辉三角计算1 + z^(-1)幂的系数
pascalTriangle(otherN, SYMBOL_SUB, tempCoef2); //利用杨辉三角计算1 - z^(-1)幂的系数
}
coefficientEquation(tempCoef1, i+1, tempCoef2, otherN+1); //两个多项式相乘，求其系数
for(int j = 0; j <= length; j++)
{
denominator[j] += pow(Fsx2, i) * (float)tempCoef1[length - j] * sb[i];
}
//分子系数
if(i == 0)
{
for(int j = 0; j <= length; j++)
{
numerator[j] = sb[0] * tempCoef2[length - j];
}
}
}
//系数归一化，分母的常数项为1
for(int i = length; i >= 0; i--)
{
numerator[i] = numerator[i] / denominator[0];
denominator[i] = denominator[i] / denominator[0];
}
if (tempCoef1 != NULL)
free(tempCoef1);
if (tempCoef2 != NULL)
free(tempCoef2);
}
/*======================================================================
* 方法名： butterPassOrStop
* 方法功能：计算巴特沃斯带通（带阻）滤波器系统方法的系数，包括分子和分母系数
*
* 变量名称：
* butterValue - 存放滤波器参数（阶数和截止频率）的结构体变量
* sb - 传入的模拟滤波器的系数，即H(s)的分母系数
* numerator - 计算后的分子系数数组
* denominator - 计算后的分母系数数组
*
* 返回值： void
*=====================================================================*/
void butterPassOrStop(
ButterFilterStruct butterValue,
double *sb,
double *numerator,
double *denominator)
{
int length = butterValue.length; //滤波器系数长度
int *tempCoef1 = (int *)malloc(sizeof(int) * length);
double *tempCoef2 = (double *)malloc(sizeof(double) * length);
double *tempCoef3 = (double *)malloc(sizeof(double) * length);
/*
int tempCoef1[length]; //定义系数数组，用于存放1 - z^(-2)、1 - 2*cos(w0)*z^(-1) + z^(-2)每项次幂（0-N）系数，最低次在第一个
double tempCoef2[length];
double tempCoef3[length];
*/
double tempCoef[3];
int otherN; //1+z^(-1)的次数（pass）,1 - 2*cos(w0)*z^(-1) + z^(-2)的次数（stop）
double Fsx2 = 1;//butterValue.fs * 2; //计算2/T = 1
for(int i = 0; i < length; i++)
{
numerator[i] = 0.0; //初始化numerator和denominator
denominator[i] = 0.0;
tempCoef1[i] = 0; //tempCoef1和tempCoef2进行初始化
tempCoef2[i] = 0.0;
tempCoef3[i] = 0.0;
}
tempCoef[0] = 1.0;
tempCoef[1] = -2.0 * butterValue.cosW0;
tempCoef[2] = 1.0;
//----------计算分子系数-----------
if(butterValue.filterType == FILTER_IIR_BUTTER_PASS) //带通滤波器
{
pascalTriangle(butterValue.N, SYMBOL_SUB, tempCoef1); //利用杨辉三角计算1 - z^(-1)幂的系数
for(int i = 0; i < length; i++) //变为1 - z^(-2)幂的系数，填充奇次幂0
{
int temp = i%2; //判断i奇偶
if(!temp) //偶次幂不为0
numerator[i] = sb[0] * tempCoef1[butterValue.N - i/2];
else
numerator[i] = 0.0;
}
}
else //带阻滤波器
{
tempCoef3[0] = 1.0; //1 - 2*cos(w0)*z^(-1) + z^(-2)的系数1,-2cos(w0),1
tempCoef3[1] = -2.0 * butterValue.cosW0;
tempCoef3[2] = 1.0;
for(int j = 1; j < butterValue.N; j++)
{
coefficientEquation2(tempCoef3, j*2+1, tempCoef, 3);
}
for(int i = 0; i < length; i++)
{
numerator[i] = sb[0] * tempCoef3[length - i - 1];
}
}
//----------计算分母系数,计算每一加数的系数----------
for(int i = 0; i <= butterValue.N; i++)
{
if(butterValue.filterType == FILTER_IIR_BUTTER_PASS)
{
otherN = butterValue.N - i;
}
else
{
otherN = i;
}
for(int j = 0; j < length; j++)
{
tempCoef1[j] = 0; //tempCoef1、tempCoef2和tempCoef3进行初始化
tempCoef2[j] = 0.0;
tempCoef3[j] = 0.0;
}
tempCoef3[0] = 1.0;
if(butterValue.N - otherN > 0) //当第0次相乘时，第一项为1，其余为0
{
tempCoef3[1] = -2.0 * butterValue.cosW0;
tempCoef3[2] = 1.0;
}
pascalTriangle(otherN, SYMBOL_SUB, tempCoef1); //利用杨辉三角计算1 - z^(-1)幂的系数
for(int j = 0; j < otherN*2+1; j++) //变为1 - z^(-2)幂的系数，填充奇次幂0
{
int temp = j%2; //判断i奇偶
if(!temp) //偶次幂不为0
{
tempCoef2[j] = (float)tempCoef1[j/2];
tempCoef1[j/2] = 0;
}
else
tempCoef2[j] = 0.0;
}
//利用多项式相乘法，计算1 - 2*cos(w0)*z^(-1) + z^(-2)幂的系数,j表示第几次相乘
for(int j = 1; j < butterValue.N - otherN; j++)
{
coefficientEquation2(tempCoef3, j*2+1, tempCoef, 3);
}
coefficientEquation2(tempCoef3, (butterValue.N - otherN)*2+1, tempCoef2, 2*otherN+1); //两个多项式相乘，求其系数
for(int j = 0; j < length; j++)
{
denominator[j] += pow(Fsx2, i) * tempCoef3[length - j - 1] * sb[i];
}
}
//系数归一化，分母的常数项为1
for(int i = length - 1; i >= 0; i--)
{
numerator[i] = numerator[i] / denominator[0];
denominator[i] = denominator[i] / denominator[0];
}
if (tempCoef1 != NULL)
free(tempCoef1);
if (tempCoef2 != NULL)
free(tempCoef2);
if (tempCoef3 != NULL)
free(tempCoef3);
}
/*======================================================================
* 函数名： pascalTriangle
* 函数功能：计算杨辉三角的第N行的值（数组），该系列值为(x+1)^N的系数，
* 加改进(x-1)^N的系数，最低次数在第一个
*
* 变量名称：
* N - 杨辉三角第N行，N=0,1,...,N
* symbol - 运算符号，0——(x+1)^N，1——(x-1)^N
* vector - 返回数组，杨辉三角的第N行的值
*
* 返回值： void
*=====================================================================*/
void pascalTriangle(
int N,
int symbol,
int *vector)
{
vector[0] = 1;
if(N == 0)
{
return;
}
else if (N == 1)
{
if(symbol == SYMBOL_ADD)
{
vector[1] = 1;
}
else
{
vector[0] = -1; //如果是减号，则第二项系数是-1
vector[1] = 1;
}
return;
}
int length = N + 1; //数组长度
int *temp = (int *)malloc(sizeof(int) * length);
/*
int temp[length]; //定义中间变量
*/
temp[0] = 1;
temp[1] = 1;
for(int i = 2; i <= N; i++)
{
vector[i] = 1;
for(int j = 1; j < i; j++)
{
vector[j] = temp[j - 1] + temp[j]; //x[m] = x[m-1][n-1] + x[m-1]
}
if(i == N) //最后一次不需要给中间变量赋值
{
if(symbol == SYMBOL_SUB) //运算符为减号
{
for(int k = 0; k < length; k++)
{
vector[k] = vector[k] * pow((float)-1, length - 1 - k);
}
}
return;
}
for(int j = 1; j <= i; j++)
{
temp[j] = vector[j];
}
}
if (temp != NULL)
free(temp);
}
/*======================================================================
* 函数名： coefficientEquation（整数）和coefficientEquation2（浮点数）
* 函数功能：计算多项式相乘的系数，最低次数在第一个
*
* 变量名称：
* originalCoef - 原来的系数数组，计算后的系数也存储在该数组内
* N - 原来数组中数据的长度，多项式最高次为N-1
* nextCoef - 与原数组相乘的数组的系数（两项）
*
* 返回值： void
*=====================================================================*/
void coefficientEquation(
int *originalCoef,
int N,
int *nextCoef,
int nextN)
{
double *tempCoef = (double *)malloc(sizeof(double) * (N+nextN-1));
/*
double tempCoef[N + nextN - 1]; //中间变量
*/
for(int i = 0; i < N + nextN - 1; i++)
{
tempCoef[i] = originalCoef[i]; //中间变量初始化
originalCoef[i] = 0;
}
for(int j = 0; j < nextN; j++)
{
for(int i = j; i < N + nextN - 1; i++)
{
originalCoef[i] += tempCoef[i-j] * nextCoef[j];
}
}
if (tempCoef != NULL)
free(tempCoef);
}
void coefficientEquation2(
double *originalCoef,
int N,
double *nextCoef,
int nextN)
{
double *tempCoef = (double *)malloc(sizeof(double) * (N+nextN-1));
/*
double tempCoef[N + nextN - 1]; //中间变量
*/
for(int i = 0; i < N + nextN - 1; i++)
{
tempCoef[i] = originalCoef[i]; //中间变量初始化
originalCoef[i] = 0;
}
for(int j = 0; j < nextN; j++)
{
for(int i = j; i < N + nextN - 1; i++)
{
originalCoef[i] += tempCoef[i-j] * nextCoef[j];
}
}
if (tempCoef != NULL)
free(tempCoef);
}
/*======================================================================
* 方法名： filter
* 方法功能：根据数字滤波器系统方法（系数），对原始信号进行滤波
*
* 变量名称：
* butterValue - 存放滤波器参数（阶数和截止频率）的结构体变量
* numerator - 系统方法，分子系数数组
* denominator - 系统方法，分母系数数组
* xVector - 输入的原始信号（数组）
* length - 原始信号的长度，也是滤波后信号的长度
* yVector - 滤波后的信号（数组）
*
* 返回值： 设计是否成功，true-成功，false-失败
*=====================================================================*/
BOOLEAN filter(
ButterFilterStruct butterValue,
double *numerator,
double *denominator,
double *xVector,
int length,
double *yVector)
{
BOOLEAN isFilterOK = FALSE;
if(!butterValue.isFOK)
{
printf("系统方法错误！");
isFilterOK = FALSE;
return isFilterOK;
}
if(butterValue.N > 10)
{
printf("失败！滤波器的阶数不能大于10。");
isFilterOK = FALSE;
return isFilterOK;
}
int N = butterValue.length; //系数数组的长度
//返回值初始化
for(int i = 0; i < length; i++)
{
yVector[i] = 0.0; //后面循环中用到y递归算法，需要提前初始化
}
//第一层循环，计算length个y的输出值
for(int i = 0; i < length; i++)
{
if(i == 0)
{
yVector[i] = numerator[i]*xVector[i];
}
else
{
yVector[i] = numerator[0]*xVector[i];
//第二层循环，计算每个y的每一项
for(int j = 1; j <= i && j < N; j++)
{
yVector[i] += numerator[j]*xVector[i-j] - denominator[j]*yVector[i-j];
}
}
yVector[i] /= denominator[0];
}
isFilterOK = TRUE;
return isFilterOK;
}
525
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
