第7章应用程序设计.ppt

《第7章应用程序设计.ppt》由会员分享，可在线阅读，更多相关《第7章应用程序设计.ppt（86页珍藏版）》请在三一办公上搜索。

1、第7章 应用程序设计,本章内容提要：,定标与溢出处理基础算术运算FIR滤波器IIR滤波器快速傅里叶变换(FFT),7.1定标与溢出处理,数的定标溢出的处理方法常用信号处理算法中的定标方法,7.1.1 数的定标,小数定标的概念设定一个16位数的小数点处于该数中的哪一位可以表示不同大小和不同精度的小数Q表示法表7-1列出了一个16位数的16种Q表示及它们所能表示的十进制数值范围,表7-1 Q表示及数值范围,不同的Q所表示的数不仅范围不同，而且精度也不相同Q越大，数值范围越小，但精度越高Q越小，数值范围越大，但精度就越低例如，Q0的数值范围是-32768到+32767，其精度为1；Q15的数值范围为

2、-1到0.9999695，精度为 1/32768=0.00003051,对定点数而言，数值范围与精度是一对矛盾。一个变量要想能够表示比较大的数值范围，必须以牺牲精度为代价；而想提高精度，则数的表示范围就相应地减小。在实际的定点算法中，应该根据具体问题进行折衷处理,以达到最佳效果。,在C55x中，16位整数采用补码形式表示。每个采用Qi定标的16位数用1个符号位、i个小数位和15-i个整数位来表示。,表7-2 同样的数在不同定标方式下所表示的具体数值,同样一个16位数，若小数点设定的位置不同，它所表示的数也就不同。,7.1.2 溢出的处理方法,如果算术运算结果超出寄存器所能表示的最大数就会出现溢

3、出因为16位定点DSP的动态范围有限，所以在使用时必须注意动态范围以防溢出溢出还与输入信号的特性和运算法则有关,1.溢出,C55x有以下几种硬件特性可以处理溢出：保护位 C55x的每个累加器都有8个保护位（3932位），允许连续256次乘加操作而累加器不溢出溢出标志位 C55x的每个累加器都有相关的溢出标志位，当累加器操作结果出现溢出时，这个标志位就会置位,2.C55x的溢出处理机制,饱和方式位SATD和SATASATD控制D单元的操作，SATA控制A单元的操作。如果SATD=1，当D单元发生溢出时，对D单元的结果进行饱和处理。不管饱和方式位的值是什么，当累加器发生溢出时，相应的溢出标志位都会

4、被置位A单元没有溢出标志位，但如果SATA=1，发生溢出时，结果也会进行饱和处理,饱和处理是用最近的边界值代替溢出结果。例如，16位寄存器的范围是8000h（最小负数）7FFFh（最大正数），饱和处理就是用7FFFh代替比7FFFh大的结果；用8000h代替比8000h小的结果。,饱和。饱和是一种处理溢出的方法，但是饱和会剪掉部分输出信号，可能会引起信号失真和引起系统非线性。输入定标。分析所要使用的系统，假定最坏的情况，然后对输入信号定标，以防止溢出。但是这种方法会极大地降低输出信号的精确度。固定定标。假定最坏的情况，对中间结果定标。这种方法可以防止溢出，同时增加了系统的信噪比。动态定标。可以

5、监测中间结果的范围，只在需要的时候对中间结果定标。这种方法可以防止溢出但会增加计算量。,3.溢出的处理方法,7.1.3 常用信号处理 算法中的定标方法,FIR滤波器的定标方法在FIR滤波器中处理溢出的最好方法是设计时使滤波器的增益小于1，这样就不需要对输入信号定标。这种方法和累加器的保护位结合起来，可以有效地防止溢出。由于对信号处理的负面影响，在FIR滤波器中不使用固定定标和输入定标。如果不考虑计算量，在FIR滤波器中可以使用动态定标。对一些类型的音频信号，饱和处理也是一种常用的方法。,IIR滤波器的定标方法IIR滤波器的定点实现推荐使用多个二阶基本节级联组成，这样可以减小高阶滤波器频率响应灵

