第7章应用程序设计

《第7章应用程序设计》由会员分享，可在线阅读，更多相关《第7章应用程序设计（86页珍藏版）》请在金锄头文库上搜索。

1、第第7章章 应用程序设计应用程序设计本章内容提要：本章内容提要：n定标与溢出处理定标与溢出处理n基础算术运算基础算术运算nFIR滤波器滤波器nIIR滤波器滤波器n快速傅里叶变换快速傅里叶变换(FFT)7.1定标与溢出处理定标与溢出处理n数的定标数的定标n溢出的处理方法溢出的处理方法n常用信号处理算法中的定标方法常用信号处理算法中的定标方法7.1.1 数的定标数的定标l小数定标的概念小数定标的概念u设定一个设定一个16位数的小数点处于该数中的哪一位位数的小数点处于该数中的哪一位u可以表示不同大小和不同精度的小数可以表示不同大小和不同精度的小数lQ表示法表示法u表表7-1列出了一个列出了一个16位

2、数的位数的16种种Q表示及它们所能表示及它们所能表示的十进制数值范围表示的十进制数值范围表表7-1 Q表示及数值范围表示及数值范围l不同的不同的Q所表示的数不仅范围不同，而且精度也不相同所表示的数不仅范围不同，而且精度也不相同uQ越大，数值范围越小，但精度越高越大，数值范围越小，但精度越高uQ越小，数值范围越大，但精度就越低越小，数值范围越大，但精度就越低例如，例如，Q0的数值范围是的数值范围是-32768到到+32767，其精度为，其精度为1；Q15的数值范围为的数值范围为-1到到0.9999695，精度为，精度为 1/32768 = 0.00003051对定点数而言，数值范围与精度是一对矛

3、盾。一个变量要对定点数而言，数值范围与精度是一对矛盾。一个变量要想能够表示比较大的数值范围，必须以牺牲精度为代价；想能够表示比较大的数值范围，必须以牺牲精度为代价；而想提高精度，则数的表示范围就相应地减小。在实际的而想提高精度，则数的表示范围就相应地减小。在实际的定点算法中，应该根据具体问题进行折衷处理定点算法中，应该根据具体问题进行折衷处理, , 以达到最以达到最佳效果。佳效果。l在在C55x中，中，16位整数采用补码形式表示。每个采用位整数采用补码形式表示。每个采用Qi定标的定标的16位数用位数用1个符号位、个符号位、i个小数位和个小数位和15-i个整个整数位来表示。数位来表示。表表7-2

4、 同样的数在不同定标方式下所表示的具体数值同样的数在不同定标方式下所表示的具体数值l同样一个同样一个16位数，若小数点设定的位置不同，它所表示位数，若小数点设定的位置不同，它所表示的数也就不同。的数也就不同。7.1.2 溢出的处理方法溢出的处理方法l如果算术运算结果超出寄存器所能表示的最大数就会如果算术运算结果超出寄存器所能表示的最大数就会出现溢出出现溢出u因为因为16位定点位定点DSP的动态范围有限，所以在使用的动态范围有限，所以在使用时必须注意动态范围以防溢出时必须注意动态范围以防溢出l溢出还与输入信号的特性和运算法则有关溢出还与输入信号的特性和运算法则有关 1. 溢出溢出 C55x有以下

5、几种硬件特性可以处理溢出：有以下几种硬件特性可以处理溢出：l保护位保护位 C55x的每个累加器都有的每个累加器都有8个保护位（个保护位（3932位），位），允许连续允许连续256次乘加操作而累加器不溢出次乘加操作而累加器不溢出l溢出标志位溢出标志位 C55x的每个累加器都有相关的溢出标志位，当累加的每个累加器都有相关的溢出标志位，当累加器操作结果出现溢出时，这个标志位就会置位器操作结果出现溢出时，这个标志位就会置位2. C55x的溢出处理机制的溢出处理机制l饱和方式位饱和方式位SATD和和SATAuSATD控制控制D单元的操作，单元的操作，SATA控制控制A单元的操作。单元的操作。u如果如果S

