基于TMS320VC5416的简易信号分析仪设计

1、 DSP原理与应用教程课程设计 DSP原理与应用教程课程设计说明书课题：基于TMS320VC5416的简易信号分析仪设计专业： 电子信息工程 班级： 电子信息工程1102班 姓名： 陈玮 学号： 3110209424 指导老师： 方卫东、李锦彬、朱悦涵 2014年 12 月 14日 - 12月 26日 摘要系统基于快速傅立叶变换(FFT)算法，以TMS320VC5416为控制与数据处理核心，通过CCS开发环境以及MATLAB模拟实现对音频信号频率成分的分析。系统由电源模块、音频输入输出模块、时钟模块等组成。通过A/DC采样并进行FFT变换，分析音频信号的频谱结构。关键词： 音频信号,CCS，F

2、FT, 频谱目录 TOC o 1-3 h z u HYPERLINK l _Toc407306637 摘要 PAGEREF _Toc407306637 h 2 HYPERLINK l _Toc407306638 一、设计任务与目标 PAGEREF _Toc407306638 h 4 HYPERLINK l _Toc407306639 1.1设计课题 PAGEREF _Toc407306639 h 4 HYPERLINK l _Toc407306640 1.2设计要求 PAGEREF _Toc407306640 h 4 HYPERLINK l _Toc407306641 二、谐波的概念 PAGE

3、REF _Toc407306641 h 4 HYPERLINK l _Toc407306642 三、谐波的分析方法 PAGEREF _Toc407306642 h 4 HYPERLINK l _Toc407306643 2.1模拟电路 PAGEREF _Toc407306643 h 4 HYPERLINK l _Toc407306644 2.2傅立叶变换 PAGEREF _Toc407306644 h 5 HYPERLINK l _Toc407306645 2.3小波变换 PAGEREF _Toc407306645 h 5 HYPERLINK l _Toc407306646 四、DFT/FFT

4、分析原理 PAGEREF _Toc407306646 h 5 HYPERLINK l _Toc407306647 4.1 DFT计算公式 PAGEREF _Toc407306647 h 5 HYPERLINK l _Toc407306648 4.2 N点DFT的计算量 PAGEREF _Toc407306648 h 6 HYPERLINK l _Toc407306649 4.3旋转因子WN的特性 PAGEREF _Toc407306649 h 6 HYPERLINK l _Toc407306650 4.4基-2 FFT算法推导 PAGEREF _Toc407306650 h 6 HYPERLI

5、NK l _Toc407306651 4.5 N点基-2 FFT算法的计算量 PAGEREF _Toc407306651 h 10 HYPERLINK l _Toc407306652 4.6 N点基-2 FFT算法的实现方法 PAGEREF _Toc407306652 h 10 HYPERLINK l _Toc407306653 4.6.1对于输入数据序列进行倒位序变换 PAGEREF _Toc407306653 h 10 HYPERLINK l _Toc407306654 4.6.2蝶形运算的循环结构 PAGEREF _Toc407306654 h 11 HYPERLINK l _Toc40

6、7306655 4.6.3浮点到定点转换需要注意的关键问题 PAGEREF _Toc407306655 h 11 HYPERLINK l _Toc407306656 4.6.4计算过程中的溢出问题 PAGEREF _Toc407306656 h 12 HYPERLINK l _Toc407306657 五、硬件系统 PAGEREF _Toc407306657 h 12 HYPERLINK l _Toc407306658 5.1硬件电路的总体框图 PAGEREF _Toc407306658 h 12 HYPERLINK l _Toc407306659 5.2各子模块的组成 PAGEREF _To

7、c407306659 h 13 HYPERLINK l _Toc407306660 5.2.1电源模块 PAGEREF _Toc407306660 h 13 HYPERLINK l _Toc407306661 5.2.2时钟模块 PAGEREF _Toc407306661 h 14 HYPERLINK l _Toc407306662 5.2.3音频输入输出模块 PAGEREF _Toc407306662 h 15 HYPERLINK l _Toc407306663 5.2.4 TMS320VC5416和CODEC接口 PAGEREF _Toc407306663 h 16 HYPERLINK l

