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

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

注：我比较喜欢使用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处理结果对比。