6、敏度。由于滤波器系数的量化引入误差，因此避免溢出对IIR滤波器非常重要。可以通过把中间结果保存在处理器累加器来避免节间数据溢出。为防止在第k阶内部发生数据溢出，需要用增益系数对滤波器的单位脉冲响应（前馈通道）定标。动态标定方法。在每个阶段滤波器内部状态都被减半，以提高指令周期换取为代价提高了结果的精度。,FFT的定标方法,在FFT操作里，每次蝶形运算后数据平均增加一位。输入定标需要移位（FFT长度为N），这会导致在计算FFT之前就衰减 6 dB。在固定定标中，每级蝶形运算输出除以2，这是最常用的FFT定标方法，因为它简单而且有比较好的信噪比。但是，对于大的FFT，这种定标可能会使信息丢失。另一

7、种方法是动态定标，即在输出溢出时再除以2。在这种情况下，会在这个过程中指定一个变量，每定标一次变量的值加1，计算结束后根据变量的值把结果乘以一个系数。动态定标的信噪比最好，但会增加FFT循环次数。,7.2 基础算术运算,加减运算乘法运算除法运算小数乘法,7.2.1 加减运算,在数字信号处理中，加减运算是常见的算术运算。一般使用16位或32位加减运算，数值分析、浮点运算和其它操作可能需要32位以上的运算。C55x有直接完成16位或32位加减运算的指令，但没有能直接完成多字加减运算的指令。要进行多字加减运算，需要通过编程方法实现。,以下指令可在单周期内完成32位加法运算：MOV40 dbl(Lme

8、m),ACxADD dbl(Lmem),ACx 64位的高32位加法要考虑低32位加法产生的进位，使用以下指令：ADD uns(Smem),CARRY,ACx 以下指令可在单周期内完成32位减法运算：MOV40 dbl(Lmem),ACxSUB dbl(Lmem),ACx 64位的高32位减法要考虑低32位减法产生的借位，使用以下指令：SUB uns(Smem),BORROW,ACx,例7-1，64位加法运算。文件名为：add64.asm。,.mmregs.model call=c55_std.model mem=large;*;64位加法 指针分配;X3 X2 X1 X0 AR1-X3(偶地

9、址);+Y3 Y2 Y1 Y0 X2;-X1;W3 W2 W1 W0 X0;AR2-Y3(偶地址);Y2;Y1;Y0;AR3-W3(偶地址);W2;W1;W0;*,.sect.text.align 4.globalstart.symstart,start,36,2,0start:MOV#0100h,AR1 MOV#0104h,AR2 MOV#0108h,AR3L1:MOV40 dbl(*AR1(#2),AC0;AC0=X1 X0 ADD dbl(*AR2(#2),AC0;AC0=X1 X0+Y1 Y0 MOV AC0,dbl(*AR3(#2);保存W1 W0.MOV40 dbl(*AR1),A