8、 _Toc407306664 5.2.5 TMS320VC5416DSP最小系统 PAGEREF _Toc407306664 h 17 HYPERLINK l _Toc407306665 六、基于MATLAB的程序仿真 PAGEREF _Toc407306665 h 18 HYPERLINK l _Toc407306666 6.1程序流程图 PAGEREF _Toc407306666 h 18 HYPERLINK l _Toc407306667 6.2 MATLAB构造音频信号 PAGEREF _Toc407306667 h 19 HYPERLINK l _Toc407306668 6.3仿真

9、步骤 PAGEREF _Toc407306668 h 19 HYPERLINK l _Toc407306669 6.4仿真结果 PAGEREF _Toc407306669 h 20 HYPERLINK l _Toc407306670 6.5结果分析 PAGEREF _Toc407306670 h 20 HYPERLINK l _Toc407306671 七、课程设计总结 PAGEREF _Toc407306671 h 21 HYPERLINK l _Toc407306672 八、参考文献 PAGEREF _Toc407306672 h 22 HYPERLINK l _Toc407306673

10、附录 PAGEREF _Toc407306673 h 22一、设计任务与目标1.1设计课题基于TMS320VC5416的简易信号分析仪设计1.2设计要求1.给出算法原理2.写出应用软件流程图3.以TMS320VC5416 DSP为核心，设计一DSP应用系统，用DSP C语言和汇编混合编程的方法设计应用软件，实现音频信号功率谱分析，并验证最终的分析结果。最后按要求写出设计说明书，绘出设计的软硬件系统.1.3设计目标 学习FFT算法在功率谱分析中的应用；以TMS320VC5416 DSP为核心CPU，设计一小型DSP应用系统，编制软件实现频谱分析，在CCS上观察分析的结果.二、谐波的概念任何周期性

11、波形均可分解为一个基频正弦波加上许多高次频率的正弦波，高次频率是基频的整倍数（N，只能为整数），直流成分称为0次谐波，基波称为1次谐波，二次以上的波形称为高次谐波，其中偶次频率的波形称为偶次谐波，奇次频率的波形称为奇次谐波。三、谐波的分析方法2.1.模拟电路 消除谐波的方法很多，即有主动型，又有被动型；既有无源的，也有有源的，还有混合型的，目前较为先进的是采用有源电力滤波器。但由于其检测环节多采用模拟电路，因而造价较高，且由于模拟带通滤波器对频率和温度的变化非常敏感，故使其基波幅值误差很难控制在10%以内，严重影响了有源滤波器的控制性能。近年来，人工神经 HYPERLINK /network/

12、 网络的研究取得了较大进展，由于神经元有自适应和自学习能力，且结构简单，输入输出关系明了，因此可用神经元替代自适应滤波器，再用一对与基波频率相同，相位相差90度的正弦向量作为神经元的输入。由神经元先得到基波电流，然后检测出应补偿的电流，从而完成谐波电流的检测。但人工神经网络的硬件目前还是一个比较薄弱的环节，限制了其应用范围。2.2.傅立叶变换利用傅立叶变换可在数字域进行谐波检测，电力系统的谐波分析，目前大都是通过该方法实现的，离散傅立叶变换所需要处理的是经过采样和A/D转换得到的数字信号，设待测信号为x(t)，采样间隔为 t秒，采样频率 =1/ t满足采样定理，即 大于信号最高频率分量的2倍，

13、则采样信号为x(n t)，并且采样信号总是有限长度的，即n=0，1N-1。这相当于对无限长的信号做了截断，因而造成了傅立叶变换的泄露现象，产生误差。此外，对于离散傅立叶变换来说，如果不是整数周期采样，那么即使信号只含有单一频率，离散傅立叶变换也不可能求出信号的准确参数，因而出现栅栏效应。通过加窗可以减小泄露现象的影响。2.3.小波变换小波变换已广泛应用于信号分析、语音识别与合成、自动控制、图象处理与分析等领域。电力谐波是由各种频率成分合成的、随机的、出现和消失都非常突然的信号，在应用离散傅立叶变换进行处理受到局限的情况下，可充分发挥小波变换的优势。即对谐波采样离散后，利用小波变换对数字信号进行

