关闭

数字信号处理公式变程序(四)——巴特沃斯滤波器(中)

标签: 算法ios开发巴特沃斯滤波器数字滤波器编程
3706人阅读 评论(0) 收藏 举报
分类:

上一篇写了巴特沃斯滤波器设计的所有理论知识,这篇文章就着手编程吧~

注:我比较喜欢使用matlab(也喜欢自己修改matlab的源码做测试,所以重装了好几次了,囧。。。),有部分参考matlab的算法进行(或者直接翻译为OC),算法的运行结果经常会跟matlab运算结果进行比较,算法只做学习用,转载请注明。

另外可能会有不足或者理解偏差的地方,路过的高人请不吝赐教。


闲谈少絮,开始!

========================================================

说明:所有流程图都是伪流程图,认真你就输了,嘿嘿。。。

具备了所有理论知识之后,编程变得简单粗暴。上一篇文章已经确定了开发的流程和思路,即:

①根据给定的滤波器的指标通带截止频率wp,阻带截止频率wst, 通带波纹δ1(Rp(dB)), 阻带波纹δ2(As(dB)),和滤波器的类型确定滤波器阶数N和3dB频率wc。

②根据N查表确定归一化滤波器的系数,并带入wc去归一化转化为H(s),求出其分式的系数数组a(N)和b(M);

③根据求得的系数和目标滤波器类型(低通、高通、带通、带阻),带入相应的s=f(z)转化公式(见上一章的思路图),求得数字滤波器的系统函数H(z)的分子分母系数数组。

④将原始信号和H(z)组成常系数线性差分方程,解出滤波后的信号(此为输出结果,可与matlab运行结果进行对比测试)。


代码分为三个模块:计算原型滤波器的参数(阶数和3dB频率)、计算数字滤波器的系统函数H(z)、滤波方法fliter。


一、计算原型滤波器的参数

首先需要输入参数,模仿matlab的实现,将传入参数定为通带截止频率、阻带截止频率、通带最大衰减、阻带最小衰减、抽样频率和滤波器的类型。流程图如下

注:传入的频率为模拟频率,单位为Hz,三种频率的关系如下表。

模拟频率     f  (Hz)                   Ω = 2π f

模拟角频率 Ω (rad/s)     =>      ω = ΩT

数字频率     ω (rad)                  ω = 2π f T

低通、高通部分求阶数和3dB的流程图如下图所示(频率不满足判断条件返回错误):


带通、带阻部分求阶数和3dB的流程图如下图所示(频率不满足判断条件返回错误):



本部分的代码如下,返回值是一个自定义的滤波器信息的结构体。

/*======================================================================
 * 方法名:  filterIIRButterLowpass
 * 方法功能:设计巴特沃斯样本低通示波器
 *
 * 变量名称:
 *          fpass - 通带截止频率(模拟频率)
 *          fstop - 阻带截止频率(模拟频率)
 *          rp    - 通带最大衰减(dB)
 *          rs    - 阻带最小衰减(dB)
 *          Fs    - 采样频率
 *
 * 返回值:  返回巴特沃斯低通滤波器的阶数N和截止频率Ws结构体
 *=====================================================================*/