6、ATD=1，当，当D单元发生溢出时，对单元发生溢出时，对D单元的结单元的结果进行饱和处理。不管饱和方式位的值是什么，当累果进行饱和处理。不管饱和方式位的值是什么，当累加器发生溢出时，相应的溢出标志位都会被置位加器发生溢出时，相应的溢出标志位都会被置位uA单元没有溢出标志位，但如果单元没有溢出标志位，但如果SATA=1，发生溢出，发生溢出时，结果也会进行饱和处理时，结果也会进行饱和处理 饱和处理是用最近的边界值代替溢出结果。饱和处理是用最近的边界值代替溢出结果。例如，例如，1616位寄存器的范围是位寄存器的范围是8000h8000h（最小负数）（最小负数）7FFFh7FFFh（最大正数），饱和处

7、理就是用（最大正数），饱和处理就是用7FFFh7FFFh代替比代替比7FFFh7FFFh大的结果；用大的结果；用8000h8000h代替比代替比8000h8000h小的结果。小的结果。l饱和。饱和是一种处理溢出的方法，但是饱和会剪掉部饱和。饱和是一种处理溢出的方法，但是饱和会剪掉部分输出信号，可能会引起信号失真和引起系统非线性。分输出信号，可能会引起信号失真和引起系统非线性。l输入定标。分析所要使用的系统，假定最坏的情况，然输入定标。分析所要使用的系统，假定最坏的情况，然后对输入信号定标，以防止溢出。但是这种方法会极大后对输入信号定标，以防止溢出。但是这种方法会极大地降低输出信号的精确度。地降

8、低输出信号的精确度。l固定定标。假定最坏的情况，对中间结果定标。这种方固定定标。假定最坏的情况，对中间结果定标。这种方法可以防止溢出，同时增加了系统的信噪比。法可以防止溢出，同时增加了系统的信噪比。l动态定标。可以监测中间结果的范围，只在需要的时候动态定标。可以监测中间结果的范围，只在需要的时候对中间结果定标。这种方法可以防止溢出但会增加计算对中间结果定标。这种方法可以防止溢出但会增加计算量。量。3. 溢出的处理方法溢出的处理方法7.1.3 常用信号处理常用信号处理 算法中的定标方法算法中的定标方法lFIR滤波器的定标方法滤波器的定标方法u在在FIR滤波器中处理溢出的最好方法是设计时使滤滤波器

9、中处理溢出的最好方法是设计时使滤波器的增益小于波器的增益小于1，这样就不需要对输入信号定标。，这样就不需要对输入信号定标。这种方法和累加器的保护位结合起来，可以有效地这种方法和累加器的保护位结合起来，可以有效地防止溢出。防止溢出。u由于对信号处理的负面影响，在由于对信号处理的负面影响，在FIR滤波器中不使滤波器中不使用固定定标和输入定标。如果不考虑计算量，在用固定定标和输入定标。如果不考虑计算量，在FIR滤波器中可以使用动态定标。对一些类型的音滤波器中可以使用动态定标。对一些类型的音频信号，饱和处理也是一种常用的方法。频信号，饱和处理也是一种常用的方法。lIIR滤波器的定标方法滤波器的定标方法

10、uIIR滤波器的定点实现推荐使用多个二阶基本节级滤波器的定点实现推荐使用多个二阶基本节级联组成，这样可以减小高阶滤波器频率响应灵敏度。联组成，这样可以减小高阶滤波器频率响应灵敏度。由于滤波器系数的量化引入误差，因此避免溢出对由于滤波器系数的量化引入误差，因此避免溢出对IIR滤波器非常重要。滤波器非常重要。u可以通过把中间结果保存在处理器累加器来避免节可以通过把中间结果保存在处理器累加器来避免节间数据溢出。为防止在第间数据溢出。为防止在第k阶内部发生数据溢出，阶内部发生数据溢出，需要用增益系数对滤波器的单位脉冲响应（前馈通需要用增益系数对滤波器的单位脉冲响应（前馈通道）定标。道）定标。u动态标定

11、方法。在每个阶段滤波器内部状态都被减动态标定方法。在每个阶段滤波器内部状态都被减半，以提高指令周期换取为代价提高了结果的精度。半，以提高指令周期换取为代价提高了结果的精度。l FFT的定标方法的定标方法u在FFT操作里，每次蝶形运算后数据平均增加一位。输入定标需要移位 （FFT长度为N），这会导致在计算FFT之前就衰减 6 dB。u在固定定标中，每级蝶形运算输出除以2，这是最常用的FFT定标方法，因为它简单而且有比较好的信噪比。但是，对于大的FFT，这种定标可能会使信息丢失。u另一种方法是动态定标，即在输出溢出时再除以2。在这种情况下，会在这个过程中指定一个变量，每定标一次变量的值加1，计算结