10、C0;AC0=X3 X2 ADD uns(*AR2(#1),CARRY,AC0;AC0=X3 X2+00 Y2+CARRY ADD*AR2#16,AC0;AC0=X3 X2+Y3 Y2+CARRY MOV AC0,dbl(*AR3);保存 W3 W2.B L1,例7-2，64位减法运算程序。文件名为：sub64.asm。,.mmregs.model call=c55_std.model mem=large;*;64位减法 指针分配;X3 X2 X1 X0 AR1-X3(偶地址);Y3 Y2 Y1 Y0 X2;-X1;W3 W2 W1 W0 X0;AR2-Y3(偶地址);Y2;Y1;Y0;AR3

11、-W3(偶地址);W2;W1;W0;*,.sect.text.align 4.globalstart.symstart,start,36,2,0start:MOV#0100h,AR1 MOV#0104h,AR2 MOV#0108h,AR3L1:MOV40 dbl(*AR1(#2),AC0;AC0=X1X0 SUB dbl(*AR2(#2),AC0;AC0=X1X0-Y1Y0 MOV AC0,dbl(*AR3(#2);保存W1W0.MOV40 dbl(*AR1),AC0;AC0=X3X2 SUB uns(*AR2(#1),BORROW,AC0;AC0=X3X2-00Y2-BORROW SUB*A

12、R2#16,AC0;AC0=X3X2-Y3Y2-BORROW MOV AC0,dbl(*AR3);保存 W3W2.B L1,7.2.2乘法运算,C55x提供了硬件乘法器，16位乘法可在一个指令周期内完成。高于16位的乘法运算可以采用下述方法实现（以32位乘法为例）。,例7-3,32位整数乘法运算。文件名:mpy32.asm,.mmregs.model call=c55_std.model mem=large;*;本子程序是两个32位整数乘法，得到一个64位结果。操作数取自数;据存储器，运算结果送回数据存储器。;数据存储:指针分配:;X1 X0 32位操作数 AR0-X1;Y1 Y0 32位操作

13、数 X0;W3 W2 W1 W0 64位结果 AR1-Y1;Y0;入口条件:AR2-W0;SXMD=1(允许符号扩展)W1;SATD=0(不做饱和处理)W2;FRCT=0(关小数模式)W3;限制条件：延迟链和输入序列必须指定为长字类型。;*,.sect.text.align 4.global start.symstart,start,36,2,0start:MOV#0100h,AR0 MOV#0102h,AR1 MOV#0104h,AR2 BSET SXMD BCLR SATD BCLR FRCTL1:AMAR*AR0+;AR0指向X0|AMAR*AR1+;AR1指向Y0MPYM uns(*A

14、R0-),uns(*AR1),AC0;ACO=X0*Y0MOV AC0,*AR2+;保存W0MACM*AR0+,uns(*AR1-),AC0#16,AC0;AC0=X0*Y016+X1*Y0MACM uns(*AR0-),*AR1,AC0;AC0=X0*Y016+X1*Y0+X0*Y1MOV AC0,*AR2+;保存W1MACM*AR0,*AR1,AC0#16,AC0;AC0=AC016+X1*Y1MOV AC0,*AR2+;保存W2MOV HI(AC0),*AR2;保存W3 B L1,7.2.3 除法运算,C55x没有提供硬件除法器，也没有提供专门的除法指令，要实现除法运算需借助于条件减法指

15、令SUBC和重复指令RPT。根据被除数绝对值与除数绝对值的大小关系，除法的实现过程略有不同：当|被除数|除数|时，商为小数。当|被除数|除数|时，商为整数。需要注意的是：SUBC指令要求被除数和除数都必须为正。下面举例说明如何在C55x DSP中实现除法运算。,例7-4，无符号16位除16位整数除法。文件名为：udiv16o16.asm。,.mmregs.model call=c55_std.model mem=large;*;指针分配;AR0-被除数;AR1-除数;AR2-商;AR3-余数;注：;无符号除法，被除数、除数均为16位;关闭符号扩展，被除数、除数均为正数;运算完成后AC0(15-

16、0)为商，AC0(31-16)为余数;*,.sect.text.align 4.globalstart.symstart,start,36,2,0start:MOV#0100h,AR0 MOV#0101h,AR1 MOV#0102h,AR2 MOV#0103h,AR3L1:BCLR SXMD;清零SXMD(关闭符号扩展)MOV*AR0,AC0;把被除数放入AC0RPT#(16-1);执行subc 16次SUBC*AR1,AC0,AC0;AR1指向除数MOV AC0,*AR2;保存商MOV HI(AC0),*AR3;保存余数 B L1,例7-5，无符号32位除16位整数除法。文件名为：udiv3

17、2o16.asm。,.mmregs.model call=c55_std.model mem=large;*;指针分配;AR0-被除数高位;被除数低位;AR1-除数;AR2-商高位;商低位;AR3-余数;注：;无符号除法，被除数为32位，除数为16位;关闭符号扩展，被除数、除数均为正数;第一次除法之前，把被除数高位存入AC0;第一次除法之后，把商的高位存入AC0(15-0);第二次除法之前，把被除数低位存入AC0;第二次除法之后，AC0(15-0)为商的低位，AC0(31-16)为余数;*,.sect.text.align 4.globalstart.symstart,start,36,2,0

18、start:MOV#0100h,AR0 MOV#0102h,AR1 MOV#0104h,AR2 MOV#0106h,AR3L1:BCLR SXMD;清零SXMD(关闭符号扩展)MOV*AR0+,AC0;把被除数高位存入AC0|RPT#(15-1);执行subc 15次SUBC*AR1,AC0,AC0;AR1指向除数SUBC*AR1,AC0,AC0;执行subc最后一次|MOV#8,AR4;把AC0_L 存储地址装入AR4MOV AC0,*AR2+;保存商的高位MOV*AR0+,*AR4;把被除数低位装入AC0_LRPT#(16-1);执行subc 16次SUBC*AR1,AC0,AC0MOV

19、AC0,*AR2+;保存商的低位MOV HI(AC0),*AR3;保存余数BSET SXMD;置位SXMD(打开符号扩展)B L1,例7-6，带符号16位除16位整数除法。文件名为：sdiv16o16.asm。,.mmregs.model call=c55_std.model mem=large;*;指针分配;AR0-被除数;AR1-除数;AR2-商;AR3-余数;注：;带符号除法，被除数为16位，除数为16位;打开符号扩展，被除数、除数可为负数;除法运算之前，商的符号存入AC0;除法运算之后，商存入AC1(15-0)，余数存入AC1(31-16);*,.sect.text.align 4.g

20、lobalstart.symstart,start,36,2,0start:MOV#0100h,AR0 MOV#0101h,AR1 MOV#0102h,AR2 MOV#0103h,AR3L1:BSET SXMD;置位SXMD(打开符号扩展)MPYM*AR0,*AR1,AC0;计算期望得到的商的符号MOV*AR1,AC1;把除数存入AC1ABS AC1,AC1;求绝对值，|除数|MOV AC1,*AR2;暂时保存|除数|MOV*AR0,AC1;把被除数存入 AC1ABS AC1,AC1;求绝对值，|被除数|RPT#(16-1);执行subc 16次,SUBC*AR2,AC1,AC1;AR2-|除

21、数|MOV HI(AC1),*AR3;保存余数MOV AC1,*AR2;保存商SFTS AC1,#16;对商移位:把符号位放在最高位NEG AC1,AC1;对商求反XCCPART label,AC0#0;如果商的符号位为负,MOV HI(AC1),*AR2;用商的负值替换原来的商label:B L1,例7-7，带符号32位除16位整数除法。文件名为：sdiv32o16.asm,.mmregs.model call=c55_std.model mem=large;*;指针分配:(被除数和商都被指定为长字);AR0-被除数高半部分(NumH)(偶地址);被除数高半部分(NumL);AR1-除数(D

22、en);AR2-商的高半部分(QuotH)(偶地址);商的低半部分(QuotL);AR3-余数(Rem);注：;带符号除法，被除数为32位，除数为16位;打开符号扩展，被除数、除数可为负数;除法运算之前，期望的商的符号存入AC0;第一次除法运算之前，把被除数的高半部分存入AC1;第一次除法运算之后，把商的高半部分存入AC1(15-0);第二次除法运算之前，把被除数的低半部分存入AC1;第二次除法运算之后，把商的低半部分存入AC1(15-0),余数存入AC1(31-16);*,.sect.text.align 4.global start.symstart,start,36,2,0start:M