14、处理，从而实现对谐波的精确测定。小波可以看作是一个双窗函数，对一信号进行小波变换相当于从这一时频窗内的信息提取信号。对于检测高频信息，时窗变窄，可对信号的高频分量做细致的观测；对于分析低频信息，这时时窗自动变宽，可对信号的低频分量做概貌分析。所以小波变换具有自动“调焦”性。其次，小波变换是按频带而不是按频点的方式处理频域信息，因此信号频率的微小波动不会对处理产生很大的影响，并不要求对信号进行整周期采样。另外，由小波变换的时间局部可知，在信号的局部发生波动时，不会象傅立叶变换那样把影响扩散到整个频谱，而只改变当时一小段时间的频谱分布，因此，采用小波变换可以跟踪时变和暂态信号。四、DFT/FFT分

15、析原理4.1 DFT计算公式其中x(n)表示输入的离散数字信号序列，WN为旋转因子，X(k)为输入序列x(n)对应的N个离散频率点的相对幅度。一般情况下，假设x(n)来自于低通采样，采样频率为fs，那么X(k)表示了从-fs/2率开始，频率间隔为fs/N,到fs/2-fs/N截至的N个频率点的相对幅度。因为DFT计算得到的一组离散频率幅度值实际上是在频率轴上从成周期变化的，即X(k+N)=X(k)。因此任意取连续的N个点均可以表示DFT的计算效果，负频率成分比较抽象，难于理解，根据X(k)的周期特性，于是我们又可以认为X(k)表示了从零频率开始，频率间隔为fs/N,到fs-fs/N截至的N个频

16、率点的相对幅度。4.2 N点DFT的计算量根据(1)式给出的DFT计算公式，我们可以知道每计算一个频率点X(k)均需要进行N次复数乘法和N-1次复数加法，计算N各点的X(k)共需要N2次复数乘法和N*(N-1)次复数加法。当x(n)为实数的情况下，计算N点的DFT需要2*N2次实数乘法，2*N*(N-1)次实数加法。4.3旋转因子WN的特性1.WN的对称性 2.WN的周期性3.WN的可约性根据以上这些性质，我们可以得到式(5)的一系列有用结果4.4基-2 FFT算法推导假设采样序列点数为N=2L，L为整数，如果不满足这个条件可以人为地添加若干个0以使采样序列点数满足这一要求。首先我们将序列x(

17、n)按照奇偶分为两组如下：于是根据DFT计算公式(1)有：至此，我们将一个N点的DFT转化为了式(7)的形式，此时k的取值为0到N-1,现在分为两段来讨论，当k为0N/2-1的时候，因为x1(r)，x2(r)为N/2点的序列，因此，此时式(7)可以写为：而当 k取值为N/2N-1时，k用k+N/2取代，k取值为0N/2-1。对式(7)化简可得：综合以上推导我们可以得到如下结论：一个N点的DFT变换过程可以用两个N/2点的DFT变换过程来表示，其具体公式如式(10)所示DFT快速算法的迭代公式：上式中X(k)为偶数项分支的离散傅立叶变换，X(k)为奇数项分支的离散傅立叶变换。 式(10)的计算过

18、程可以用图1的蝶形算法流图直观地表示出来。图1 时间抽取法蝶形运算流图在图1中，输入为两个N/2点的DFT输出为一个N点的DFT结果，输入输出点数一致。运用这种表示方法，8点的DFT可以用图2来表示：图2 8点DFT的4点分解根据公式(10)，一个N点的DFT可以由两个N/2点的DFT运算构成，再结合图1的蝶形信号流图可以得到图2的8点DFT的第一次分解。该分解可以用以下几个步骤来描述： 1.将N点的输入序列按奇偶分为2组分别为N/2点的序列 2.分别对1中的每组序列进行DFT变换得到两组点数为N/2的DFT变换值X1和X2 3.按照蝶形信号流图将2的结果组合为一个N点的DFT变换结果根据式(

19、10)我们可以对图2中的4点DFT进一步分解，得到图3的结果，分解步骤和前面一致。图3 8点DFT的2点分解最后对2点DFT进一步分解得到最终的8点FFT信号计算流图：图4 8点DFT的全分解从图2到图4的过程中关于旋转系数的变化规律需要说明一下。看起来似乎向前推一级，在奇数分组部分的旋转系数因子增量似乎就要变大，其实不是这样。事实上奇数分组部分的旋转因子指数每次增量固定为1，只是因为每向前推进一次，该分组序列的数据个数变少了，为了统一使用以原数据N为基的旋转因子就进行了变换导致的。每一次分组奇数部分的系数WN,这里的N均为本次分组前的序列点数。以上边的8点DFT为例，第一次分组N=8,第二次