+ (ButterFilterStruct)filterIIRButterLowpass:(float *)passF andStopF:(float *)stopF andPassRipple:(float)rp andStopRipple:(float)rs andFs:(float)fs andFilterType:(int)filterType
{
    ButterFilterStruct nAndFc;      //返回滤波器的阶数N和截止频率fc
    nAndFc.filterType = filterType; //滤波器类型
    float nOfN = 0.0;
    float passW = 0.0, stopW = 0.0, wa, wc; //wa = stopW/passW或其导数
    float passF1 = 0.0, passF2 = 0.0, stopF1 = 0.0, stopF2 = 0.0, w0 = 0.0;//w0 - 中心频率
    float 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;
                NSLog(@"错误!应满足:passF < stopF");
                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;
                NSLog(@"错误!应满足:passF > stopF");
                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;
                NSLog(@"错误!应满足: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;
                NSLog(@"错误!应满足: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;
}

结构体的定义为:

//巴特沃斯滤波器参数
typedef struct
{
    int N;          //巴特沃斯滤波器阶数
    int length;     //滤波器系统函数系数数组的长度
    float fc;       //巴特沃斯滤波器截止频率
    float cosW0;    //中心频率,带通带阻时用到
    float fs;       //采样频率
    int filterType; //需要设计的数字滤波器类型
    Boolean isFOK;
}ButterFilterStruct;

二、计算数字滤波器的系统函数H(z)


1.总框图

我习惯先说框架(整体的内容),然后分开讲解各部分的具体实现。滤波器设计的总框架流程图如下图所示:


2.各部分的具体实现

下面把总框图中的部分单元(模块)摘出来一一介绍。

(1)首先去归一化的流程图如下:


看上去流程图很复杂的样子,起始代码很简单。。。其中,g_butterPb是全局变量二维数组,存放归一化巴特沃斯系统函数的分母的系数数组。

/*======================================================================
 * 方法名:  butterSbValue
 * 方法功能:计算巴特沃斯滤波器分母多项式H(s)的系数Sb,注意:分子为Wc^N
 * 说明:   Sb[k] = Wc^(N-k) * Pb,其中Pb是归一化的分母多项式的根,可查表得到
 *         系数由低次向高次排列
 *
 * 变量名称:
 *          butterValue   - 存放滤波器参数(阶数和截止频率)的结构体变量
 *          returnSb      - 计算结果
 *
 * 返回值:  void
 *=====================================================================*/

+ (void)butterSbValue:(ButterFilterStruct)butterValue andReturnSb:(float *)returnSb
{
    int length = butterValue.N;        //滤波器阶数
    float 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
}
附上全局变量的值

//定义巴特沃斯滤波器pb系数列表(b0,b1,...,bn)
static float g_butterPb[10][10] = {{1.0,0,0,0,0,0,0,0,0,0},
    {1.0, 1.4142136, 0,0,0,0,0,0,0,0},
    {1.0, 2.0, 2.0, 0,0,0,0,0,0,0},
    {1.0, 2.6131259, 3.4142136, 2.6131259, 0,0,0,0,0,0},
    {1.0, 3.236068, 5.236068, 5.236068, 3.236068, 0,0,0,0,0},
    {1.0, 3.8637033, 7.4641016, 9.1416202, 7.4641016, 3.8637033, 0,0,0,0},
    {1.0, 4.4939592, 10.0978347, 14.5917939, 14.5917939, 10.0978347, 4.4939592, 0,0,0},
    {1.0, 5.1258309, 13.1370712, 21.8461510, 25.6883559, 21.8461510, 13.1370712, 5.1258309, 0,0},
    {1.0, 5.7587705, 16.5817187, 31.1634375, 41.9863857, 41.9863857, 31.1634375, 16.5817187, 5.7587705, 0},
    {1.0, 6.3924532, 20.4317291, 42.8020611, 64.8823963, 74.2334292, 64.8823963, 42.8020611, 20.4317291, 6.3924532}};


由于高通与低通的滤波器设计思路和设计代码近似(仅有变化公式不同),所以将二者代码合并为一个方法,中间用滤波器的类型参数区别。同理,带通和带阻滤波器也做相应的处理。因此在总框图中就出现“低高通处理”和“带通、带阻处理”的模块,下面先介绍低高通处理方法。

(2)高低通处理方法

高低通处理方法的流程图如下图所示:


代码如下

#pragma mark - 滤波器数字化s = f(z)


/*======================================================================
 * 方法名:  butterLowOrHigh
 * 方法功能:计算巴特沃斯低通(高通)滤波器系统方法的系数,包括分子和分母系数
 *
 * 变量名称:
 *          butterValue   - 存放滤波器参数(阶数和截止频率)的结构体变量
 *          sb            - 传入的模拟滤波器的系数,即H(s)的分母系数
 *          numerator     - 计算后的分子系数数组
 *          denominator   - 计算后的分母系数数组
 *
 * 返回值:  void
 *=====================================================================*/

+ (void)butterLowOrHigh:(ButterFilterStruct)butterValue andSb:(float *)sb andNumerator:(float *)numerator andDenominator:(float *)denominator
{
    int length = butterValue.N;    //滤波器阶数
    
    int tempCoef1[length + 1];     //定义系数数组,用于存放1 - z^(-1)、1 + z^(-1)每项次幂(0-N)系数,最低次在第一个
    int tempCoef2[length + 1];
    int otherN;                    //1+z^(-1)的次数
    
    float 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)
        {
            [MyMath pascalTriangle:i andSymbol:SYMBOL_SUB andReturnVector:tempCoef1];      //利用杨辉三角计算1 - z^(-1)幂的系数
            [MyMath pascalTriangle:otherN andSymbol:SYMBOL_ADD andReturnVector:tempCoef2]; //利用杨辉三角计算1 + z^(-1)幂的系数
        }
        else
        {
            [MyMath pascalTriangle:i andSymbol:SYMBOL_ADD andReturnVector:tempCoef1];      //利用杨辉三角计算1 + z^(-1)幂的系数
            [MyMath pascalTriangle:otherN andSymbol:SYMBOL_SUB andReturnVector:tempCoef2]; //利用杨辉三角计算1 - z^(-1)幂的系数
        }
        
        [MyMath coefficientEquation:tempCoef1 andOriginalN:i+1 andNextCoef:tempCoef2 andNextN: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];
    }
}

(3)带通、带阻处理方法

带通阻处理方法的流程图如下图所示:


代码如下