23、OV#0100h,AR0 MOV#0102h,AR1 MOV#0104h,AR2 MOV#0106h,AR3 MOV#0108h,AR4L1:BSET SXMD;置位SXMD(打开符号扩展)MPYM*AR0,*AR1,AC0;除法结果的符号位(NumH x Den)MOV*AR1,AC1;AC1=DenABS AC1,AC1;AC1=abs(Den)MOV AC1,*AR3;Rem=abs(Den)MOV40 dbl(*AR0),AC1;AC1=NumH NumLABS AC1,AC1;AC1=abs(Num)MOV AC1,dbl(*AR2);QuotH=abs(NumH);QuotL=ab

24、s(NumL),MOV*AR2,AC1;AC1=QuotHRPT#(15-1);执行subc 15次SUBC*AR3,AC1,AC1SUBC*AR3,AC1,AC1;最后一次执行subc|MOV#11,AR4;把AC1_L存储地址装入AR4MOV AC1,*AR2+;保存 QuotHMOV*AR2,*AR4;AC1_L=QuotHRPT#(16-1);执行subc 16次SUBC*AR3,AC1,AC1MOV AC1,*AR2-;保存 QuotLMOV HI(AC1),*AR3;保存 RemBCC skip,AC0=#0;如果实际结果应该为正数,跳到skip.MOV40 dbl(*AR2),A