12、束后根据变量的值把结果乘以一个系数。动态定标的信噪比最好，但会增加FFT循环次数。7.2 基础算术运算基础算术运算n加减运算加减运算n乘法运算乘法运算n除法运算除法运算n小数乘法小数乘法7.2.1 加减运算加减运算l在数字信号处理中，加减运算是常见的算术运算。在数字信号处理中，加减运算是常见的算术运算。一般使用一般使用16位或位或32位加减运算，数值分析、浮点位加减运算，数值分析、浮点运算和其它操作可能需要运算和其它操作可能需要32位以上的运算。位以上的运算。lC55x有直接完成有直接完成16位或位或32位加减运算的指令，位加减运算的指令，但没有能直接完成多字加减运算的指令。要进行但没有能直接

13、完成多字加减运算的指令。要进行多字加减运算，需要通过编程方法实现。多字加减运算，需要通过编程方法实现。l 以下指令可在单周期内完成以下指令可在单周期内完成32位加法运算：位加法运算：MOV40 dbl(Lmem),ACxADD dbl(Lmem),ACxl 64位的高位的高32位加法要考虑低位加法要考虑低32位加法产生的进位，位加法产生的进位，使用以下指令：使用以下指令：ADD uns(Smem), CARRY, ACxl 以下指令可在单周期内完成以下指令可在单周期内完成32位减法运算：位减法运算：MOV40 dbl(Lmem),ACxSUB dbl(Lmem),ACxl 64位的高位的高32

14、位减法要考虑低位减法要考虑低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 (偶地址); + Y3 Y2 Y1 Y0 X2; - X1; W3 W2 W1 W0 X0; AR2 - Y3 (偶地址); Y2; Y1; Y0; AR3 - W3 (偶地址); W2; W1; W0;*.sect