20、分组N为4，为了统一根据式(4)进行了变换将N变为了8，但指数相应的需要乘以2。4.5 N点基-2 FFT算法的计算量从图4可以看到N点DFT的FFT变换可以转为log2(N)级级联的蝶形运算，每一级均包含有N/2次蝶形计算。而每一个蝶形运算包含了1次复数乘法，2次复数加法。因此N点FFT计算的总计算量为：复数乘法N/2log2(N) 复数加法Nlog2(N)。假设被采样的序列为实数序列，那么也只有第一级的计算为实数与复数的混合计算，经过一次迭代后来的计算均变为复数计算，在这一点上和直接的DFT计算不一致。因此对于输入序列是复数还是实数对FFT算法的效率影响较小。一次复数乘法包含了4次实数乘法

21、，2次实数加法，一次复数加法包含了2次复数加法。因此对于N点的FFT计算需要总共的实数乘法数量为：2Nlog2(N)；总的复数加法次数为：2xNxlog2(N)。4.6 N点基-2 FFT算法的实现方法从图4我们可以总结出对于点数为N=2L的DFT快速计算方法的流程：4.6.1对于输入数据序列进行倒位序变换该变换的目的是使输出能够得到X(0)X(N-1)的顺序序列，同样以8点DFT为例，该变换将顺序输入序列x(0)x(7)变为如图4的x(0),x(4),x(2),x(6),x(1),x(5),x(3),x(7)序列。其实现方法是：假设顺序输入序列一次村在A(0)A(N-1)的数组元素中，首先我

22、们将数组下标进行二进制化(例：对于点数为8的序列只需要LOG2(8) = 3位二进制序列表示，序号6就表示为110)。二进制化以后就是将二进制序列进行倒位，倒位的过程就是将原序列从右到左书写一次构成新的序列，例如序号为6的二进制表示为110，倒位后变为了011，即使十进制的3。第三步就是将倒位前和倒位后的序号对应的数据互换。依然以序号6为例，其互换过程如下：Temp = A(6); A(6) = A(3); A(3) = A(6);实际上考虑到执行效率，如果对于每一次输入的数据都需要这个处理过程是非常浪费时间的。我们可以采用指向指针的指针来实现该过程，或者是采用指针数组来实现该过程。4.6.2

23、蝶形运算的循环结构从图4中我们可以看到对于点数为N = 2L的fft运算，可以分解为L阶蝶形图级联，每一阶蝶形图内又分为M个蝶形组，每个蝶形组内包含K个蝶形。根据这一点我们就可以构造三阶循环来实现蝶形运算。编程过程需要注意旋转因子与蝶形阶数和蝶形分组内的蝶形个数存在关联。4.6.3浮点到定点转换需要注意的关键问题上边的分析都是基于浮点运算来得到的结论，事实上大多数嵌入式系统对浮点运算支持甚微，因此在嵌入式系统中进行离散傅里叶变换一般都应该采用定点方式。对于简单的DFT运算从浮点到定点显得非常容易。根据式(1)，假设输入x(n)是经过AD采样的数字序列，AD位数为12位，则输入信号范围为0409

24、6。为了进行定点运算我们将旋转因子实部虚部同时扩大212倍，取整数部分代表旋转因子。之后，我们可以按照(1)式计算，得到的结果与原结果成比例关系，新的结果比原结果的212倍。但是，对于使用蝶形运算的fft我们不能采用这种简单的放大旋转因子转为整数计算的方式。因为fft是一个非对称迭代过程，假设我们对旋转因子进行了放大，根据蝶形流图我们可以发现其最终的结果是，不同的输入被放大了不同的倍数，对于第一个输入x(0)永远也不会放大。举一个更加形象的例子，还是以图4为例。从图中可以看出右侧的X(0)可以直接用下式表示：从上式我们可以看到不同输入项所乘的旋转因子个数(注意这里是个数，就算是wn0，也被考虑

25、进去了，因为在没有放大时wn0等于1，放大后所有旋转因子指数模均不为1，因此需要考虑)。这就导致输入不平衡，运算结果不正确。经查阅相关资料，比较妥善的做法是，首先对所有旋转因子都放大2Q倍，Q必须要大于等于L，以保证不同旋转因子的差异化。旋转因子放大，为了保证其模为1，在每一次蝶形运算的乘积运算中我们需要将结果右移Q位来抵消这个放大，从而得到正确的结果。之所以采用放大倍数必须是2的整数次幂的原因也在于此，我们之后可以通过简单的右移位运算将之前的放大抵消，而右移位又代替了除法运算，大大节省了时间。4.6.4计算过程中的溢出问题最后需要注意的一个问题就是计算过程中的溢出问题。在实际应用中，AD虽然