25、C1;否则,对商取反.NEG AC1,AC1MOV AC1,dbl(*AR2)skip:B L1,7.2.4 小数乘法,在定点DSP的某些应用中，整数运算很难满足要求。这是因为它自身存在缺陷：两个16位整数相乘，乘积总是“向左增长”（即小数点左侧的位数增加），这意味着多次相乘后，乘积将很快超出定点器件的数据范围。保存32位乘积到存储器，要占用2个CPU周期和2个字的存储器空间。由于乘法器都是16位相乘，因此将32位乘积再作为乘法器的输入时就显得较繁琐，不能胜任递归运算。,为了克服这些缺陷，在实际应用中更多采用的是小数运算。小数运算具有如下优点：乘积总是“向右增长”。这就意味着超出定点器件数据范

26、围的将是不太感兴趣的部分。既可以存储32位乘积，也可以近存储高16位乘积，这就允许用较少的资源保存结果。可以用于递归运算。,例7-8，两个Q31格式有符号小数相乘，得到一个Q31格式结果。文件名为：mpyQ31.asm。,需要注意的是：两个带符号小数相乘，所得乘积带有2个符号位。为了解决冗余符号位的问题，需要在程序中设定状态寄存器ST1中的FRCT（小数方式）为1，这样当乘法器将结果传送至累加器时就会自动左移1位。,.mmregs.model call=c55_std.model mem=large;*;操作数取自数据存储器，运算结果送回数据存储器。;数据存储:指针分配:;X1 X0 Q31操

27、作数 AR0-X1;Y1 Y0 Q31操作数 X0;W1 W0 Q31结果 AR1-Y1;Y0;入口条件:AR2-W0（偶地址）;SXMD=1(允许符号扩展)W1;SATD=0(不做饱和处理);FRCT=1(运算结果左移一位);限制条件：W1被指定为偶地址;*,.sect.text.align 4.globalstart.symstart,start,36,2,0start:MOV#0100h,AR0 MOV#0102h,AR1 MOV#0104h,AR2 BSET SXMD BCLR SATD BSET FRCTL1:AMAR*AR0+;AR0指向X0 MPYM uns(*AR0-),*AR

28、1+,AC0;AC0=X0*Y1 MACM*AR0,uns(*AR1-),AC0;AC0=X0*Y1+X1*Y0 MACM*AR0,*AR1,AC0#16,AC0;AC0=AC016+X1*Y1 MOV AC0,dbl(*AR2);保存W1W0 B L1,7.3 FIR滤波器,数字滤波器是DSP的基本应用，有2种基本类型：有限冲激响应滤波器FIR无限冲激滤波器IIR一般来说，如果需要线性相位则选择用FIR滤波器，对于相位要求不敏感的场合可以选用IIR滤波器。本节主要讨论FIR滤波器的DSP实现方法，有关IIR滤波器的实现将在下一节中介绍。,7.3.1 FIR滤波器的基本结构,一个FIR滤波器的

29、输出序列和输入序列之间的关系，满足差分方程：传递函数为,FIR滤波器的结构：FIR滤波器的单位冲激响应是一个有限长序列。若为实数，且满足偶对称或奇对称的条件，则FIR滤波器具有线性相位特性。偶对称线性相位FIR滤波器的差分方程为：,7.3.2 FIR滤波器的C语言编程实现,例7-9，直接型FIR滤波器的C语言编程实现。,/*fir.c 该程序用于实现FIR滤波器*L滤波器的阶数*bi滤波器的系数,i=0,1,L-1*xi输入信号向量,i=0,1,L-1;x0对应于当前值,x1对应于上一采样值*x_in 输入信号的当前值*y_out 输出信号的当前值*/,float fir(float x_in

30、,float*x,float*b,int L)float y_out;int i;/*-*/*把上一个采样时间的输入信号向量延迟一个单元,得到当前采样时间的输入信号向量*/for(i=L-1;i0;i-)xi=xi-1;x0=x_in;/*-*/*完成FIR滤波*/y_out=0.0;for(i=0;iL;i+)y_out=y_out+bi*xi;return y_out;,直接型FIR滤波器的实现涉及到两个基本操作，一个是输入信号向量与滤波器系数向量的内积计算，另一个是输入信号向量的更新处理。在每个采样周期信号缓冲器都要更新一次，最老的采样被抛弃，而其他的信号则向缓冲器的右方移动一个单元，一