15、.text.align 4.globalstart.symstart,start, 36, 2, 0 start: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), AC0 ; AC0 = X3 X2 ADD uns(*AR2(#1),CARRY,AC0 ; AC0 = X3 X2 + 00

16、Y2 + CARRY ADD *AR2 X3 (偶地址); 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,0 start:MOV #0100h,AR1 MOV #0104h,AR2 MOV #0108h,AR3L1: MOV40 dbl(*AR1(#2),AC0 ;AC0=X1X0 SUB dbl(*AR2(#2),AC0 ;AC0=X1X0-Y1Y0

17、MOV AC0,dbl(*AR3(#2) ;保存W1W0. MOV40 dbl(*AR1),AC0 ;AC0=X3X2 SUB uns(*AR2(#1),BORROW,AC0;AC0=X3X2-00Y2-BORROW SUB *AR2 X1; Y1 Y0 32位操作数 X0; W3 W2 W1 W0 64位结果 AR1 - Y1; Y0; 入口条件: AR2 - W0; SXMD = 1 (允许符号扩展) W1; SATD = 0 (不做饱和处理) W2; FRCT = 0 (关小数模式) W3; ;限制条件：延迟链和输入序列必须指定为长字类型。;*.sect.text.align 4.glo

18、bal 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(*AR0-),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*Y0

19、16+X1*Y0+X0*Y1MOV AC0,*AR2+ ;保存W1MACM *AR0,*AR1,AC0#16,AC0 ;AC0=AC016+X1*Y1MOV AC0,*AR2+ ;保存W2MOV HI(AC0),*AR2 ;保存W3 B L17.2.3 除法运算除法运算lC55x没有提供硬件除法器，也没有提供专门的除法指令，没有提供硬件除法器，也没有提供专门的除法指令，要实现除法运算需借助于条件减法指令要实现除法运算需借助于条件减法指令SUBC和重复指令和重复指令RPT。l根据被除数绝对值与除数绝对值的大小关系，除法的实根据被除数绝对值与除数绝对值的大小关系，除法的实现过程略有不同：现过程略有

20、不同：u当当|被除数被除数|被除数; AR1-除数; AR2-商; AR3-余数; 注：; 无符号除法，被除数、除数均为16位; 关闭符号扩展，被除数、除数均为正数; 运算完成后AC0(15-0)为商，AC0(31-16)为余数;* .sect.text.align 4.global start.symstart,start,36,2,0 start:MOV #0100h,AR0 MOV #0101h,AR1 MOV #0102h,AR2 MOV #0103h,AR3L1:BCLR SXMD ; 清零SXMD (关闭符号扩展)MOV *AR0,AC0 ; 把被除数放入AC0RPT #(16-1

21、) ; 执行subc 16次SUBC *AR1,AC0,AC0 ; AR1指向除数MOV AC0,*AR2 ; 保存商MOV HI(AC0),*AR3 ; 保存余数 B L1例例7-5，无符号，无符号32位除位除16位整数除法。位整数除法。文件名为：文件名为：udiv32o16.asm。.mmregs.model call=c55_std.model mem=large;*; 指针分配; AR0-被除数高位;被除数低位; AR1-除数; AR2-商高位; 商低位; AR3-余数;注：;无符号除法，被除数为32位，除数为16位;关闭符号扩展，被除数、除数均为正数;第一次除法之前，把被除数高位存入

22、AC0;第一次除法之后，把商的高位存入AC0(15-0);第二次除法之前，把被除数低位存入AC0; 第二次除法之后，AC0(15-0)为商的低位，AC0(31-16)为余数;*.sect.text.align 4.globalstart.symstart,start,36,2,0start: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,A

23、C0,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 AC0,*AR2+ ; 保存商的低位MOV HI(AC0),*AR3 ; 保存余数BSET SXMD ; 置位SXMD (打开符号扩展) B L1例例7-6，带符号，带符号16位除位除16位整数除法。位整数除法。文件名为：文件名为：sdiv16o16.

24、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.globalstart.symstart,start,36,2,0 start:MOV #0100h,AR0 MOV #0101h,AR1 MOV #0102h,AR2 MOV #0103h,AR3

25、L1: 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 - |除数|MOV HI(AC1),*AR3 ; 保存余数MOV AC1,*AR2 ; 保存商SFTS AC1,#16 ; 对商移位: 把符号位

26、放在最高位NEG AC1,AC1 ; 对商求反XCCPART label,AC0被除数高半部分(NumH)(偶地址);被除数高半部分(NumL); AR1-除数(Den);AR2-商的高半部分(QuotH)(偶地址);商的低半部分(QuotL);AR3-余数(Rem); 注：; 带符号除法，被除数为32位，除数为16位;打开符号扩展，被除数、除数可为负数;除法运算之前，期望的商的符号存入AC0;第一次除法运算之前，把被除数的高半部分存入AC1;第一次除法运算之后，把商的高半部分存入AC1(15-0);第二次除法运算之前，把被除数的低半部分存入AC1;第二次除法运算之后，把商的低半部分存入AC1

27、(15-0),余数存入AC1(31-16);*.sect.text.align 4.global start.symstart,start,36,2,0start:MOV #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 =

28、 abs(Den) MOV40 dbl(*AR0),AC1 ; AC1 = NumH NumLABS AC1,AC1 ; AC1 = abs(Num)MOV AC1,dbl(*AR2) ; QuotH = abs(NumH) ; QuotL = abs(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 ;

29、 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),AC1 ; 否则,对商取反.NEG AC1,AC1MOV AC1,dbl(*AR2)skip: B L17.2.4 小数乘法小数乘法l在定点在定点DSP的某些应用中，整数运算很难满足要求。这的某些应用中，整数运算很难满足要求。这是因为它自身存在缺陷：是因为它自身存在缺陷：u两个两个1

30、6位整数相乘，乘积总是位整数相乘，乘积总是“向左增长向左增长”（即小数（即小数点左侧的位数增加），这意味着多次相乘后，乘积将点左侧的位数增加），这意味着多次相乘后，乘积将很快超出定点器件的数据范围。很快超出定点器件的数据范围。u保存保存32位乘积到存储器，要占用位乘积到存储器，要占用2个个CPU周期和周期和2个个字的存储器空间。字的存储器空间。u由于乘法器都是由于乘法器都是16位相乘，因此将位相乘，因此将32位乘积再作为位乘积再作为乘法器的输入时就显得较繁琐，不能胜任递归运算。乘法器的输入时就显得较繁琐，不能胜任递归运算。l为了克服这些缺陷，在实际应用中更多采用的是小数为了克服这些缺陷，在实际

31、应用中更多采用的是小数运算。小数运算具有如下优点：运算。小数运算具有如下优点： u乘积总是乘积总是“向右增长向右增长”。这就意味着超出定点器件。这就意味着超出定点器件数据范围的将是不太感兴趣的部分。数据范围的将是不太感兴趣的部分。u既可以存储既可以存储32位乘积，也可以近存储高位乘积，也可以近存储高16位乘积，位乘积，这就允许用较少的资源保存结果。这就允许用较少的资源保存结果。u可以用于递归运算。可以用于递归运算。 例例7-8，两个，两个Q31格式有符号小数相乘，得到一格式有符号小数相乘，得到一个个Q31格式结果。文件名为：格式结果。文件名为：mpyQ31.asm。 需要注意的是：需要注意的是

32、：u两个带符号小数相乘，所得乘积带有两个带符号小数相乘，所得乘积带有2个符号位。个符号位。u为了解决冗余符号位的问题，需要在程序中设定状为了解决冗余符号位的问题，需要在程序中设定状态寄存器态寄存器ST1中的中的FRCT（小数方式）为（小数方式）为1，这样当，这样当乘法器将结果传送至累加器时就会自动左移乘法器将结果传送至累加器时就会自动左移1位。位。.mmregs.model call=c55_std.model mem=large;*;操作数取自数据存储器，运算结果送回数据存储器。; 数据存储: 指针分配:; X1 X0 Q31操作数 AR0 - X1; Y1 Y0 Q31操作数 X0; W1

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

34、*AR1+,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 L17.3 FIR滤波器滤波器l数字滤波器是数字滤波器是DSP的基本应用，有的基本应用，有2种基本类型：种基本类型：u有限冲激响应滤波器有限冲激响应滤波器FIRu无限冲激滤波器无限冲激滤波器IIRl一般来说，如果需要线性相位则选择用一般来说，如果需要线性相位则选择用FIR滤波器，滤波器，对于相位要求不敏感的场合可以选用对于相

35、位要求不敏感的场合可以选用IIR滤波器。滤波器。l本节主要讨论本节主要讨论FIR滤波器的滤波器的DSP实现方法，有关实现方法，有关IIR滤波器的实现将在下一节中介绍。滤波器的实现将在下一节中介绍。7.3.1 FIR滤波器的基本结构滤波器的基本结构l一个一个FIR滤波器的输出序列和输入序列之间的关系，滤波器的输出序列和输入序列之间的关系，满足差分方程：满足差分方程：l传递函数为传递函数为 l FIR滤波器的结构：l FIR滤波器的单位冲激响应是一个有限长序列。若为实数，且满足偶对称或奇对称的条件，则FIR滤波器具有线性相位特性。 偶对称线性相位FIR滤波器的差分方程为：7.3.2 FIR滤波器的

36、滤波器的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, float *x,float *b,int L)float y_out;int i;/* - */* 把上一个采样时间的输入信号向量延迟一个单元把上一个采样时间的输入信号向量延迟一个单元,得到当前

37、得到当前采样时间的输入信号向量采样时间的输入信号向量 */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;u 直接型直接型FIRFIR滤波器的实现涉及到两个基本操作，一个是滤波器的实现涉及到两个基本操作，一个是输入信号向量与滤波器系数向量的内积计算，另一个是输入信号向量与滤波器系数向量的内积计算，另一个是输入信号向量的更新处理。输入信号向量的更新处理。u 在每个采样周期信号缓冲器都要更新一次，最老在每个采样

38、周期信号缓冲器都要更新一次，最老的采样被抛弃，而其他的信号则向缓冲器的右方移动的采样被抛弃，而其他的信号则向缓冲器的右方移动一个单元，一个新的采样被插入存储单元，并被标记。一个单元，一个新的采样被插入存储单元，并被标记。 如果这个操作过程不用如果这个操作过程不用DSPDSP硬件完成，那么它需要硬件完成，那么它需要很多的时间。很多的时间。7.3.3 FIR滤波器的汇编语言编程实现滤波器的汇编语言编程实现l处理信号缓冲器的最有效方法，是把信号采样加载到循环处理信号缓冲器的最有效方法，是把信号采样加载到循环缓冲器中。缓冲器中。l在循环缓冲器中，采取数据保持固定、反时针方向移动地在循环缓冲器中，采取数

39、据保持固定、反时针方向移动地址的方式，代替保持缓冲器地址固定且正方向移动数据。址的方式，代替保持缓冲器地址固定且正方向移动数据。l信号采样的起点由指针信号采样的起点由指针x(n) 指定，其它诸采样则沿着顺指定，其它诸采样则沿着顺时针方向，从起点开始依次顺序加载。时针方向，从起点开始依次顺序加载。图图7-4 FIR滤波器的循环缓冲区滤波器的循环缓冲区(a)信号循环缓冲区信号循环缓冲区 (b)系数循环缓冲区系数循环缓冲区例例7-10，FIR滤波器的滤波器的C55x汇编语言实现。汇编语言实现。(1) 主程序 fir_test.c/* fir_test.c */#includemath.h#defin

40、e 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*

41、f2*T) /* 2*pi*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,

42、31,-8,-58,-71,-15,83,139,76,-90,-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; /*

43、Output buffer */main() unsigned int i; float 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)(0x7fff*signal)+0.5); outindex = fir(in,coeff,L,index); index-; if(in

44、dex=-1) index=L-1; (2) 汇编语言整数汇编语言整数fir滤波器函数：滤波器函数：fir.asm; fir.asm 该程序用于实现FIR滤波器,可被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_5

