"); //-->
总体说来,定点数的优点在于运行速度比浮点数快,缺点是设计时比较麻烦,要通过malab在stm32上实现IIR的定点滤波程序,需要按照以下步骤来实现。
1.通过仿真确定内部状态值的范围,范围定的误差越小,精度越高,但容易出现饱和现象。
首先使用fdatool创建一个IIR滤波器,这里我们仍然设计一个2阶高通滤波,类型选切比雪夫I型,截止频率300Hz,量化方法采用浮点数,然后将该滤波器实例化,生成仿真模型,在仿真模型中添加一个正弦信号发生器、一个常量、一个示波器,然后将信号线连接好,如下图所示。
双击图中的sine wave模块,弹出设置对话框,由于我们的stm32的ADC输入电压范围在0~3.3V,输入信号是一个带直流偏置的正弦波信号,为了获得最大的摆幅,将直流信号设置到12位ADC的中间值2048,交流信号的幅值设置到2047,这样就可以满足ADC采集到最大摆幅值。因此Amplitude输入2047,频率先设置为1KHz。
点击OK,返回仿真主画面,进入Filter模块,添加一个示波器,用于观察各状态值,如图所示。
完成以上设置之后,再点击仿真运行,仿真完成之后可以看到内部各状态的最大值,由于输入的频率会影响内部状态值,因此需要将输入正弦波的频率从0改变到4KHz,这里只需要输入0、10、1000几个典型值即可,由于定点数表示的范围都是2的n次幂,因此通过观察得到状态值大概的范围如下。
在清楚了各状态的范围之后,就可以设置matlab量化参数了,回到fdatool界面,点击图中红色圈标记的图标,这个图标是滤波器数值量化方法。
点击后选择图中红色的”Fixed-point”,可以看到右侧金黄色圈住的三个选项,分别是Coefficients、Input/Output、FilterInternals,其中Coefficients表示IIR系数的定标,Input/Output表示输入、输出以及段输入输出定标,Filter Internals表示内部状态的定标。在c语言程序中IIR系数为16位有符号数,因此Coefficient word lenght 填入16根据分子、分母以及缩放因子值得大小,选择合适的定标位数,如分母的值是1、-2、1,分子的值是1、-1.7465710178967648、0.79152631529839501,因此分母和分子小数点位数填入14(s16,14的范围为[-2 2)),缩放因子值为0.039834401095698677和19.790259679335225,因此缩放值得小数点位填入10(s16,10的范围为[-32 32))。注:在生成的c头文件中会发现matlab计算出的分子、分母定标是按照缩放因子小数位定标的。
完成系数的定标之后,点击Input/Output选项,设置输入输出的定标。如图所示,我们已经清楚了输入输出、段输入输出的数值范围,根据这一范围设定各个定标值。其中输入范围[0,4096),但我们采用有符号数计算,因此实际输入范围扩大到[-4096 4096),因此将输入定标为3;同理输出范围[-2048 2048),因此输出定标4,段输入定标为23,段输出定标为24。
完成了输入输出的定标,接下来设置内部状态值的定标,如下图所示。根据各参数值的具体的范围,设置对应的定标数。由于和积部分采用32位数,因此在Product word length以及Accum word length两项填入32。
到此完成了定点数量化定标的工作,最后将每个参数的定标值再列在下面的图中。注意途中红色圈住的三个参数与前面定标设置菜单中填入的”s16,14”不同,在生成的c头文件中实际是按照缩放因子的定标计算的。
2.修改官方程序。官方给定程序只是演示,具体定标需要自己操作。
再修改程序之前,先看一下matlab导出的c头文件(注:由于IIR系数设置为16位有符号数,故在导出头文件时选择数据类型为Signed 16-bit integer,如下图)的注释。
/* 头文件注释
* Discrete-Time IIR Filter (real)
* -------------------------------
* Filter Structure : Direct-Form II, Second-Order Sections直接II型
* Number of Sections : 1 定义一段二阶
* Stable : Yes
* Linear Phase : No
* Arithmetic : fixed
* Numerator : s16,14 -> [-2 2) IIR系数分母定标
* Denominator : s16,14 -> [-2 2) IIR系数分子定标
* Scale Values : s16,10 -> [-32 32) 伸缩因子定标
* Input : s16,3 -> [-4096 4096) 输入定标
* Section Input : s32,23 -> [-256 256) 段输入定标
* Section Output : s32,24 -> [-128 128) 段输出定标
* Output : s16,4 -> [-2048 2048) 输出定标
* State : s32,19 -> [-4096 4096) 状态值定标
* Numerator Prod : s32,17 -> [-16384 16384) 分母系数积的定标
* Denominator Prod : s32,17 -> [-16384 16384) 分子系数积的定标
* Numerator Accum : s32,17 -> [-16384 16384) 分母系数和的定标
* Denominator Accum : s32,17 -> [-16384 16384) 分子系数和的定标
* Round Mode : convergent
* Overflow Mode : wrap
* Cast Before Sum : true
*/
#define MWSPT_NSEC 3
const int NL[MWSPT_NSEC] = { 1,3,1 };
const int16_T NUM[MWSPT_NSEC][3] = {
{
41, 0, 0
},
{
1024, -2048, 1024
},
{
20265, 0, 0
}
};
const int DL[MWSPT_NSEC] = { 1,3,1 };
const int16_T DEN[MWSPT_NSEC][3] = {
{
1024, 0, 0
},
{
1024, -1789, 811
},
{
1024, 0, 0
}
};
结合matlab生成的直接II型滤波器的内部结构,可以发现这些定标值是与定点设置中的量化小数值相同的。根据分子与分母的值(头文件中红色的值有效),可以定义IIR系数:
int16_t IIRCoeff[5] = {41,-1789,811,-2048,20265};
再看官方给定的程序,如下所示,这里只列举了一级二阶滤波。
void iir_biquad_stm32(u16 *y, u16 *x, int16_t *IIRCoeff, u16 ny)
{
s32 w1_2 = 0, w1_1 = 0, w1=0;
for (i=0; i<ny-2; i++)
{
w1 = x[i+2] - IIRCoeff[0]*w1_1 - IIRCoeff[1]*w1_2;
y[i+2] = IIRCoeff[2]*w1 + IIRCoeff[3]*w1_1 + IIRCoeff[3]*w1_2;
w1_2 = w1_1;
w1_1 = w1;
}
}
根据官方代码,转换为结构流图应该是这样的:
为了配合matlab生成的IIR系数,需要修改官方程序,这里我们仍然采用ADC采集一次数据处理一次的办法,所以不需要输入数组,另外由于是高通滤波,输出y具有正负数特性,因此输出y应该改为int类型,修改后的代码如下:
void iir_biquad_stm32(int16_t *y, u16 *x, int16_t *IIRCoeff, u16 ny)
{
static s32 w1_2 = 0, w1_1 = 0, w1=0;
w1 = IIRCoeff[0]*x[0]- IIRCoeff[1]*w1_1 - IIRCoeff[2]*w1_2;
y[0] = IIRCoeff[4]*(w1+ IIRCoeff[3]*w1_1 + w1_2);
w1_2 = w1_1;
w1_1 = w1;
}
测试该定点数程序,发现输出波形与我们预想的结果完全不同!甚至是一个很乱的波形,这就涉及到了定点数定标计算方法。再简单的回顾一下定点数计算规则:两定点数相乘,得到结果的定标是两乘数定标值相加;两定点数相加减,前提是两数的定标值相等。用数学表达式表达出来就是:
x*y=z,其中x定标a,y定标b,那么z定标为a+b;
x-y=z,其中x定标a,那么y定标也必须为a,若两者定标不同,则要通过移位来让两者定标相等,最后z定标也为a。因此代码应该这样修改:
void iir_biquad_stm32(int16_t *y, u16 *x, int16_t *IIRCoeff, u16 ny)
{
static s32 w1_2 = 0, w1_1 = 0, w1=0;
s32 ytmp; //计算y[0]时使用的临时变量,防止编译器优化造成的精度损失。
w1 = ((s32)((IIRCoeff[0]*x[0])<<3) - (s32)(IIRCoeff[1]*((s16)(w1_1>>16))) - (s32)(IIRCoeff[2]*((s16)(w1_2>>16))))<<6;
ytmp = (w1>>6) + IIRCoeff[3]*((s16)(w1_1>>16)) + (w1_2>>6);
ytmp = ytmp>>3; //
ytmp = IIRCoeff[4]*ytmp; //
ytmp = ytmp>>20;
y[0] = ytmp;
w1_2 = w1_1;
w1_1 = w1;
}
w1的计算过程
看起来是不是有点乱?不急,我们先看w1的计算,这个计算对了,输出也就好计算了。
式中”w1 = ((s32)((IIRCoeff[0]*x[0])<<3) - (s32)(IIRCoeff[1]*((s16)(w1_1>>16))) - (s32)(IIRCoeff[2]*((s16)(w1_2>>16))))<<6“,由于w1定标为19位,所以等号后面的部分也应该按照定标19位来计算。
第一部分”((s32)((IIRCoeff[0]*x[0])<<3)”,由于x[0]定标3,因此要左移3位,这样第一部分的定标为10+3=13。
第二部分“(s32)(IIRCoeff[1]*((s16)(w1_1>>16)))”,这里先将w1_1右移16位,是为了让16位数IIRCoeff与右移后的w1_1相乘得到一个32位数,否则将会溢出。右移之后再将w1_1转化为16位有符号数,然后再与IIRCoeff相乘之后得到一个定标为13(10+19-16)的定点数。
第三部分“(s32)(IIRCoeff[2]*((s16)(w1_2>>16)))”,转化方法与第二部分相同。
最终这三个部分都是13位的定点数,因此可以直接相加减,得到一个13位定点数,由于整个式子等号的左边w1_1为19位定点数,因此等号右边的式子需要整体再左移6(19-13)位。
y[0] 的计算过程
由于y[0]是一个16位有符号数,而计算过程中的w1、w1_1、w1_2为32位有符号数,为了提高计算精度,这里引进一个临时变量ytmp(32位有符号数)。
先看”ytmp = (w1>>6) + IIRCoeff[3]*((s16)(w1_1>>16)) + (w1_2>>6)“,第一部分“(w1>>6)”与第三部分“(w1_2>>6)”右移的位数都是为了配合第二部分的定标数,因为第二部分的计算方法与w1的第二三部分相同,都是为了得到一个32位有符号数。因此ytmp定标为13。
再看下面几个式子:
ytmp = ytmp>>3; //q10
ytmp = IIRCoeff[4]*ytmp; //q20
ytmp = ytmp>>20; //q0
y[0] = ytmp;
为什么ytmp 要先右移3位,与IIRCoeff[4]相乘之后又右移20位,这里需要倒着推回去。虽然y[0]定标为Q4,但我们最后在使用这个变量时需要定标为Q0,所以最终输出y[0]的范围是[-2048 2048)。算法的最后一步由ytmp给y[0]赋值,为了达到最高运算精度,ytmp的十六进制数最大值是0x7fffffff,而这个值右移20位正好是2047,因此最后ytmp右移20位。
再看“ytmp = IIRCoeff[4]*ytmp”,因最终右移20位的ytmp定标为Q0,所以这个式子左边的ytmp定标为Q20,而IIRCoeff[4]定标Q10,所以右边的ytmp定标为10(20-10),再倒推”ytmp = ytmp>>3“,因左边的ytmp定标为Q13,为了与后一步的定标Q10吻合,因此右移3位。
到此完成了y[0]的计算过程。
2.测试
为了简便,程序使用一级二阶直接II型IIR滤波来做测试。测试的目的一是为了检验各种输入情况下都不会造成内部状态的溢出,不同于浮点数计算,定点数首先是根据运算范围来定的小数点标,因此在不同输入下有可能造成内部溢出,这样就不能达到预期滤波的效果。但只要在前期matlab 仿真时对内部状态范围有很清楚的掌握,那么在编写程序时就不会造成数据溢出。
为了对比方便,将定点数和浮点数放在一张图来比较。
f = 50
f = 100
f = 200
f = 300
f = 500
f = 800
通过对比,直观上我们发现定点IIR滤波和浮点IIR滤波效果基本上是一致的,不仅是幅度的大小还是相位的变化,两者几乎没有太大的差别。
另一方面,使用定点数的目的是为了在没有FPU的单片机上加快数据的处理速度,还需要对比使用浮点计算的IIR滤波程序的耗时。根据测试结果,采用一级二阶的IIR滤波的程序,用c语言定点数的耗时为1.14us,浮点数耗时为9.84us,定点数计算速度明显快于浮点数,但定点精度不及浮点。
本人水平有限,若有错误之处,敬请指出,谢谢!转自www.mcukey.com
*博客内容为网友个人发布,仅代表博主个人观点,如有侵权请联系工作人员删除。