26、有12位的位宽，但是采样得到的信号可能较小，例如可能在08之间波动，也就是说实际可能只有3位的情况。这种情况下为了在计算过程中不丢失信息，一般都需要先将输入数据左移P位进行放大处理，数据放大可能会导致溢出，从而使计算错误，而溢出的极限情况是这样：假设我们数据位宽为D位(不包括符号位)，AD采样位数B位，数字放大倍数P位，旋转因此放大倍数Q位，FFT级联运算带来的最大累加倍数L位。我们得到:假设AD位宽12，数据位宽32，符号位1位，因此有效位宽31位，采样点数N，那么我们可以得到log2(N)+P+Q=L可以得到放大倍数P VECT PAGE 0.sysregs: BIOSREGS PAGE

27、1 .trcinit: EPROG PAGE 0 .gblinit: EPROG PAGE 0 frt: EPROG PAGE 0 .text: EPROG PAGE 0/* 已初始化段，包含所有的可执行代码、字符串和编译器生成的常数 */ .cinit: EPROG PAGE 0/* 已初始化段，包含初始化变量和常数表 */ .pinit: EPROG PAGE 0/* 已初始化段，包含运行时调用的全局目标构造器列表 */ .sysinit: EPROG PAGE 0 .bss: IDATA PAGE 1/* 未初始化段 */ .far: IDATA PAGE 1 .const: IDATA

28、 PAGE 1/* 已初始化段，包含由C限定词const定义的字符串常量和数据 */ .switch: IDATA PAGE 1/* 已初始化段 */ .sysmem: IDATA PAGE 1/* 未初始化段，为动态存储器函数malloc等分配存储器空间 */ .cio: IDATA PAGE 1 .MEM$obj: IDATA PAGE 1 .sysheap: IDATA PAGE 1 .stack: IDATA PAGE 1/* 未初始化段，为系统堆栈分配存储器 */ fft.c#include #include #include #include #include fft.h#incl

29、ude sin.incconst float dataR256=/y=100*sin(2*pi/N)*2*n)+200*sin(2*pi/N)*4*n)+300*sin(2*pi/N)*6*n)0,68.529,135.91,201,262.72,320.05,372.05,417.89,456.85,488.37,511.99,527.43,534.57,533.42,524.17,507.16,482.84,451.84,414.88,372.78,326.45,276.88,225.08,172.11,119,66.802,16.49,-31.003,-74.826,-114.22,-1

30、48.55,-177.27,-200,-216.48,-226.58,-230.34,-227.9,-219.56,-205.74,-186.96,-163.84,-137.1,-107.51,-75.891,-43.102,-10,22.563,53.768,82.843,109.08,131.86,150.65,165.02,174.66,179.4,179.16,174.01,164.13,149.82,131.49,109.64,84.883,57.869,29.322,7.3479e-014,-29.322,-57.869,-84.883,-109.64,-131.49,-149.8

31、2,-164.13,-174.01,-179.16,-179.4,-174.66,-165.02,-150.65,-131.86,-109.08,-82.843,-53.768,-22.563,10,43.102,75.891,107.51,137.1,163.84,186.96,205.74,219.56,227.9,230.34,226.58,216.48,200,177.27,148.55,114.22,74.826,31.003,-16.49,-66.802,-119,-172.11,-225.08,-276.88,-326.45,-372.78,-414.88,-451.84,-48

32、2.84,-507.16,-524.17,-533.42,-534.57,-527.43,-511.99,-488.37,-456.85,-417.89,-372.05,-320.05,-262.72,-201,-135.91,-68.529,1.6111e-012,68.529,135.91,201,262.72,320.05,372.05,417.89,456.85,488.37,511.99,527.43,534.57,533.42,524.17,507.16,482.84,451.84,414.88,372.78,326.45,276.88,225.08,172.11,119,66.8