45、5 pshm ST3_55or #0x340,mmap(ST1_55) ;设置FRCT,SXMD,SATDbset SMUL ;置位SMULmov mmap(AR0),BSA01 ;AR0=输入信号循环缓冲区的起始地址mov mmap(AR1),BSA23 ;AR1=滤波器系数循环缓冲区的起始地址mov mmap(T0),BK03 ;设置循环缓冲区大小or #0x5,mmap(ST2_55) ;AR0和AR2为循环缓冲区指针mov T1,AR0 ; AR0从index偏移量开始mov #0,AR2 ; AR2从0偏移量开始sub #2,T0 ;T0=L-2 mov T0,CSR ;设置外部循环

46、次数为L-1 mpym *AR0+,*AR2+,AC0 ;执行第一次运算| rpt CSR ;启动循环macm *AR0+,*AR2+,AC0mov hi(AC0),T0 ; 用Q15格式存放结果popm ST3_55 ; 恢复ST1, ST2和 ST3popm ST2_55 popm ST1_55 ret.end7.4 IIR滤波器滤波器lIIR滤波器的优点：结构简单，运算量小，可以用滤波器的优点：结构简单，运算量小，可以用较少的阶数获得很高的选择性。较少的阶数获得很高的选择性。lIIR滤波器的缺点：具有相位特性差，存在稳定性滤波器的缺点：具有相位特性差，存在稳定性问题。问题。l高阶高阶II