/*======================================================================
 * 方法名:  butterPassOrStop
 * 方法功能:计算巴特沃斯带通(带阻)滤波器系统方法的系数,包括分子和分母系数
 *
 * 变量名称:
 *          butterValue   - 存放滤波器参数(阶数和截止频率)的结构体变量
 *          sb            - 传入的模拟滤波器的系数,即H(s)的分母系数
 *          numerator     - 计算后的分子系数数组
 *          denominator   - 计算后的分母系数数组
 *
 * 返回值:  void
 *=====================================================================*/

+ (void)butterPassOrStop:(ButterFilterStruct)butterValue andSb:(float *)sb andNumerator:(float *)numerator andDenominator:(float *)denominator
{
    int length = butterValue.length;      //滤波器系数长度
    
    int tempCoef1[length];                //定义系数数组,用于存放1 - z^(-2)、1 - 2*cos(w0)*z^(-1) + z^(-2)每项次幂(0-N)系数,最低次在第一个
    float tempCoef2[length];
    float tempCoef3[length], tempCoef[3];
    int otherN;                           //1+z^(-1)的次数(pass),1 - 2*cos(w0)*z^(-1) + z^(-2)的次数(stop)
    
    float 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) //带通滤波器
    {
        [MyMath pascalTriangle:butterValue.N andSymbol:SYMBOL_SUB andReturnVector: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++)
        {
            [MyMath coefficientEquation2:tempCoef3 andOriginalN:j*2+1 andNextCoef:tempCoef andNextN: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;
        }
        
        [MyMath pascalTriangle:otherN andSymbol:SYMBOL_SUB andReturnVector: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++)
        {
            [MyMath coefficientEquation2:tempCoef3 andOriginalN:j*2+1 andNextCoef:tempCoef andNextN:3];
        }
        
        [MyMath coefficientEquation2:tempCoef3 andOriginalN:(butterValue.N - otherN)*2+1 andNextCoef:tempCoef2 andNextN: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];
    }
}

(4)辅助程序(工具)

在(2)(3)流程图中用到了杨辉三角和多项式相乘的方法,这两个方法有助于多项式的展开时系数的计算,流程图就不贴了,说明如下

①杨辉三角是一个数列,当一个表达式是类似(s+a)(s+a)(s+a)表达式时,其展开结果就是杨辉三角的一行数。当符号为减时,只需要吧a换成(-a)即可。因此,需要传递一下分式的符号作为参数,进行计算。

②当表达式是两个任意的式子时,需要将表达式写成幂连续下降的形式(如a0+a1*s+0*s^2+0*s^3+a4*s^4),当中间的某项没有时,说明其系数为0。

两部分的代码如下:首先是杨辉三角

/*======================================================================
 * 函数名:  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 andSymbol:(int)symbol andReturnVector:(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[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][n] = x[m-1][n-1] + x[m-1][n]
        }
        if(i == N) //最后一次不需要给中间变量赋值
        {
            if(symbol == SYMBOL_SUB) //运算符为减号
            {
                for(int k = 0; k < length; k++)
                {
                    vector[k] = vector[k] * pow(-1, length - 1 - k);
                }
            }
            return;
        }
        for(int j = 1; j <= i; j++)
        {
            temp[j] = vector[j];
        }
    }
}
然后是多项式相乘:

/*======================================================================
 * 函数名:  coefficientEquation(整数)和coefficientEquation2(浮点数)
 * 函数功能:计算多项式相乘的系数,最低次数在第一个
 *
 * 变量名称:
 *          originalCoef - 原来的系数数组,计算后的系数也存储在该数组内
 *          N            - 原来数组中数据的长度,多项式最高次为N-1
 *          nextCoef     - 与原数组相乘的数组的系数(两项)
 *
 * 返回值:  void
 *=====================================================================*/

+ (void)coefficientEquation:(int *)originalCoef andOriginalN:(int)N andNextCoef:(int *)nextCoef andNextN:(int)nextN
{
    float 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];
        }
    }
}
+(void)coefficientEquation2:(float *)originalCoef andOriginalN:(int)N andNextCoef:(float *)nextCoef andNextN:(int)nextN
{
    float 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];
        }
    }
}

现在已经把数字滤波器的系统函数(其实是系统函数的分子分母的系数)求出来了,但是只有系统函数不能处理信号也没有什么用,也就是说对于滤波还差一步:对信号的处理——fliter方法。


累残了,fliter部分放在下一节吧。。。

=============================================================

OK,具体求滤波器系统函数的代码随笔就此结束,下一节讲fliter方法和与matlab处理结果对比。

2
0

查看评论
* 以上用户言论只代表其个人观点,不代表CSDN网站的观点或立场
    个人资料
    • 访问:65744次
    • 积分:883
    • 等级:
    • 排名:千里之外
    • 原创:20篇
    • 转载:8篇
    • 译文:0篇
    • 评论:29条
    文章分类
    最新评论