33、02,16.49,-31.003,-74.826,-114.22,-148.55,-177.27,-200,-216.48,-226.58,-230.34,-227.9,-219.56,-205.74,-186.96,-163.84,-137.1,-107.51,-75.891,-43.102,-10,22.563,53.768,82.843,109.08,131.86,150.65,165.02,174.66,179.4,179.16,174.01,164.13,149.82,131.49,109.64,84.883,57.869,29.322,-3.1247e-013,-29.322,-5

34、7.869,-84.883,-109.64,-131.49,-149.82,-164.13,-174.01,-179.16,-179.4,-174.66,-165.02,-150.65,-131.86,-109.08,-82.843,-53.768,-22.563,10,43.102,75.891,107.51,137.1,163.84,186.96,205.74,219.56,227.9,230.34,226.58,216.48,200,177.27,148.55,114.22,74.826,31.003,-16.49,-66.802,-119,-172.11,-225.08,-276.88

35、,-326.45,-372.78,-414.88,-451.84,-482.84,-507.16,-524.17,-533.42,-534.57,-527.43,-511.99,-488.37,-456.85,-417.89,-372.05,-320.05,-262.72,-201,-135.91,-68.529;float re_inputK_FFT_SIZE,im_inputK_FFT_SIZE;/ FFT/DFT运算输入float am_outputK_FFT_SIZE;/ FFT/DFT输出谐波幅值float re_outputK_FFT_SIZE,im_outputK_FFT_SIZ

36、E;/ DFT运算输出谐波的实部和虚部void clr_ram(void) unsigned int i;for (i=0;iK_FFT_SIZE;i+) am_outputi=0.0;re_inputi=im_inputi=0.0;re_outputi=im_outputi=0.0;/ 用于DFTvoid put_in(unsigned int N) int i;for(i=0;iN;i+) re_inputi=dataR256i/1.0;/ 实部im_inputi=0.0;/ 虚部/ 位倒序运算void re_bit(unsigned int N)int i,j,M,s;float tem

37、p;M=N/2;i=0;for(j=1;j=s) i=i-s;s=s/2;i=i+s;if(ji) temp=re_inputj;re_inputj=re_inputi;re_inputi=temp;temp=im_inputj;im_inputj=im_inputi;im_inputi=temp;/ 蝶形运算void butterfly(unsigned int l) int i,j,r;/分别是级号、蝶群号、运算蝶号int m1,m2,m3,m4;int N,k1,k2;float u,v;int x;N=1l;for(i=1;il+1;i+)/第i级/m1=pow(2,i-1);/m1=