47、R滤波器经常以串联或并联二阶环节的形式滤波器经常以串联或并联二阶环节的形式予以实现。予以实现。7.4.1 二阶二阶IIR滤波器的结构滤波器的结构l二阶二阶IIR滤波器，又称为二阶基本节，分为直接型、滤波器，又称为二阶基本节，分为直接型、标准型和变换型。标准型和变换型。l二阶二阶IIR滤波器传递函数为：滤波器传递函数为：图图7-5 二阶二阶IIR滤波器的直接型滤波器的直接型的实现的实现图图7-6 的信号流图的信号流图图图7-7 二阶二阶IIR滤波器的直接型滤波器的直接型的实现的实现由于直接型由于直接型对于给定的传递函数具有最小可能的延迟对于给定的传递函数具有最小可能的延迟数、加法器数和乘法器数，

48、所以被称为标准型。数、加法器数和乘法器数，所以被称为标准型。7.4.2 高阶高阶IIR滤波器的结构滤波器的结构l高阶高阶IIR滤波器的差分方程和系统函数分别为：滤波器的差分方程和系统函数分别为：图图7-8 高阶高阶IIR滤波器的直接型滤波器的直接型实现（实现（L=M+1）图图7-9 7-9 高阶高阶IIRIIR滤波器的串联型结构滤波器的串联型结构图图7-10 第第k个二阶节个二阶节对于 ，有其中，注意：7.4.3 IIR滤波器的滤波器的C语言实现语言实现例例7-11,采用二维数组编写的采用二维数组编写的IIR滤波器滤波器C语言程序。语言程序。temp=xin; /*xin为IIR滤波器的输入

49、*/for(k=0;k w1i AR7 - a1k; w1j a2k; : b2k; w2i b0k; w2j b1k; : : ; 标定：Q14 格式; .global _iir .sect iir_code _iir pshm ST1_55 ; 保存ST1,ST2,ST3 pshm ST2_55 pshm ST3_55 psh T3 ; 保存T3 pshboth XAR7 ; 保存AR7 or #0x340,mmap(ST1_55) ;设置FRCT,SXMD,SATD bset SMUL ;置位SMUL sub #1,T0 ; 样本数 - 1 mov T0,BRC0 ; 设置外循环计数器

50、sub #1,T1,T0; 二阶节个数 - 1 mov T0,BRC1; 设置内循环计数器 mov T1,T0 ;设置循环缓冲区大小 sfts T0,#1 mov mmap(T0),BK03;BK03=2*二阶节个数 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 #0x88,mmap(ST2_55) mov #1,T0 ;

51、用于左移| rptblocal sample_loop-1 ;启动IIR滤波器环 mov *AR0+ #14,AC0 ;AC0=x(n)/2(即Q14)| rptblocal filter_loop-1 masm *(AR3+T1),*AR7+,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

52、) 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 ST3_55 ;恢复ST1,ST2,ST3 popm ST2_55 popm ST1_55 ret .end */图图7-11 IIR滤波器系数和信号缓冲区配置滤波器系数和信号缓冲区配置7.5 快速傅里叶变换快速傅里叶变换FFTnFFT算法原理算法原理n库利一图基算法库利一图基算法nFFT算法的实现算法的实现7.5.1 FFT