31、个新的采样被插入存储单元，并被标记。如果这个操作过程不用DSP硬件完成，那么它需要很多的时间。,7.3.3 FIR滤波器的汇编语言编程实现,处理信号缓冲器的最有效方法，是把信号采样加载到循环缓冲器中。在循环缓冲器中，采取数据保持固定、反时针方向移动地址的方式，代替保持缓冲器地址固定且正方向移动数据。信号采样的起点由指针x(n)指定，其它诸采样则沿着顺时针方向，从起点开始依次顺序加载。,图7-4 FIR滤波器的循环缓冲区,(a)信号循环缓冲区(b)系数循环缓冲区,例7-10，FIR滤波器的C55x汇编语言实现。,(1)主程序 fir_test.c/*fir_test.c*/#includemat

32、h.h#define L 64/*Number of FIR filter coefficients*/#define Fs 8000/*8000 Hz sampling frequency*/#define T 1/Fs#define f1 800/*800 Hz frequency*/#define f2 1800/*1800 Hz frequency*/#define f3 3300/*3300 Hz frequency*/#define PI 3.1415926#define w1(2*PI*f1*T)/*2*pi*f1/Fs*/#define w2(2*PI*f2*T)/*2*pi*

33、f2/Fs*/#define w3(2*PI*f3*T)/*2*pi*f3/Fs*/#define a1 0.333/*Magnitude for wave 1*/#define a2 0.333/*Magnitude for wave 2*/#define a3 0.333/*Magnitude for wave 3*/,extern int fir(int*,int*,unsigned int,int);/*Low-pass FIR filter coefficients*/int coeffL=-26,-13,14,36,31,-8,-58,-71,-15,83,139,76,-90,-

34、231,-194,50,331,383,78,-405,-654,-347,403,1024,863,-228,-1577,-1972,-453,2910,6836,9470,9470,6836,2910,-453,-1972,-1577,-228,863,1024,403,-347,-654,-405,78,383,331,50,-194,-231,-90,76,139,83,-15,-71,-58,-8,31,36,14,-13,-26;int inL;/*input buffer*/int outL;/*Output buffer*/,main()unsigned int i;float

35、 signal;unsigned int n=0;int index=0;for(i=0;iL;i+)ini=0;outi=0;while(1)signal=a1*cos(float)w1*n);signal+=a2*cos(float)w2*n);signal+=a3*cos(float)w3*n);n+;inindex=(int)(0 x7fff*signal)+0.5);outindex=fir(in,coeff,L,index);index-;if(index=-1)index=L-1;,(2)汇编语言整数fir滤波器函数：fir.asm,;fir.asm 该程序用于实现FIR滤波器,

36、可被C语言程序调用;int fir(int*,int*,unsigned int，int);参数0:AR0 输入信号缓冲区指针;参数1:AR1-FIR滤波器系数向量指针;参数2:T0-FIR 滤波器的阶数L；参数3:T1-输入信号当前值在循环缓冲区的序数;返回值:T0-输出信号当前值.def _fir _fir pshm ST1_55;现场ST1,ST2和ST3入栈 pshm ST2_55 pshm ST3_55or#0 x340,mmap(ST1_55);设置FRCT,SXMD,SATDbset SMUL;置位SMUL,mov mmap(AR0),BSA01;AR0=输入信号循环缓冲区的起始

37、地址mov mmap(AR1),BSA23;AR1=滤波器系数循环缓冲区的起始地址mov mmap(T0),BK03;设置循环缓冲区大小or#0 x5,mmap(ST2_55);AR0和AR2为循环缓冲区指针mov T1,AR0;AR0从index偏移量开始mov#0,AR2;AR2从0偏移量开始sub#2,T0;T0=L-2 mov T0,CSR;设置外部循环次数为L-1 mpym*AR0+,*AR2+,AC0;执行第一次运算|rpt CSR;启动循环macm*AR0+,*AR2+,AC0mov hi(AC0),T0;用Q15格式存放结果popm ST3_55;恢复ST1,ST2和 ST3p