38、1,2,4,8,16,32,64,128,.m1=1(i-1);m2=2*m1;/第i级的蝶群长度m3=N/m2;/第i级包含的蝶形单元，即“群”个数for(j=0;jm3;j+)/第j号蝶群m4=j*m2;/m4=j*(2i) 蝶群址for(r=0;rm1;r+)/第r号运算蝶k1 = (m4+r);/s=j*(2i)+r 蝶址k2 = (k1+m1);/u=Rek2*cos(2*p*r/m2)+Imk2*sin(2*p*r/m2);/v=Imk2*cos(2*p*r/m2)-Rek2*sin(2*p*r/m2);x = (1(l-i);x = x*r;u = (re_inputk2*cos

39、_tab256x)+(im_inputk2*sin_tab256x);v= (im_inputk2*cos_tab256x)-(re_inputk2*sin_tab256x);re_inputk2 = (re_inputk1-u);im_inputk2 = (im_inputk1-v);re_inputk1 = (re_inputk1+u);im_inputk1 = (im_inputk1+v);/ FFT运算/ input: n=采样点数N（必须是2的l次幂）void fft(unsigned int n,unsigned int l) re_bit(n);/ 位倒序butterfly(l)

40、;/ DFT运算，计算结果为X(0),X(1),.X(N-1)/ sign=1:正变换;sign=-1:反变换void dft(int sign) int i,k;double c,q,w,s;q=(2.0*PI)/K_FFT_SIZE;for (k=0; kK_FFT_SIZE; k+) w=k*q;re_outputk=im_outputk=0.0;for (i=0; iK_FFT_SIZE; i+) unsigned int l;/d=i*w;/c=cos(d);/s=sin(d)*sign;l=(k*i)%K_FFT_SIZE;c=cos_tab256l;s=sin_tab256l;r

41、e_outputk += (c*re_inputi+s*im_inputi);/ 计算X(k)的实部im_outputk += (c*im_inputi-s*re_inputi);/ 计算X(k)的虚部if (sign=-1) c=1.0/K_FFT_SIZE;for (k=0; kK_FFT_SIZE; k+) re_outputk=c*re_outputk;im_outputk=c*im_outputk;/ 计算单次X(k)，即只计算某一次谐波的值/ input:k=谐波次数void dft_k(unsigned int k) int i;double c,s;re_outputk=im_

42、outputk=0.0;for (i=0; iK_FFT_SIZE; i+) unsigned int l;l=(k*i)%K_FFT_SIZE;c=cos_tab256l;s=sin_tab256l;re_outputk += (c*re_inputi+s*im_inputi);/ 计算X(k)的实部 im_outputk += (c*im_inputi-s*re_inputi);/ 计算X(k)的虚部/ 求FFT运算结果谐波的幅值/ input:n=采样点数=谐波个数/ output:am_output=谐波的幅值void get_am(unsigned int N) unsigned i

43、nt i;for (i=0; iN; i+)am_outputi=sqrt(re_inputi*re_inputi+im_inputi*im_inputi);am_outputi /= (float)K_FFT_SIZE;am_outputi *= 2;if (am_outputi ZERO_VALVE)/ 屏蔽计算带来的0的误差am_outputi = 0.0;void get_am_dft(void) unsigned int i;for (i=0; iK_FFT_SIZE; i+)am_outputi=sqrt(re_outputi*re_outputi+im_outputi*im_ou

44、tputi);am_outputi /= (float)K_FFT_SIZE;am_outputi *= 2;if (am_outputi ZERO_VALVE)am_outputi = 0.0;#defineEXE_FFT0/ FFT运算#defineEXE_DFT1/ DFT计算所有谐波#defineEXE_DFT_K2/ DFT计算指定次谐波#defineEXE_MODEEXE_FFTvoid main() int m;clr_ram();/ 清内存put_in(K_FFT_SIZE);/ 伪A/D采样#if EXE_MODE=EXE_FFT/ 快速傅里叶变换FFTfft(K_FFT_S

45、IZE, K_LOGN);/ 得X(k),k=0,1,.,N-1get_am(K_FFT_SIZE);/ 求谐波幅值/ 快速傅里叶逆变换IFFTfor (m=0; mK_FFT_SIZE; m+)/ 求X(k)共轭复数im_inputm=-im_inputm;fft(K_FFT_SIZE, K_LOGN);for (m=0; mK_FFT_SIZE; m+)/ 最后取共轭并乘以1/N得到序列x(n),n=0,1,.,N-1im_inputm=-im_inputm/K_FFT_SIZE;re_inputm/=K_FFT_SIZE;if (abs(re_inputm) ZERO_VALVE)re_

46、inputm = 0.0;if (abs(im_inputm) ZERO_VALVE)im_inputm = 0.0;asm( NOP);#endif#if EXE_MODE=EXE_DFTdft(1);get_am_dft();asm( NOP);#endif#if EXE_MODE=EXE_DFT_Kdft_k(5);get_am_dft();asm( NOP);#endifasm( NOP);while(1);基于MATLAB音频构造代码function st=mstgN=256Fs=2800;T=1/Fs;Tp=N*T;t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;fc1=