53、算法原理算法原理l快速傅里叶变换（快速傅里叶变换（FFT）是离散傅里叶变换（）是离散傅里叶变换（DFT）的一种快速算法。通过的一种快速算法。通过FFT算法，算法，DFT的计算量大大的计算量大大减少，运算时间缩短减少，运算时间缩短12个数量级。个数量级。lDFT的变换公式为的变换公式为 其中其中 为旋转因子。为旋转因子。l FFT之所以减少运算量，主要是利用了旋转因子的以之所以减少运算量，主要是利用了旋转因子的以下下3点特性：点特性：u对称性对称性u周期性周期性u可约性可约性利用这些特性可以使利用这些特性可以使DFT运算中有些项进行合并，将长序列的运算中有些项进行合并，将长序列的DFT分解为短序

54、列的分解为短序列的DFT。DFT从算法上分为按时间抽选（从算法上分为按时间抽选（DIT）和按频率抽选）和按频率抽选(DIF）。）。基基2的的DIT又被称为库利一图基算法。基又被称为库利一图基算法。基2的的DIF又称为桑德又称为桑德图基算法。图基算法。 7.5.2 库利一图基算法库利一图基算法n信号流图信号流图n比特反转比特反转n蝶形运算蝶形运算1. 信号流图信号流图2比特反转比特反转l图图7-11的输入信号的顺序是按照比特反转排列的，输出的输入信号的顺序是按照比特反转排列的，输出序列是按照自然顺序的。比特反转就是将序列下标用二序列是按照自然顺序的。比特反转就是将序列下标用二进制表示，然后将二进

55、制数按照相反的方向排列，即得进制表示，然后将二进制数按照相反的方向排列，即得到这个序列的实际位置。到这个序列的实际位置。l按照自然排序的时域信号数据是按照自然排序的时域信号数据是x(0)、x(1)、x(2)、x(3)、x(4)、x(5)、x(6)、x(7)，其序号写成二进制数，其序号写成二进制数分别为分别为000b、001b、010b、011b、100b、101b、110b、111b，将这些二进制数前后倒转，即得到进行，将这些二进制数前后倒转，即得到进行FFT前数据所对应的实际二进制数地址：前数据所对应的实际二进制数地址：000b、100b、010b、110b、001b、101b、011b、1

56、11b，对应的，对应的十进制数是：十进制数是：0、4、2、6、1、5、3、7。序号为。序号为3的存的存储单元，按照自然排序应该存放储单元，按照自然排序应该存放x(3)，但由于，但由于FFT计算规计算规则的要求，现在应该存放则的要求，现在应该存放x(6)。 3.蝶形运算蝶形运算l 基基2DIT FFT算法，共由算法，共由M级构成，每级计算由级构成，每级计算由N/2个蝶形运算构成。个蝶形运算构成。l 基本运算单元为以下所谓蝶形运算：基本运算单元为以下所谓蝶形运算：l 蝶形运算中上下两个节点蝶形运算中上下两个节点p、q的间距为：的间距为： 7.5.3 FFT算法的实现算法的实现l为了叙述简单，本书给

57、出采用为了叙述简单，本书给出采用C语言编写的语言编写的FFT程序，程序，相应的汇编程序请读者自行完成。相应的汇编程序请读者自行完成。（1）主程序：fft_test.c#include #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语言实

58、现。语言实现。extern void generator(int *,unsigned int)#define N 128 /* FFT的数据个数 */#define M 7 /* M=log2(N) */#define PI 3.1415926 complex XN; /* 说明输入信号数组,为复数 */complex WM; /* 说明旋转因子数组e(-j2PI/N),为复数 */ complex temp; /* 说明临时复数变量 */float xinN;float spectrumN; /* 说明功率谱信号数组,为实数 */float re1N,im1N; /* 说明临时变量数组,为

59、实数 */void main() unsigned int i,j,L,LE,LE1; /*- */* 产生旋转因子表 */ for (L=1; L=M; L+) LE=11; /* 子FFT中的蝶形运算数目*/ WL-1.re = cos(PI/LE1); WL-1.im = -sin(PI/LE1); /* - */ generator(xin,N); for (;) for (i=0; iN; i+) /* 构造输入信号样本构造输入信号样本 */ Xi.re =xini; Xi.im = 0; /* 复制到参考缓冲器复制到参考缓冲器 */ re1i = Xi.re; im1i = Xi.

60、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 #include fcomplex.hvoid fft(complex *X,unsigne

61、d int M,complex *W,unsigned int SCALE) complex temp; /* 复变量临时存储器复变量临时存储器 */ complex U; /* 旋转因子旋转因子Wk */ unsigned int i,j; unsigned int id; /* 蝶形运算中下位节点的序号蝶形运算中下位节点的序号 */ unsigned int N=1M; /* FFT 的点数的点数*/ unsigned int L; /* FFT 的级序号的级序号 */ unsigned int LE; /* L级子级子FFT的点数的点数 */ unsigned int LE1; /* L

62、级子级子FFT蝶形运算的个数蝶形运算的个数 */ float scale; scale = 0.5; for(L=1;L=M;L+) LE=11; U.re = 1.0; U.im = 0.; for(j=0;jLE1;j+) for(i=j;iN;i+=LE) /*进行蝶形计算进行蝶形计算*/ id=i+LE1; temp.re=(Xid.re*U.re-Xid.im*U.im)*scale; temp.im=(Xid.im*U.re+Xid.re*U.im)*scale; Xid.re=Xi.re*scale-temp.re; Xid.im=Xi.im*scale-temp.im; Xi.

