DSP原理与应用 第7章 应用程序设计 刘忠国：liuzhg@sdu.edu.cn 电话:18764171197 ; 微信: jnliuzg 第7章 应用程序设计 刘忠国：liuzhg@sdu.edu.cn 电话:18764171197 ; 微信: jnliuzg 山东大学生物医学工程
网站: http://course.sdu.edu.cn/G2S/dsp.cc http://www.ti.com 任课教师:刘忠国 网站: http://course.sdu.edu.cn/G2S/dsp.cc http://www.ti.com TMS320C55x DSP Programmer’s Guide Preliminary Draft (SPRU376A, 2001年) TMS320C55x DSP Library Programmer’s Reference (SPRU422J, 2009年) TMS320C55x Assembly Language Tools User’s Guide (SPRU280H, 2002年) TMS320C55x DSP Mnemonic Instruction Set Reference Guide (SPRU374G, 2002年)
第7章 应用程序设计 内容提要： 7.1 定标与溢出处理 7.2 基础算术运算 7.3 FIR滤波器 7.4 IIR滤波器 第7章 应用程序设计 内容提要： 7.1 定标与溢出处理 7.2 基础算术运算 7.3 FIR滤波器 7.4 IIR滤波器 7.5 快速傅里叶变换(FFT) 7.6 DSPLIB的使用
7.1 定标与溢出处理 7.1.1 数的定标 7.1.2 溢出的处理方法 7.1.3 常用信号处理算法中的定标方法
7.1.1 数的定标 1.小数定点数的定标(scale) 计算机有两种数据表示方法： 定点数表示法 浮点数表示法 定点数就是小数点位置固定的数。数据的表示形式是： 符号位 整数部分.小数部分 C55x DSP是定点芯片，采用补码形式来表示无符号数和有符号数。
7.1.1 数的定标 2. Q表示法 CPU在执行指令时，并不知道处理的数据是整数还是小数，也不能指出小数点的位置。在编程时必须由程序员指出一个数的小数点处于哪一位，这就是定标(scale)的概念。 常用的定标方法是Q表示法, Q表示法用Q0, Q1,… Qi, …, Q15 分别表示小数点在D0位之后、D1位之后、…Di位之后、…D15位之后。 严格地, 是Qm.n 格式: n位小数, m位整数, 1位符号位, 共m+n+1位二进制数。 也称S表示法Sm.n
表7-1 Q表示及数值范围 Q表示 小数点位置 整数位 小数位 精度--幂 Q15 在D15之后 15 -1~0.9999695 2 -15 十进制表示范围 精度--幂 Q15 在D15之后 15 -1~0.9999695 2 -15 Q14 在D14之后 1 14 -2~1.9999390 2 -14 Q13 在D13之后 2 13 -4~3.9998779 2 -13 Q12 在D12之后 3 12 -8~7.9997559 2 -12 Q11 在D11之后 4 11 -16~15.9995117 2 -11 Q10 在D10之后 5 10 -32~31.9990234 2 -10 Q9 在D9之后 6 9 -64~63.99804375 2 -9 Q8 在D8之后 7 8 -128~127.9960938 2 -8 Q7 在D7之后 -256~255.9921875 2 -7 Q6 在D6之后 -512~511.9804375 2 -6 Q5 在D5之后 -1024~1023.96875 2 -5 Q4 在D4之后 -2048~2047.9375 2 -4 Q3 在D3之后 -4096~4095.875 2 -3 Q2 在D2之后 -8192~8191.75 2 -2 Q1 在D1之后 14 1 -16384~16383.5 2 -1 Q0 在D0之后 15 -32768~32767 2 0
不同的Q所表示的数不仅范围不同, 而且精度也不相同: 1/32768 = 0.00003051 对定点数而言, 数值范围与精度是一对矛盾。一个变 量要想能够表示比较大的数值范围, 必须以牺牲精度 为代价; 而想提高精度, 则数的表示范围就相应地减 小。 在实际的定点算法中，应该根据具体问题进行折衷处 理, 以达到最佳效果。
在C55x中, 16位整数采用补码形式表示。每个采用Qi定标的16位数用1个符号位、i个小数位和15-i个整数位来表示。 . . .
表7-2 同样的数在不同定标方式下所表示的具体数值 同样一个16位数，若小数点设定的位置不同，它所表示的数也就不同。 定标不同 表7-2 同样的数在不同定标方式下所表示的具体数值 Q格式 十六进制数 二进制数 十进制数 Q0 2000H 0010 0000 0000 0000b 8192 Q14 0.5 Q15 0.25 E000H 1110 0000 0000 0000b -8192 -0.5 -0.25 . . . . . .
7.1 定标与溢出处理 7.1.1 数的定标 7.1.2 溢出的处理方法 7.1.3 常用信号处理算法中的定标方法 7.1 定标与溢出处理 7.1.1 数的定标 7.1.2 溢出的处理方法 7.1.3 常用信号处理算法中的定标方法 TMS320C55x DSP Programmer’s Guide Preliminary Draft (SPRU376A, 2001年) 5.5 Methods of Handling Overflows
7.1.2 溢出的处理方法 1. 溢出 如果算术运算结果超出寄存器所能表示的最大 数就会出现溢出; 因为16位定点DSP的动态范围有限，所以在使 用时必须注意动态范围以防溢出; 溢出还与输入信号的特性和运算法则有关 。
2. C55x的溢出处理机制 C55x有以下几种硬件特性可以处理溢出： 保护位 溢出标志位 C55x的每个累加器都有相关的溢出标志位, 当累加器操作结果出现溢出时,这个标志位就会置位。 M40=1
SATD控制D单元的操作, SATA控制A单元的操作。 饱和方式位SATD和SATA SATD控制D单元的操作, SATA控制A单元的操作。 若SATD=1, 当D单元发生溢出时, 对D单元结果进行饱和处理。不管饱和方式位值是什么, 当累加器发生溢出时, 相应溢出标志位(ACOVx)都会被置位; A单元没有溢出标志位, 但如果SATA=1, 发生溢出时, 结果也会进行饱和处理。 ☼ 饱和处理是用最近的边界值代替溢出结果。 例, 16位寄存器范围是8000h(最小负数)～7FFFh(最大正 数), 饱和处理就是: 用7FFFh代替比7FFFh大的结果; 用8000h代替比8000h小的结 果。
3. 溢出的处理方法 饱和。饱和是一种处理溢出的方法，但是饱和会剪掉 部分输出信号，可能会引起信号失真和引起系统非线 性。 输入定标scaling。分析所要使用的系统，假定最坏 的情况，然后对输入信号定标，以防止溢出。但是这 种方法会极大地降低输出信号的精确度。 固定定标scaling。假定最坏情况,对中间结果定标。 这种方法可防止溢出, 同时增加了系统的信噪比。 动态定标scaling。可以监测中间结果的范围，只在 需要的时候对中间结果定标。这种方法可以防止溢出 但会增加计算量(计算时间)。
7.1 定标与溢出处理 7.1.1 数的定标 7.1.2 溢出的处理方法 7.1.3 常用信号处理算法中的定标方法 7.1 定标与溢出处理 7.1.1 数的定标 7.1.2 溢出的处理方法 7.1.3 常用信号处理算法中的定标方法 TMS320C55x DSP Programmer’s Guide Preliminary Draft (SPRU376A, 2001年) 5.5 Methods of Handling Overflows
7.1.3 常用信号处理算法中的定标方法 FIR滤波器的定标scaling方法 在FIR滤波器中一般不用固定定标和输入定标,因其影响信号分辨率(基本是乘累加1次对应1位)。在FIR滤波器中可使用动态定标, 前提是由此导致的计算时间增加可以承受。饱和处理也是某些音频信号一种常用的方法。
IIR滤波器的定标scaling方法 由于滤波器系数的量化引入舍入误差,用定点法实现IIR滤波器时, 推荐使用多个二阶基本节级联组成, 这样可以减小高阶滤波器对误差的频率响应灵敏度。另外避免溢出对IIR滤波器也非常重要。 可以通过把中间结果保存在处理器累加器来避免 节间数据溢出。为防止在第k阶内部发生数据溢出， 需要用增益系数对滤波器的单位脉冲响应(前馈通 道)定标。 动态定标scaling方法。在每个阶段若检测到溢出, 滤波器内部状态,定标都被减半, 以提高指令周期换取为代价提高了结果的精度。
在FFT操作里，每次蝶形运算后数据平均增加一位。 输入定标需要移位 (FFT长度为N), 这会导致在 计算FFT之前就衰减 6 dB。 (右移1位即÷2, ) 在FFT操作里，每次蝶形运算后数据平均增加一位。 输入定标需要移位 (FFT长度为N), 这会导致在 计算FFT之前就衰减 6 dB。 在固定定标中, 每级蝶形运算输出除以2, 这是最常用 的FFT定标方法, 因为它简单而且有比较好的信噪比。 但是, 对于大的FFT, 这种定标可能会使信息丢失。 另一种方法是动态定标, 即在输出溢出时再除以2。在 这种情况下, 会在这个过程中指定一个变量, 每定标一 次变量的值加1, 计算结束后根据变量的值把结果乘以 一个系数。动态定标的信噪比最好, 但会增加FFT循环 次数。
第7章 应用程序设计 内容提要： 7.1 定标与溢出处理 7.2 基础算术运算 7.3 FIR滤波器 7.4 IIR滤波器 第7章 应用程序设计 内容提要： 7.1 定标与溢出处理 7.2 基础算术运算 7.3 FIR滤波器 7.4 IIR滤波器 7.5 快速傅里叶变换(FFT) 7.6 DSPLIB的使用
7.2 基础算术运算 7.2.1 加减运算 7.2.2 乘法运算 7.2.3 除法运算 7.2.4 小数乘法
7.2.1 加减运算 在数字信号处理中，加减运算是常见的算术运算。 一般使用16位或32位加减运算，数值分析、浮 点运算和其它操作可能需要32位以上的运算。 C55x有直接完成16位或32位加减运算的指令， 但没有能直接完成多字加减运算的指令。要进行 多字加减运算，需要通过编程方法实现。
MOV40 dbl(Lmem),ACx ADD dbl(Lmem),ACx 以下指令可在单周期内完成32位加法运算： MOV40 dbl(Lmem),ACx ADD dbl(Lmem),ACx 64位的高32位加法要考虑低32位加法产生的进位，高32位使用以下指令： ADD uns(Smem), CARRY, ACx 以下指令可在单周期内完成32位减法运算： MOV40 dbl(Lmem), ACx SUB dbl(Lmem), ACx 64位的高32位减法要考虑低32位减法产生的借位，使用以下指令： SUB uns(Smem), BORROW,ACx 可进行单独32位或64位的低32位加 32~47位 48~63位 ADD (Smem-1)<< #16, ACx 可进行单独32位或64位的低32位加 SUB (Smem-1)<< #16, ACx
例7-1，64位加法运算。文件名为：add64.asm。 .mmregs为存储器映射寄存器定义符号名. 与对所有存储器映射寄存器执行.set功能一样。C54x 汇编器(手册)支持.mmregs指令。 .mmregs .model call=c55_std .model mem=large C函数调用时参量与寄存器(10个)的排列关系按表6-11执行(默认: c55_std) 。 c55_compat只能传递3个参数到寄存器 。 指针分配 AR1 -> X3(偶地址) X2 X1 X0 AR2 -> Y3(偶地址) Y2 Y1 Y0 AR3 -> W3(偶地址) W2 W1 W0 1000h ;***************** ; 64位加法 ; ; X3 X2 X1 X0 ; + Y3 Y2 Y1 Y0 ; --------------------- ; W3 W2 W1 W0 地址值依次增加 1004h 1008h
MOV40 dbl(*AR1(#2)), AC0 ;AC0 = X1 X0 .sym-- COFF Symbolic Debugging Directives已废弃 .sect ".text" .align 4 .global start ;.sym start,start, 36, 2, 0 start: MOV #0100h,AR1 MOV #0104h,AR2 MOV #0108h,AR3 L1: 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 Y2 + CARRY ADD *AR2<< #16, AC0 ;AC0= X3 X2 +Y3 Y2 + CARRY MOV AC0, dbl(*AR3) ;保存 W3 W2. B L1 The .sym directive defines a global variable, a local variable, or a function. Several parameters allow you to associate various debugging information with the variable or function. 格式: .sym name, value[, type, storage class, size, tag, dims] 用数表示对应C变量类型,占存储单元位数 标号,地址,表示函数,外部,不占存储 参考文献: TMS320C28X系列DSP指令和编程指南-附录B 符号调试伪指令.sym(需看附录A:难懂) 茶苑老吴的博客:”关于TMS320C55x的汇编语言中的.sym伪指令 “(此文献非常详细明了)
MOV40 dbl(*AR1(#2)), AC0 ;AC0 = X1 X0 AR1 -> X3(偶地址) X2 X1 X0 AR2 -> Y3(偶地址) Y2 Y1 Y0 AR3 -> W3(偶地址) W2 W1 W0 1000h 地址值依次增加 1004h 1008h 指针分配 .sect ".text" .align 4 .global start ;.sym start,start, 36, 2, 0 start: MOV #0100h,AR1 MOV #0104h,AR2 MOV #0108h,AR3 L1: 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=X3X2+00Y2 +CARRY ADD *AR2<< #16, AC0 ;两行完成AC0=X3X2+Y3Y2+CARRY MOV AC0, dbl(*AR3) ;保存W3 W2. B L1 1002h 1005h 1006h 100Ah 该指令执行后AR1的值不变,该指令寻址:(AR1)+2 无带进位的32位加指令 加法AC0=AC0+Y300 此指令不是带进位加, 也不需要, 因X2加Y2的进位就在AC0(16)中.
MOV40 dbl(*AR1(#2)), AC0 ;AC0 = X1 X0 .sect ".text" .align 4 .global start ;.sym start,start,36,2,0 start: MOV #0100h,AR1 MOV #0104h,AR2 MOV #0108h,AR3 L1: 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=X3X2+00Y2 +CARRY ADD *AR2<< #16, AC0 ;两行完成AC0=X3X2+Y3Y2+CARRY MOV AC0, dbl(*AR3) ;保存 W3 W2. B L1 存储单元0100h~0107h中的数据, 例题中没有给出,可添加数据测试: RPT #7 MOV #8888H,*AR1+ MOV #1288H,*(#0100h) MOV #1288H,*(#0104h) MOV #0100h,AR1 无带进位的32位加指令
例7-2 , 64位减法运算程序。文件名为: sub64.asm .mmregs .model call=c55_std .model mem=large 指针分配 AR1 -> X3(偶地址) X2 X1 X0 AR2 -> Y3(偶地址) Y2 Y1 Y0 AR3 -> W3(偶地址) W2 W1 W0 1000h ;***************** ; 64位减法 ; ; X3 X2 X1 X0 ; − Y3 Y2 Y1 Y0 ; --------------------- ; W3 W2 W1 W0 地址值依次增加 1004h 1008h
MOV40 dbl(*AR1(#2)),AC0 ;AC0=X1X0 SUB dbl(*AR2(#2)),AC0 ;AC0=X1X0-Y1Y0 AR1 -> X3(偶地址) X2 X1 X0 AR2 -> Y3(偶地址) Y2 Y1 Y0 AR3 -> W3(偶地址) W2 W1 W0 1000h 地址值依次增加 1004h 1008h 指针分配 .sect ".text" .align 4 .global start ;.sym start, start,36,2,0 start: MOV #0100h,AR1 MOV #0104h,AR2 MOV #0108h,AR3 L1: 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 *AR2<<#16,AC0 ;AC0=X3X2-Y3Y2-BORROW MOV AC0,dbl(*AR3) ;保存 W3W2. B L1 除加法换成减法, CARRY换成BORROW外, 步骤相同.
7.2 基础算术运算 7.2.1 加减运算 7.2.2 乘法运算 7.2.3 除法运算 7.2.4 小数乘法
7.2.2 乘法运算 C55x提供了硬件乘法器，16位乘法可在一个 指令周期内完成。 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位操作数 X0 ; W3 W2 W1 W0 64位结果 AR1 -> Y1 ; Y0 ; 入口条件: AR2 -> W0 ; SXMD = 1 (允许符号扩展) W1 ; SATD = 0 (不做饱和处理) W2 ; FRCT = 0 (关小数模式) W3 ;;限制条件: 延迟链和输入序列必须指定为长字类型。 ;*************************************************************
例7-3, 32位整数乘法运算。 .sect ".text" .align 4 .global start 有符号数的数据存储: X1 X0 32位 × Y1 Y0 32位 无符号乘 X0 ×Y0 有符号与无符号乘 X1×Y0 有符号与无符号乘 X0×Y1 有符号乘 X1×Y1 结果累加 W3 W2 W1 W0 64位 .sect ".text" .align 4 .global start ;.sym start,start,36,2,0 start: MOV #0100h, AR0 MOV #0102h, AR1 MOV #0104h, AR2 BSET SXMD; 入口条件 BCLR SATD BCLR FRCT 32位 32位 32位 指针分配: AR0 -> X1 X0 AR1 -> Y1 Y0 AR2 -> W0 W1 W2 W3 0100h 地址值依次增加 0102h 0104h 存储单元0100h~0103h中的数据, 例题中没有给出,可添加数据测试: RPT #3 MOV #8888H,*AR0+ MOV #1288H,*(#0100h) ;MOV #1288H,*(#0102h) ;正数1288使乘积为正, 否则一正一负,乘积为负 MOV #0100h,AR0
MPYM uns(*AR0-),uns(*AR1),AC0 ;ACO=X0*Y0 例7-3, 32位整数乘法运算。 AR2->W0 W1 W2 W3 0104h AR0->X1 X0 AR1->Y1 Y0 0100h 0102h 指针分配 L1: AMAR *AR0+ ;AR0指向X0 ||AMAR *AR1+ ;AR1指向Y0 MPYM uns(*AR0-),uns(*AR1),AC0 ;ACO=X0*Y0 MOV AC0,*AR2+ ;保存W0 (是AC0的低16位) MACM *AR0+,uns(*AR1-),AC0>>#16,AC0 ;AC0=X0*Y0>>16+X1*Y0 MACM uns(*AR0-),*AR1,AC0 ;AC0=X0*Y0>>16+X1*Y0+X0*Y1 MOV AC0,*AR2+ ;保存W1 MACM *AR0,*AR1,AC0>>#16,AC0 ;AC0=AC0>>16+X1*Y1 MOV AC0,*AR2+ ;保存W2 MOV HI(AC0),*AR2 ;保存W3 B L1 1005h 1001h 1003h 带符号乘 ;即累加前AC0的高16位 X1 X0 32位 × Y1 Y0 32位 X0 ×Y0 X1×Y0 X0×Y1 X1×Y1 累加W3 W2 W1 W064位
7.2 基础算术运算 7.2.1 加减运算 7.2.2 乘法运算 7.2.3 除法运算 7.2.4 小数乘法
7.2.3 除法运算 C55x没有提供硬件除法器，也没有提供专门的除法指令，要实现除法运算需借助于条件减法指令SUBC和重复指令RPT。 7.2.3 除法运算 C55x没有提供硬件除法器，也没有提供专门的除法指令，要实现除法运算需借助于条件减法指令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-0)为商, AC0(31-16)为余数。
BCLR SXMD ; 清零SXMD(关闭符号扩展) MOV *AR0,AC0 ; 把被除数放入AC0 .sect ".text" .align 4 .global start ;.sym start,start,36,2,0 start: MOV #0100h,AR0 MOV #0101h,AR1 MOV #0102h,AR2 MOV #0103h,AR3 L1: BCLR SXMD ; 清零SXMD(关闭符号扩展) MOV *AR0,AC0 ; 把被除数放入AC0 RPT #(16-1) ; 执行subc 16次 SUBC *AR1,AC0,AC0 ; AR1指向除数 MOV AC0,*AR2 ; 保存商 MOV HI(AC0),*AR3 ; 保存余数 B L1 例7-4, 无符号16位除16位整数除法。 指针分配 AR0->被除数 AR1->除数 AR2->商 AR3->余数 0100h 0102h 0101h 0103h 存储单元0100h~0101h中的数据, 例题中没有给出,可添加数据测试: MOV #0608H,*(#0100h) MOV #020H,*(#0101h)
例7-5，无符号32位除16位整数除法。 文件名为：udiv32o16.asm。 .mmregs .model call=c55_std .model mem=large ;*********************************************** ; 指针分配 ; AR0->被除数高位 ; 被除数低位 ; AR1->除数 ; AR2->商高位 ; 商低位 ; AR3->余数 ;注： ;无符号除法，被除数为32位，除数为16位 ;关闭符号扩展，被除数、除数均为正数 ;第一次除法之前, 把被除数高16位存入AC0(15-0) ;第一次除法之后, 商的高位在AC0(15-0)中, 需转存到AR2寻址单元 ;第二次除法之前,余数在AC0(31-16),把被除数低16位存入AC0(15-0) ;第二次除法之后, AC0(15-0)为商的低16位, AC0(31-16)为余数 ;************************************************
L1: BCLR SXMD ;清零SXMD(关闭符号扩展) MOV *AR0+,AC0 ;把被除数高位存入AC0 指针分配 AR0->被除数高位 被除数低位 AR1->除数 AR2->商高位 商低位 AR3->余数 0100h 0104h 0102h 0106h .sect ".text" .align 4 .global start start: MOV #0100h,AR0 MOV #0102h,AR1 MOV #0104h,AR2 MOV #0106h,AR3 L1: 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 存储地址(08h)装入AR4 MOV AC0,*AR2+ ;保存商的高位 MOV *AR0+,*AR4 ;把被除数低位装入AC0_L(08h) RPT #(16-1) SUBC *AR1,AC0,AC0 ;执行subc 16次 MOV AC0,*AR2+ ;保存商的低位 MOV HI(AC0),*AR3 ;保存余数 BSET SXMD ;置位SXMD (打开符号扩展) B L1 或者用MOV *AR0+,*(#AC0L) ;不能用MOV *AR0+,AC0(高16位清0)
例7-6， 带符号16位除16位整数除法。 文件名为：sdiv16o16.asm。 .mmregs .model call=c55_std .model mem=large .sect ".text" .align 4 ;***************************************** ; 指针分配 ; AR0-> 被除数 ; AR1-> 除数 ; AR2-> 商 ; AR3-> 余数 ; 注： ; 带符号除法，被除数为16位，除数为16位 ;打开符号扩展，被除数、除数可为负数 ;除法运算之前，商的符号存入AC0 ;除法运算之后，商存入AC1(15-0), 余数存入AC1(31-16)
指针分配 AR0->被除数 AR1->除数 AR2->商 AR3->余数 0100h 0102h 0101h 0103h .global start start: MOV #0100h,AR0 MOV #0101h,AR1 MOV #0102h,AR2 MOV #0103h,AR3 L1: BSET SXMD ; 置位号SXMD (打开符扩展) MPYM *AR0,*AR1,AC0 ; 计算期望得到的商的符号 MOV *AR1,AC1 ; 把除数存入AC1 ABS AC1,AC1 ; 求绝对值，|除数| MOV AC1,*AR2 ; 暂时保存 |除数| MOV *AR0,AC1 ; 把被除数存入 AC1 ABS AC1,AC1 ; 求绝对值，|被除数| RPT #(16-1) ; 执行subc 16次 把除数 取绝对值 把被除数 取绝对值 无符号数除法 SUBC *AR2,AC1,AC1 ;AR2 -> |除数| MOV HI(AC1),*AR3 ; 保存余数 MOV AC1,*AR2 ; 保存商(AC1的低16位,是正数) SFTS AC1,#16 ;对正数商左移16位: 使符号位在最高位(31) NEG AC1,AC1 ; 对商求反2补码(商变成负数) 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-> 被除数高半部分(高16位)(NumH)(偶地址) ; 被除数低半部分(低16位)(NumL) ; AR1-> 除数(Den) ; AR2-> 商的高半部分(高16位)(QuotH)(偶地址) ; 商的低半部分(低16位)(QuotL) ; AR3-> 余数(Rem) ; 注： ; 带符号除法，被除数为32位，除数为16位 ;打开符号扩展，被除数、除数可为负数 ;除法运算之前，期望的商的符号存入AC0 ;第一次除法运算之前，把被除数的高16位分存入AC1 ;第一次除法运算之后，把商的高16位存入AC1(15-0) ;第二次除法运算之前，把被除数的低16位存入AC1 ;第二次除法运算之后，把商的低16位存入AC1(15-0),余数存入AC1(31-16) ;*******************************************************
.sect ".text" .align 4 .global start .sym start,start,36,2,0 start: MOV #0100h,AR0 MOV #0102h,AR1 MOV #0104h,AR2 MOV #0106h,AR3 MOV #0108h,AR4 L1: BSET SXMD ; 置位SXMD (打开符号扩展) MPYM *AR0,*AR1,AC0 ; 除法结果的符号位( NumH x Den ) MOV *AR1,AC1 ; AC1 = Den (16位除数) ABS AC1,AC1 ; AC1 = abs(Den) (除数取绝对值) MOV AC1,*AR3 ; Rem = abs(Den) (暂存除数绝对值) MOV40 dbl(*AR0),AC1 ; AC1 = NumH NumL (32位被除数) ABS AC1,AC1 ; AC1 = abs(Num) (被除数取绝对值) MOV AC1,dbl(*AR2) ;QuotH = abs(NumH) (暂存被除数绝对值) ;QuotL = abs(NumL) 指针分配 AR0->被除数高位(NumH) 被除数低位(NumL) AR1->除数 AR2->商高位(QuotH) 商低位(QuotL) AR3->余数 AR4->AC1_L存储地址 0100h 0104h 0102h 0106h 0108h 11h
MOV *AR2,AC1 ; AC1 = QuotH (被除数取绝对值高16位) RPT #(15-1) ; 执行subc 15次 SUBC *AR3,AC1,AC1 ; (被除数取绝对值高16位除以除数) SUBC *AR3,AC1,AC1 ; 最后一次执行subc || MOV #11,AR4 ; 把AC1_L存储地址装入AR4 MOV AC1,*AR2+ ; 保存 QuotH (存储结果商的高16位) MOV *AR2,*AR4 ;AC1_L=QuotH(应该是被除数绝对值低16位) RPT #(16-1) ;执行subc 16次 SUBC *AR3,AC1,AC1 MOV AC1,*AR2- ; 保存 QuotL MOV HI(AC1),*AR3 ; 保存 Rem BCC skip,AC0>= #0 ; 如果实际结果应该为正数,跳到skip. MOV40 dbl(*AR2),AC1 ; 否则,对商取反. NEG AC1,AC1 MOV AC1,dbl(*AR2) skip: B L1 (QuotL)的存储单元内容
7.2 基础算术运算 7.2.1 加减运算 7.2.2 乘法运算 7.2.3 除法运算 7.2.4 小数乘法
7.2.4 小数乘法 在定点DSP的某些应用中，整数运算很难满足要求。这是因为它自身存在缺陷： 两个16位整数相乘, 乘积总是“向左增长” (即小数点左侧的位数增加), 这意味着多次相乘后, 乘积将很快超出定点器件的数据范围。 保存32位乘积到存储器, 要占用2个CPU周期和2个字的存储器空间。 由于乘法器都是16位相乘, 因此将32位乘积再作为乘法器的输入时就显得较繁琐, 不能胜任递归运算。
7.2.4 小数乘法 为了克服这些缺陷，在实际应用中更多采用的是 小数运算。小数运算具有如下优点： 乘积总是“向右增长”。这就意味着超出定点器 件数据范围的将是不太感兴趣的部分。 既可以存储32位乘积, 也可近似存储高16位乘积， 这就允许用较少的资源保存结果。 可以用于递归运算。 小数部分越靠近小数点, 其权值越大。
例7-8, 两个Q31格式的有符号小数相乘, 得到一个Q31格式结果。 文件名为: mpyQ31.asm。 需要注意的是： 两个带符号小数相乘, 所得乘积带有2个符号位。 为解决冗余符号位的问题, 需在程序中设定状态寄 存器ST1中的FRCT(小数方式)为1, 这样当乘法器 将结果传送至累加器时就会自动左移1位。
;******************************************* ;操作数取自数据存储器，运算结果送回数据存储器。 例7-8, 两个Q31格式有符号小数相乘。 X1 X0 Q31 × Y1 Y0 Q31 X0 ×Y0 X1×Y0 X0×Y1 X1×Y1 累加 W1 W0 w’1 w’0(忽略) .mmregs .model call=c55_std .model mem=large ;******************************************* ;操作数取自数据存储器，运算结果送回数据存储器。 ; 数据存储: 指针分配: ; X1 X0 Q31操作数 AR0 -> X1 ; ×Y1 Y0 Q31操作数 X0 ; W1 W0 Q31结果 AR1 -> Y1 ; Y0 ; 入口条件: AR2 ->W0(偶地址) ; SXMD = 1 (允许符号扩展) W1 ; SATD = 0 (不做饱和处理) ; FRCT = 1 (运算结果左移一位) ;限制条件：W1被指定为偶地址
MPYM uns(*AR0-),*AR1+,AC0 ;AC0=X0*Y1 AR0->X1 Q31乘数高位 X0 Q31乘数低位 AR1->Y1 Q31乘数高位 Y0 Q31乘数低位 AR2->W1 Q31乘结果高位 W0 Q31乘结果低位 0100h 0104h 0102h .sect ".text" .align 4 .global start ;.sym start,start,36,2,0 start:MOV #0100h,AR0 MOV #0102h,AR1 MOV #0104h,AR2 BSET SXMD;入口条件 BCLR SATD BSET FRCT L1: AMAR *AR0+ ;AR0指向X0 MPYM uns(*AR0-),*AR1+,AC0 ;AC0=X0*Y1 MACM *AR0,uns(*AR1-),AC0 ;AC0=X0*Y1+X1*Y0 MACM *AR0,*AR1,AC0>>#16,AC0;AC0=AC0>>16+X1*Y1 MOV AC0,dbl(*AR2) ;保存W1W0 B L1 1001h Y1 Y0 Q31 × X1 X0 Q31 X0 ×Y0(忽略,不乘) X0×Y1 X1×Y0 X1×Y1 累加 W1 W0 w’1 w’0(忽略) (忽略,不乘) 存储单元0100h~0103h中的数据, 例题中没有给出,可添加数据测试: RPT #3 MOV #8888H,*AR0+ MOV #1288H,*(#0100h) ;MOV #1288H,*(#0102h) ;正数1288使乘积为正, 否则一正一负,乘积为负 MOV #0100h,AR0
第7章作业 7.7 (仿真调试)利用第3章中C语言工程所用命 令文件Exam3_2.cmd ，在CCS下建立关于例 7-8的语言工程，并对工程进行构建、调试。 7.1, 7.2, 7.6 (参考习题看书)
第7章 应用程序设计 内容提要： 7.1 定标与溢出处理 7.2 基础算术运算 7.3 FIR滤波器 7.4 IIR滤波器 第7章 应用程序设计 内容提要： 7.1 定标与溢出处理 7.2 基础算术运算 7.3 FIR滤波器 7.4 IIR滤波器 7.5 快速傅里叶变换(FFT) 7.6 DSPLIB的使用
7.3 FIR滤波器 数字滤波器是DSP的基本应用，有2种基本类型： 有限冲激响应滤波器FIR 无限冲激滤波器IIR 本节主要讨论FIR滤波器的DSP实现方法，有关 IIR滤波器的实现将在下一节中介绍。 7.3.1 FIR滤波器的基本结构 7.3.2 FIR滤波器的C语言编程实现 7.3.3 FIR滤波器的汇编语言编程实现
7.3.1 FIR滤波器的基本结构 一个FIR滤波器的输出序列和输入序列之间的关系，满足差分方程： 传递函数为
FIR滤波器的结构： FIR滤波器的单位冲激响应是一个有限长序列。若为实数，且满足偶对称或奇对称的条件，则FIR滤波器具有线性相位特性。 偶对称线性相位FIR滤波器的差分方程为：
7.3.2 FIR滤波器的C语言编程实现 例7-9，直接型FIR滤波器的C语言编程实现。 /****************************************** * fir.c 该程序用于实现FIR滤波器 * L——滤波器的阶数 * b[i]——滤波器的系数,i=0,1,…,L-1 * x[i]——输入信号向量,i=0,1,…,L-1;x[0]对应于当前值,x[1]对应于上一采样值 * x_in —— 输入信号的当前值 * y_out—— 输出信号的当前值 ********************************************/
float fir(float x_in, float *x,float *b,int L) { float y_out; int i; /* -------------------------------------------------------- */ /* 把上一个采样时间的输入信号向量延迟一个单元,得到当前 采样时间的输入信号向量 */ for(i=L-1;i>0;i--) x[i] = x[i-1]; } x[0]=x_in; /* ------------------------------------------------------- */ /*完成FIR滤波 */ y_out = 0.0; for(i=0;i<L;i++) y_out = y_out + b[i]*x[i]; return y_out;
直接型FIR滤波器的实现涉及到两个基本操作，一个是输入信号向量与滤波器系数向量的内积计算，另一个是输入信号向量的更新处理。 在每个采样周期信号缓冲器都要更新一次，最老的采样被抛弃，而其他的信号则向缓冲器的右方移动一个单元，一个新的采样被插入存储单元，并被标记。 ☼ 若此操作过程不用DSP硬件完成, 那么它需很多时间。
7.3.3 FIR滤波器的汇编语言编程实现 处理信号缓冲器的最有效方法，是把信号采样加 载到循环缓冲器中。 在循环缓冲器中，采取数据保持固定、反时针方 向移动地址的方式，代替保持缓冲器地址固定且 正方向移动数据。 信号采样的起点由指针x(n) 指定，其它诸采样则 沿着顺时针方向，从起点开始依次顺序加载。
(a)信号循环缓冲区 (b)系数循环缓冲区 图7-4 FIR滤波器的循环缓冲区
例7-10，FIR滤波器的C55x汇编语言实现。 (1) 主程序 fir_test.c /* fir_test.c */ #include"math.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*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 coeff[L]={ -26,-13,14,36,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 in[L]; /* input buffer */ int out[L]; /* Output buffer */
main() { unsigned int i; float signal; unsigned int n=0; int index=0; for(i=0;i<L;i++) { in[i]=0; out[i]=0; } while(1) signal = a1*cos((float)w1*n); signal += a2*cos((float)w2*n); signal += a3*cos((float)w3*n); n++; in[index] = (int)((0x7fff*signal)+0.5); out[index] = fir(in,coeff,L,index); index--; if(index==-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_55 pshm ST3_55 or #0x340,mmap(ST1_55) ;设置FRCT,SXMD,SATD bset SMUL ;置位SMUL
mov mmap(AR0),BSA01 ;AR0=输入信号循环缓冲区的起始地址 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 ;设置外部循环次数为L-1 mpym *AR0+,*AR2+,AC0 ;执行第一次运算 || rpt CSR ;启动循环 macm *AR0+,*AR2+,AC0 mov hi(AC0),T0 ; 用Q15格式存放结果 popm ST3_55 ; 恢复ST1, ST2和 ST3 popm ST2_55 popm ST1_55 ret .end
第7章 应用程序设计 内容提要： 7.1 定标与溢出处理 7.2 基础算术运算 7.3 FIR滤波器 7.4 IIR滤波器 第7章 应用程序设计 内容提要： 7.1 定标与溢出处理 7.2 基础算术运算 7.3 FIR滤波器 7.4 IIR滤波器 7.5 快速傅里叶变换(FFT) 7.6 DSPLIB的使用
7.4 IIR滤波器 IIR滤波器的优点：结构简单，运算量小，可以用 较少的阶数获得很高的选择性。 7.4.3 IIR滤波器的C语言实现 7.4.4 IIR滤波器的汇编语言实现
7.4.1 二阶IIR滤波器的结构 二阶IIR滤波器，又称为二阶基本节，分为直接型、标准型和变换型。 二阶IIR滤波器传递函数为：
图7-5 二阶IIR滤波器的直接型Ⅰ的实现
图7-6 的信号流图
由于直接型Ⅱ对于给定的传递函数具有最小可能的延迟数、加法器数和乘法器数，所以被称为标准型。 图7-7 二阶IIR滤波器的直接型Ⅱ的实现 由于直接型Ⅱ对于给定的传递函数具有最小可能的延迟数、加法器数和乘法器数，所以被称为标准型。
7.4.2 高阶IIR滤波器的结构 高阶IIR滤波器的差分方程和系统函数分别为：
图7-8 高阶IIR滤波器的直接型Ⅱ实现（L=M+1）
图7-9 高阶IIR滤波器的串联型结构 图7-10 第k个二阶节
对于 ，有 其中， 注意：
7.4.3 IIR滤波器的C语言实现 例7-11,采用二维数组编写的IIR滤波器C语言程序。 a[i][k]（i=1,2）和b[i][k] （i=0,1,2）分别为第k（k=0,1,K）个二阶IIR滤波器环节的系数。 例7-11,采用二维数组编写的IIR滤波器C语言程序。 temp=xin; /*xin为IIR滤波器的输入 */ for(k=0;k<N_IIR;k++) /*N_IIR为IIR滤波器二阶节的个数 { w[0][k]=temp-a[1][k]*w[1][k]-a[2][k]*w[2][k]; /*这里,temp为本二阶节的输入,也是上一级二阶节的输出*/ temp=b[0][k]*w[0][k]+b[1][k]*w[1][k]+b[2][k]*w[2][k]; /*这里,temp为本二阶节的输出,也是下一级二阶节的输入*/ w[2][k]=w[1][k]; w[1][k]=w[0][k]; } xoutput=temp;/*xoutput为IIR滤波器的输出 */ w[i][k] (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 *,unsigned 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 _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 ; 设置外循环计数器 sub #1,T1,T0 ; 二阶节个数 - 1 mov T0,BRC1 ; 设置内循环计数器 mov T1,T0 ;设置循环缓冲区大小 sfts T0,#1 mov mmap(T0),BK03 ;BK03=2*二阶节个数 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 ;用于左移 || 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(AC0<<T0)),*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 ST3_55 ;恢复ST1,ST2,ST3 popm ST2_55 popm ST1_55 ret .end */
图7-11 IIR滤波器系数和信号缓冲区配置
第7章 应用程序设计 内容提要： 7.1 定标与溢出处理 7.2 基础算术运算 7.3 FIR滤波器 7.4 IIR滤波器 第7章 应用程序设计 内容提要： 7.1 定标与溢出处理 7.2 基础算术运算 7.3 FIR滤波器 7.4 IIR滤波器 7.5 快速傅里叶变换(FFT) 7.6 DSPLIB的使用
7.5 快速傅里叶变换FFT FFT算法原理 库利一图基算法 FFT算法的实现
7.5.1 FFT算法原理 快速傅里叶变换（FFT）是离散傅里叶变换（DFT）的一种快速算法。通过FFT算法，DFT的计算量大大减少，运算时间缩短1~2个数量级。 DFT的变换公式为 其中 为旋转因子。
FFT之所以减少运算量, 主要利用旋转因子的3点特性: 对称性 周期性 可约性 利用这些特性可以使DFT运算中有些项进行合并，将长序列的DFT分解为短序列的DFT。 DFT从算法上分为按时间抽选（DIT）和按频率抽选(DIF）。 基2的DIT又被称为库利一图基算法。基2的DIF又称为桑德—图基算法。
7.5 快速傅里叶变换FFT FFT算法原理 库利一图基算法 FFT算法的实现
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，将这些二进制数前后倒转，即得到进行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 快速傅里叶变换FFT FFT算法原理 库利一图基算法 FFT算法的实现
7.5.3 FFT算法的实现 为了叙述简单，本书给出采用C语言编写的FFT程序，相应的汇编程序请读者自行完成。 例7-13，基2DIT FFT算法的C语言实现。 （1）主程序：fft_test.c #include <math.h> #include "fcomplex.h" /*包含浮点复数结构体定义头文件fcomplex.h*/ extern void bit_rev(complex *,unsigned int);/*位反转函数声明*/ extern void fft(complex *,unsigned int,complex *,unsigned int); /* fft函数声明 */
extern void generator(int *,unsigned int) #define N 128 /* FFT的数据个数 */ #define M 7 /* M=log2(N) */ #define PI 3.1415926 complex X[N]; /* 说明输入信号数组,为复数 */ complex W[M]; /* 说明旋转因子数组e^(-j2PI/N),为复数 */ complex temp; /* 说明临时复数变量 */ float xin[N]; float spectrum[N]; /* 说明功率谱信号数组,为实数 */ float re1[N],im1[N]; /* 说明临时变量数组,为实数 */
void main() { unsigned int i,j,L,LE,LE1; /*----------------------------------------------------- */ /* 产生旋转因子表 */ for (L=1; L<=M; L++) LE=1<<L; /* 子FFT中的点数LE=2^L */ LE1=LE>>1; /* 子FFT中的蝶形运算数目*/ W[L-1].re = cos(PI/LE1); W[L-1].im = -sin(PI/LE1); } /* ----------------------------------------------------- */ generator(xin,N);
for (;;) { for (i=0; i<N; i++) /* 构造输入信号样本 */ X[i].re =xin[i]; X[i].im = 0; /* 复制到参考缓冲器 */ re1[i] = X[i].re; im1[i] = X[i].im; } /* 启动 FFT */ bit_rev(X,M); /* 以倒位次序排列X[] */ fft(X,M,W,1); /* 执行 FFT */ /* 计算功率谱，验证FFT结果 */ temp.re = X[i].re*X[i].re; temp.im = X[i].im*X[i].im; spectrum[i] = (temp.re + temp.im)*4;
（2）浮点复数基2 DIT FFT函数：fft_float.c #include "fcomplex.h" void fft(complex *X,unsigned int M,complex *W,unsigned int SCALE) { complex temp; /* 复变量临时存储器 */ complex U; /* 旋转因子W^k */ unsigned int i,j; unsigned int id; /* 蝶形运算中下位节点的序号 */ unsigned int N=1<<M; /* FFT 的点数*/ unsigned int L; /* FFT 的级序号 */ unsigned int LE; /* L级子FFT的点数 */ unsigned int LE1; /* L级子FFT蝶形运算的个数 */ float scale; scale = 0.5; for(L=1;L<=M;L++) LE=1<<L; LE1=LE>>1; U.re = 1.0; U.im = 0.;
for(j=0;j<LE1;j++) { for(i=j;i<N;i+=LE) /*进行蝶形计算*/ id=i+LE1; temp.re=(X[id].re*U.re-X[id].im*U.im)*scale; temp.im=(X[id].im*U.re+X[id].re*U.im)*scale; X[id].re=X[i].re*scale-temp.re; X[id].im=X[i].im*scale-temp.im; X[i].re=X[i].re*scale+temp.re; X[i].im=X[i].im*scale+temp.im; } /*递推计算W^k*/ temp.re=U.re*W[L-1].re-U.im*W[L-1].im; U.im=U.re*W[L-1].im+U.im*W[L-1].re; U.re=temp.re;
（3）位反转函数：bit_rev.c #include "fcomplex.h" void bit_rev(complex *X,unsigned int M) { complex temp; unsigned int i,j,k; unsigned int N=1<<M; /* FFT 的点数*/ unsigned int N2=N>>1; for (j=0,i=1; i<N-1; i++) k=N2; while(k<=j){ j-=k; k>>=1; } j+=k; if(i<j) { temp=X[j]; X[j]=X[i];X[i]=temp;} }
（4）信号发生器函数：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;i<N;i++) X[i]=a1*cos((float)w1*i); }
（5）复数结构定义头文件：fcomplex.h struct cmpx { float re; float im; }; typedef struct cmpx complex;
第7章 应用程序设计 内容提要： 7.1 定标与溢出处理 7.2 基础算术运算 7.3 FIR滤波器 7.4 IIR滤波器 第7章 应用程序设计 内容提要： 7.1 定标与溢出处理 7.2 基础算术运算 7.3 FIR滤波器 7.4 IIR滤波器 7.5 快速傅里叶变换(FFT) 7.6 DSPLIB的使用
7.6 DSPLIB的使用 7.6.1 DSPLIB简介 7.6.2 CCS下DSPLIB的安装 7.6.3 DSPLIB的数据类型
7.6.1 DSPLIB简介 DSPLIB是TI公司提供的一个优化的DSP函数库，它包括50多个可在C程序中调用的汇编优化的通用信号处理程序。 采用DSPLIB 可以获得执行速度显著快于采用标准C语言编写的相应程序，大大减少DSP应用系统的开发时间。 DSPLIB 为免费软件，可到下列网址下载： http://www.ti.com/sc/docs/products/dsp/c5000/index.Htm。
7.6.2 CCS下DSPLIB的安装 DSPLIB库的相关文件被压缩在一个名字为 55xdsplib.exe的压缩包中 安装方法: (1)运行55xdsplib.exe，得到名字为 dsplib_2.40.00的文件夹 (2)将该文件夹拷贝到CCS5.4下的相应子目录下。 (以下设此子目录为D:\ccs54\ccsv5 \tools \compiler\c5500_5.4.1)
7.6.2 CCS下DSPLIB的安装 dsplib_2.40.00文件夹中的主要内容 55x_src (文件夹) : 汇编源程序文件。 examples(文件夹)：调用DSPLIB函数的范例。 include(文件夹)：DSPLIB函数的头函数文件。 55xdsp.lib :小存储模式下的dsplib程序库。 55xdspx.lib :大存储模式下的dsplib程序库。 blt55x.bat :用于重新生成基于55xdsp.src的 55xdsp.lib。 blt55xx.bat :用于重新生成基于55xdspx.src 的55xdspx.lib。 redme.txt：关于DSPLIB库的简单说明。
7.6.3 DSPLIB的数据类型 DSPLIB新定义了以下数据类型类型： Q.15(DATA): 一个Q.15操作数表达为short数据类型(16 bit), 在dsplib.h 头文件中预定义为DATA。 Q.31(LDATA): 一个Q.31操作数表达为long数据类型(32 bit), 在dsplib.h 头文件中预定义为LDATA。 Q.3.12: 含有3个整数位和12个小数位。 除非特别说明，DSPLIB采用Q15小数数据类型。 为了提高效率，DSPLIB库函数采用的操作数通常为向量形式。标量可视为维数为1的向量。
7.6.4 DSPLIB的参量 向量操作数中的各元素连续地存放在存储器空间中。 复数元素以实-虚格式存放。除非特别注明，源操作数 与目的操作数使用同一存储器空间，以节约内存，称 之为同址（In-place）计算。 DSPLIB库函数参量采用的符号约定： x,y：输入数据向量。 r ： 输出数据向量。 nx,ny,nr：向量 x,y,r的大小。若nx = nr = nr, 则只使用nx。 h: 滤波器系数向量。 Nh: h的大小。
7.6.5 DSPLIB函数简介 DSPLIB函数可以分为以下8类: FFT函数 滤波和卷积 自适应滤波 相关 数学 三角 矩阵 其他
7.6.5 DSPLIB函数简介 FFT函数 void cfft (DATA *x, ushort nx, type):基2复数FFT。 void cifft (DATA *x, ushort nx, type):基2复数IFFT。 void cbrev (DATA *x, DATA *r, ushort n):复数位反转。 void rfft (DATA *x, ushort nx, type):基2实数FFT。 void rifft (DATA *x, ushort nx, type):基2实数IFFT。 滤波和卷积函数 ushort fir (DATA *x, DATA *h, DATA *r, DATA *dbuffer,ushort nx, ushort nh) :FIR直接型。 ushort fir2 (DATA *x, DATA *h, DATA *r, DATA *dbuffer,ushort nx, ushort nh): FIR直接型(采用DUAL–MAC优化)。 ushort firs (DATA *x, DATA *h, DATA *r, DATA *dbuffer,ushort nx, ushort nh2) :对称FIR直接型。 ushort cfir (DATA *x, DATA *h, DATA *r, DATA *dbuffer,ushort nx, ushort nh):复数FIR直接型。
7.6.5 DSPLIB函数简介 ushort convol (DATA *x, DATA *h, DATA *r, ushort nr,ushort nh):卷积。 ushort convol1 (DATA *x, DATA *h, DATA *r, ushort nr,ushort nh):卷积 (采用DUAL–MAC优化)。 ushort convol2 (DATA *x, DATA *h, DATA *r, ushort nr,ushort nh):卷积 (采用DUAL–MAC优化)。 ushort iircas4 (DATA *x, DATA *h, DATA *r, DATA *dbuffer,ushort nbiq, ushort nx):IIR级联直 接2型，每节4个系数。 ushort iircas51 (DATA *x, DATA *h, DATA *r, DATA *dbuffer,ushort nbiq, ushort nx): IIR级联直 接1型，每节5个系数。 ushort iirlat (DATA *x, DATA *h, DATA *r, DATA *pbuffer,int nx, int nh):格型逆IIR滤波器。
7.6.5 DSPLIB函数简介 ushort firlat (DATA *x, DATA *h, DATA *r, DATA *pbuffer,int nx, int nh):格型前向FIR滤波器。 ushort firdec (DATA *x, DATA *h, DATA *r, DATA *dbuffer,ushort nh, ushort nx, ushort D):抽样FIR滤波器。 ushort firinterp (DATA *x, DATA *h, DATA *r, DATA *dbuffer,ushort nh, ushort nx, ushort I):插值FIR滤波器。 ushort hilb16 (DATA *x, DATA *h, DATA *r, DATA *dbuffer,ushort nx, ushort nh):FIR Hilbert变换器。 ushort iir32 (DATA *x, LDATA *h, DATA *r, LDATA *dbuffer,ushort nbiq, ushort nr):双精度IIR。 自适应滤波 ushort dlms (DATA *x, DATA *h, DATA *r, DATA *des,DATA *dbuffer, DATA step, ushort nh, ushort nx)： LMS 滤波器 (延迟版)
7.6.5 DSPLIB函数简介 相关 ushort acorr (DATA *x, DATA *r, ushort nx, ushort nr, type)：自相关 (只保留正边)。 ushort corr (DATA *x, DATA *y, DATA *r, ushort nx, ushort ny, type)：相关 (双边)。 数学 ushort add (DATA *x, DATA *y, DATA *r, ushort nx,ushort scale)：优化后的向量加法。 ushort expn (DATA *x, DATA *r, ushort nx) :指数。 short bexp (DATA *x, ushort nx) :指数。 ushort logn (DATA *x, LDATA *r, ushort nx) :自然对数。 ushort log_2 (DATA *x, LDATA *r, ushort nx):基2对数。 ushort log_10 (DATA *x, LDATA *r, ushort nx) :基10对数 short maxidx (DATA *x, ushort ng, ushort ng_size)：求向量中最大值的序号。 short maxidx34 (DATA *x, ushort nx) :求向量中最大值的序号(nx≤34)。
7.6.5 DSPLIB函数简介 short maxval (DATA *x, ushort nx) :求向量中最大值。 void maxvec (DATA *x, ushort nx, DATA *r_val, DATA *r_idx)：求向量中最大值及其序号。 short minidx (DATA *x, ushort nx) :求向量中最小值的序号。 short minval (DATA *x, ushort nx) :求向量中的最小值。 void minvec (DATA *x, ushort nx, DATA *r_val,DATA *r_idx):求向量中最小值及其序号。 ushort mul32 (LDATA *x, LDATA *y, LDATA *r, ushort nx) :32位向量乘法。 short neg (DATA *x, DATA *r, ushort nx)：16位取负。 short neg32 (LDATA *x, LDATA *r, ushort nx) :32位取负。 short power (DATA *x, LDATA *r, ushort nx) :向量的平方和(功率)。 void recip16 (DATA *x, DATA *r, DATA *rexp, ushort nx) :返回Q15格式向量x各元素的尾数和指数。 void ldiv16 (LDATA *x, DATA *y, DATA *r, DATA *rexp,ushort nx) :32位除以16位的除法。 ushort sqrt_16 (DATA *x, DATA *r, short nx): 向量平方根； short sub (DATA *x, DATA *y, DATA *r, ushort nx,ushort scale)：向量减法。
7.6.5 DSPLIB函数简介 三角 ushort sine (DATA *x, DATA *r, ushort nx)：向量的正弦。 ushort atan2_16 (DATA *i, DATA *q, DATA *r, short nx)：四象限反正切。 ushort atan16 (DATA *x, DATA *r, ushort nx) ：反正切。 矩阵 ushort mmul (DATA *x1, short row1, short col1,DATA *x2, short row2, short col2, DATA *r): 矩阵乘。 ushort mtrans (DATA *x, short row, short col, DATA *r) ): 矩阵转置 其他 ushort fltoq15 (float *x, DATA *r, ushort nx) :浮点转Q15。 ushort q15tofl (DATA *x, float *r, ushort nx) :Q15转浮点。 ushort rand16 (DATA *r, ushort nr):随机数据产生。 void rand16init(void):初始化随机数据。
7.6.6 DSPLIB函数的调用 从C程序中调用一个DSPLIB函数的方法如下： 在C源程序中包含dsplib.h。 在CCS中则修改链接器选项，把你的代码与 55xdsp.lib或55xdspx.lib链接。 例7-14，DSPLIB函数的调用。通过调用 q15tofl(x, r, NX)函数，将向量x[8]和标量y[1]由 Q15小数格式转化为浮点格式。 源程序(设文件名为Ex7_14.c)如下：
7.6.6 DSPLIB函数的调用 #include <dsplib.h> #define NX 8 #define NY 1 DATA x[8]={-17621,7002,-919,25644,17176,-2853,-1556,21063}; DATA y[1]={17621}; float r[NX]; float s[NY]; void main(void) { q15tofl(x,r,NX); //把Q15小数格式的Nx=8维向量x转化为浮点格式向量r q15tofl(y,s,NY); //把Q15小数格式的Ny=1维向量x转化为浮点格式向量s return; }
7.6.6 DSPLIB函数的调用 在CCS54中调试该程序的步骤如下： ①建立工程, 设工程名为Ex7_14; ②建立C源程序Ex7_14.c; ③打开工程属性对话框, 修改链接器选项,见图7-13;
7.6.6 DSPLIB函数的调用 ④单击工具按钮 对工程进行构建，并进入调试环境; ⑤打开表达式显示窗口，添加变量x、r、y、s； ⑥在 “q15tofl(y, s, NY);”语句行号处双击鼠标左键, 设置断点； ⑦单击工具按钮 运行程序至断点处,可看到图7-14所示运行结果。
第7章作业 7.7 (仿真调试)利用第3章中C语言工程所用命 令文件Exam3_2.cmd ，在CCS下建立关于例 7-8的语言工程，并对工程进行构建、调试。 7.1, 7.2, 7.6 (参考习题看书)