47、Fs/256*2;fc2=Fs/256*4;fc3=Fs/256*6;xt1=100*sin(2*pi*fc1*t);xt2=200*sin(2*pi*fc2*t);xt3=300*sin(2*pi*fc3*t);st=xt1+xt2+xt3;fxt=fft(st,N);dlmwrite(dsp.txt,st);subplot(3,1,1)plot(t,st);grid;xlabel(t/s);ylabel(s(t);axis(0,Tp,min(st),max(st);title(a)st的波形);subplot(3,1,2)stem(f,abs(fxt)/max(abs(fxt),.);gr

48、id;title(b)st的频谱);axis(0,Fs/5,0,1.2);xlabel(f/Hz);ylabel(幅度);基于MATLAB的SIN/COS的查表法代码fidc=fopen(sin.inc,wt);fprintf(fidc,%s,const float cos_tab256=);fprintf(fidc,%s,1.0000,);for n=1:255a(n)=cos(2*pi/255)*1*n);fprintf(fidc,%.4f,a(n);endfprintf(fidc,%s,);fprintf(fidc,nnn);fprintf(fidc,%s,const float sin

49、_tab256=);fprintf(fidc,%s,0,);for n=1:255a(n)=sin(2*pi/255)*1*n);fprintf(fidc,%.4f,a(n);endfprintf(fidc,%s,);fclose(fidc);sin.inc: const floatcos_tab256=1.0000,0.9997,0.9988,0.9973,0.9951,0.9924,0.9891,0.9852,0.9806,0.9755,0.9698,0.9635,0.9566,0.9491,0.9411,0.9325,0.9233,0.9135,0.9032,0.8924,0.8810,

50、0.8691,0.8566,0.8437,0.8302,0.8162,0.8017,0.7867,0.7713,0.7554,0.7390,0.7222,0.7049,0.6872,0.6691,0.6506,0.6317,0.6124,0.5928,0.5727,0.5524,0.5317,0.5106,0.4893,0.4677,0.4457,0.4235,0.4011,0.3784,0.3555,0.3324,0.3090,0.2855,0.2618,0.2379,0.2139,0.1898,0.1656,0.1412,0.1168,0.0923,0.0677,0.0431,0.0185

51、,-0.0062,-0.0308,-0.0554,-0.0800,-0.1045,-0.1290,-0.1534,-0.1777,-0.2019,-0.2260,-0.2499,-0.2737,-0.2973,-0.3207,-0.3439,-0.3670,-0.3898,-0.4124,-0.4347,-0.4567,-0.4785,-0.5000,-0.5212,-0.5421,-0.5626,-0.5828,-0.6026,-0.6221,-0.6412,-0.6599,-0.6782,-0.6961,-0.7136,-0.7307,-0.7473,-0.7634,-0.7791,-0.

52、7943,-0.8090,-0.8233,-0.8370,-0.8502,-0.8629,-0.8751,-0.8868,-0.8979,-0.9085,-0.9185,-0.9280,-0.9369,-0.9452,-0.9529,-0.9601,-0.9667,-0.9727,-0.9781,-0.9830,-0.9872,-0.9908,-0.9939,-0.9963,-0.9981,-0.9993,-0.9999,-0.9999,-0.9993,-0.9981,-0.9963,-0.9939,-0.9908,-0.9872,-0.9830,-0.9781,-0.9727,-0.9667

53、,-0.9601,-0.9529,-0.9452,-0.9369,-0.9280,-0.9185,-0.9085,-0.8979,-0.8868,-0.8751,-0.8629,-0.8502,-0.8370,-0.8233,-0.8090,-0.7943,-0.7791,-0.7634,-0.7473,-0.7307,-0.7136,-0.6961,-0.6782,-0.6599,-0.6412,-0.6221,-0.6026,-0.5828,-0.5626,-0.5421,-0.5212,-0.5000,-0.4785,-0.4567,-0.4347,-0.4124,-0.3898,-0.

54、3670,-0.3439,-0.3207,-0.2973,-0.2737,-0.2499,-0.2260,-0.2019,-0.1777,-0.1534,-0.1290,-0.1045,-0.0800,-0.0554,-0.0308,-0.0062,0.0185,0.0431,0.0677,0.0923,0.1168,0.1412,0.1656,0.1898,0.2139,0.2379,0.2618,0.2855,0.3090,0.3324,0.3555,0.3784,0.4011,0.4235,0.4457,0.4677,0.4893,0.5106,0.5317,0.5524,0.5727,

55、0.5928,0.6124,0.6317,0.6506,0.6691,0.6872,0.7049,0.7222,0.7390,0.7554,0.7713,0.7867,0.8017,0.8162,0.8302,0.8437,0.8566,0.8691,0.8810,0.8924,0.9032,0.9135,0.9233,0.9325,0.9411,0.9491,0.9566,0.9635,0.9698,0.9755,0.9806,0.9852,0.9891,0.9924,0.9951,0.9973,0.9988,0.9997,1.0000,const float sin_tab256=0,0.0246,0.0493,0.0739,0.0984,0.1229,0.1473,0.1716,0.1958,0.2199,0.2439,0.