63、re=Xi.re*scale+temp.re; Xi.im=Xi.im*scale+temp.im; /*递推计算递推计算Wk*/ temp.re=U.re*WL-1.re-U.im*WL-1.im; U.im=U.re*WL-1.im+U.im*WL-1.re; U.re=temp.re; （3）位反转函数：bit_rev.c #include fcomplex.hvoid bit_rev(complex *X,unsigned int M) complex temp; unsigned int i,j,k; unsigned int N=11; for (j=0,i=1; iN-1; i+

64、) k=N2; while(k=1; j+=k; if(ij) temp=Xj; Xj=Xi;Xi=temp;（4）信号发生器函数：generator.c/* generator.c -该程序用于产生一组信号样本 */#include math.h#define PI=3.14159265358972 #define Fs=8000 ;采样频率设为8000Hz#define T=1/Fs ;采样时间为0.25ms#define f1=500 ;信号源频率1取为500Hz#define a1=0.5 ;信号源幅度1取为0.5#define w1=2*PI*f1*T void generator(int *x,unsigned int N) unsigned int i; for(i=0;iN;i+) Xi=a1*cos(float)w1*i);（5）复数结构定义头文件：）复数结构定义头文件：fcomplex.hstruct cmpxfloat re;float im;typedef struct cmpx complex;