38、opm ST2_55 popm ST1_55 ret.end,7.4 IIR滤波器,IIR滤波器的优点：结构简单，运算量小，可以用较少的阶数获得很高的选择性。IIR滤波器的缺点：具有相位特性差，存在稳定性问题。高阶IIR滤波器经常以串联或并联二阶环节的形式予以实现。,7.4.1 二阶IIR滤波器的结构,二阶IIR滤波器，又称为二阶基本节，分为直接型、标准型和变换型。二阶IIR滤波器传递函数为：,图7-5 二阶IIR滤波器的直接型的实现,图7-6 的信号流图,图7-7 二阶IIR滤波器的直接型的实现,由于直接型对于给定的传递函数具有最小可能的延迟数、加法器数和乘法器数，所以被称为标准型。,7.4

39、.2 高阶IIR滤波器的结构,高阶IIR滤波器的差分方程和系统函数分别为：,图7-8 高阶IIR滤波器的直接型实现（L=M+1）,图7-9 高阶IIR滤波器的串联型结构,图7-10 第k个二阶节,对于，有,其中，,注意：,7.4.3 IIR滤波器的C语言实现,例7-11,采用二维数组编写的IIR滤波器C语言程序。,temp=xin;/*xin为IIR滤波器的输入*/for(k=0;kN_IIR;k+)/*N_IIR为IIR滤波器二阶节的个数w0k=temp-a1k*w1k-a2k*w2k;/*这里,temp为本二阶节的输入,也是上一级二阶节的输出*/temp=b0k*w0k+b1k*w1k+b

40、2k*w2k;/*这里,temp为本二阶节的输出,也是下一级二阶节的输入*/w2k=w1k;w1k=w0k;xoutput=temp;/*xoutput为IIR滤波器的输出*/,aik（i=1,2）和bik（i=0,1,2）分别为第k（k=0,1,K）个二阶IIR滤波器环节的系数。,wik(i=0,1,2)是第k个二阶IIR滤波器环节中对应于图7-10中(m=0,1,2)的中间信号。,7.4.4 IIR滤波器的汇编语言实现,例7-12，采用指针编写的IIR滤波器汇编语言程序。,;iir.asm;函数原型:;void iir(int*,unsigned int,int*,int*,unsigne

41、d int,int*);入口参数：;参数0:AR0-输入信号缓冲区指针;参数1:T0-输入信号缓冲区的样本数;参数2:AR2-输出信号缓冲区指针;参数3:AR1 IIR滤波器系数数组指针;参数4:T1-二阶IIR节的个数;参数5:AR3 延迟线指针,;第k节IIR滤波器：;w(n)=x(n)-a1k*w(n-1)-a2k*w(n-2);y(n)=b0k*w(n)+b1k*w(n-1)+b2k*w(n-2);存储器分配:;暂存单元2*N_IIR 系数 5*N_IIR;AR3-w1i AR7-a1k;w1j a2k;:b2k;w2i b0k;w2j b1k;:;标定：Q14 格式;,.global

42、 _iir.sect iir_code _iir pshm ST1_55;保存ST1,ST2,ST3 pshm ST2_55 pshm ST3_55 psh T3;保存T3 pshboth XAR7;保存AR7 or#0 x340,mmap(ST1_55);设置FRCT,SXMD,SATD bset SMUL;置位SMUL sub#1,T0;样本数-1 mov T0,BRC0;设置外循环计数器 sub#1,T1,T0;二阶节个数-1 mov T0,BRC1;设置内循环计数器 mov T1,T0;设置循环缓冲区大小 sfts T0,#1 mov mmap(T0),BK03;BK03=2*二阶节个

43、数 sfts T0,#1 add T1,T0,mov mmap(T0),BK47;BK47=5*二阶节个数 mov mmap(AR3),BSA23;初始化延迟线基地址 mov mmap(AR2),BSA67;初始化系数基地址 amov#0,AR3;初始化延迟缓冲区入口 amov#0,AR7;初始化系数入口 or#0 x88,mmap(ST2_55)mov#1,T0;用于左移|rptblocal sample_loop-1;启动IIR滤波器环 mov*AR0+#14,AC0;AC0=x(n)/2(即Q14)|rptblocal filter_loop-1 masm*(AR3+T1),*AR7+,

44、AC0;AC0-=a1k*wk(n-1)masm T3=*AR3,*AR7+,AC0;AC0-=a2k*wk(n-2)mov rnd(hi(AC0T0),*AR3;wk(n-2)=wk(n)|mpym*AR7+,T3,AC0;AC0+=b2k*wk(n-2)macm*(AR3+T1),*AR7+,AC0;AC0+=b0k*wk(n-1)macm*AR3+,*AR7+,AC0;AC0+=b1k*wk(n),filter_loop mov rnd(hi(AC0#2),*AR1+;按Q15格式存放结果sample_loop popboth XAR7;恢复AR7 pop T3;恢复T3 popm ST

45、3_55;恢复ST1,ST2,ST3 popm ST2_55 popm ST1_55 ret.end*/,图7-11 IIR滤波器系数和信号缓冲区配置,7.5 快速傅里叶变换FFT,FFT算法原理库利一图基算法FFT算法的实现,7.5.1 FFT算法原理,快速傅里叶变换（FFT）是离散傅里叶变换（DFT）的一种快速算法。通过FFT算法，DFT的计算量大大减少，运算时间缩短12个数量级。DFT的变换公式为 其中 为旋转因子。,FFT之所以减少运算量，主要是利用了旋转因子的以下3点特性：对称性周期性可约性,利用这些特性可以使DFT运算中有些项进行合并，将长序列的DFT分解为短序列的DFT。DFT从

46、算法上分为按时间抽选（DIT）和按频率抽选(DIF）。基2的DIT又被称为库利一图基算法。基2的DIF又称为桑德图基算法。,7.5.2 库利一图基算法,信号流图比特反转蝶形运算,1.信号流图,2比特反转,图7-11的输入信号的顺序是按照比特反转排列的，输出序列是按照自然顺序的。比特反转就是将序列下标用二进制表示，然后将二进制数按照相反的方向排列，即得到这个序列的实际位置。按照自然排序的时域信号数据是x(0)、x(1)、x(2)、x(3)、x(4)、x(5)、x(6)、x(7)，其序号写成二进制数分别为000b、001b、010b、011b、100b、101b、110b、111b，将这些二进制数

47、前后倒转，即得到进行FFT前数据所对应的实际二进制数地址：000b、100b、010b、110b、001b、101b、011b、111b，对应的十进制数是：0、4、2、6、1、5、3、7。序号为3的存储单元，按照自然排序应该存放x(3)，但由于FFT计算规则的要求，现在应该存放x(6)。,3.蝶形运算,基2DIT FFT算法，共由M级构成，每级计算由N/2个蝶形运算构成。基本运算单元为以下所谓蝶形运算：蝶形运算中上下两个节点p、q的间距为：,7.5.3 FFT算法的实现,为了叙述简单，本书给出采用C语言编写的FFT程序，相应的汇编程序请读者自行完成。,（1）主程序：fft_test.c#inc

48、lude#include fcomplex.h/*包含浮点复数结构体定义头文件fcomplex.h*/extern void bit_rev(complex*,unsigned int);/*位反转函数声明*/extern void fft(complex*,unsigned int,complex*,unsigned int);/*fft函数声明*/,例7-13，基2DIT FFT算法的C语言实现。,extern void generator(int*,unsigned int)#define N 128/*FFT的数据个数*/#define M 7/*M=log2(N)*/#define

49、PI 3.1415926 complex XN;/*说明输入信号数组,为复数*/complex WM;/*说明旋转因子数组e(-j2PI/N),为复数*/complex temp;/*说明临时复数变量*/float xinN;float spectrumN;/*说明功率谱信号数组,为实数*/float re1N,im1N;/*说明临时变量数组,为实数*/,void main()unsigned int i,j,L,LE,LE1;/*-*/*产生旋转因子表*/for(L=1;L1;/*子FFT中的蝶形运算数目*/WL-1.re=cos(PI/LE1);WL-1.im=-sin(PI/LE1);/

50、*-*/generator(xin,N);,for(;)for(i=0;iN;i+)/*构造输入信号样本*/Xi.re=xini;Xi.im=0;/*复制到参考缓冲器*/re1i=Xi.re;im1i=Xi.im;/*启动 FFT*/bit_rev(X,M);/*以倒位次序排列X*/fft(X,M,W,1);/*执行 FFT*/*计算功率谱，验证FFT结果*/for(i=0;iN;i+)temp.re=Xi.re*Xi.re;temp.im=Xi.im*Xi.im;spectrumi=(temp.re+temp.im)*4;,（2）浮点复数基2 DIT FFT函数：fft_float.c#in