信号的频谱分析及DSP实现.doc

1、现代信号处理课程设计课程设计报告设计题目 信号的频谱分析及 DSP 实现 目 录第 1 章 绪论11.1 设计背景.11.2 设计要求.31.3 设计思路简介.3第 2 章 系统开发平台与环境.32.1 CCS 开发环境.32.2 SEED-DEC2812 开发实验箱3第 3 章 频谱分析仪设计过程 .43.1 FFT 算法实现 .4 3.2 频谱分析系统设计 8 3.3 源程序代码 9 第 4 章 系统仿真14 4.1 仿真设置 .14 4.2 仿真图 16第 6 章 总结20参考文献.201第 1 章 绪论1.1 设计背景随着数字技术与计算机技术的发展，数字信号处理(DSP)技术已深入到各

2、个学科领域。近些年来，数字信号处理技术同数字计算器、大规模集成电路等，有了突飞猛进的发展。在数字信号处理中，离散傅里叶变换(DiscreteTime Fourier Transform，DFT)是常用的变换方法，它在数字信号处理系统中扮演着重要角色。由离散傅里叶变换发现频率离散化，可以直接用来分析信号的频谱、计数滤波器的频率响应，以及实现信号通过线系统的卷积运算等，因而在信号的频谱分析方面有很大的作用。由于 DFT 的运算量太大，即使是采用计算机也很难对问题进行实时处理，所以经过很多学者的不懈努力，便出现了通用的快速傅里叶变换（FFT） 。快速傅里叶变换（Fast Fourier Transf

3、orm，FFT）并不是与离散傅里叶变换不同的另一种变换，而是为了减少 DFT 计算次数的一种快速有效的算法。对 FFT 算法及其实现方式的研究是很有意义的。目前，FFT 己广泛应用在频谱分析、匹配滤波、数字通信、图像处理、语音识别、雷达处理、遥感遥测、地质勘探和无线保密通讯等众多领域。在不同应用场合，需要不同性能要求的 FFT 处理器。在很多应用领域都要求 FFT 处理器具有高速度、高精度、大容量和实时处理的性能。因此，如何更快速、更灵活地实现 FFT 变得越来越重要。数字信号处理器（DSP）是一种可编程的高性能处理器。它不仅是一种适用于数字信号处理，而且在图像处理、语音处理、通信等领域得到广

4、泛的应用。DSP处理器中集成有高速的乘法硬件，能快速的进行大量的乘法加法运算。频谱分析在生产实践和科学研究中获得日益广泛的应用。例如，对汽车、飞机、轮船、汽轮机等各类旋转机械、电机、机床等机器的主体或部件进行实际运行状态下的谱分析，可以提供设计数据和检验设计效果，或者寻找振源和诊断故障，保证设备的安全运行等；在声纳系统中，为了寻找海洋水面船只或2潜艇，需要对噪声信号进行谱分析，以提供有用信息，判断舰艇运动速度、方向、位置、大小等。因此对谱分析方法的研究，受到普遍注意和重视，是当前信号处理技术中一个十分活跃的课题。1965 年库利首次提出了快速傅里叶变换(FFT)算法，FFT 和频谱分析很快发展

5、成为机械设备故障诊断、振动分析、无线电通信、信息图象处理和自动控制等多种学科重要的理论基础。然而长期的应用和近年来的理论分析表明：经快速傅立叶变换得到的离散频谱，频率、幅值和相位均可能产生较大误差，单谐波加矩形窗时最大误差从理论上分析可达 36.4%；即使加其他窗时，也不能完全消除此影响，在加汉宁（Hanning）窗时，只进行幅值恢复时的最大幅值误差仍高达 15.3%，相位误差高达 90 度。因此，频谱分析的结果在许多领域只能定性而不能精确的定量分析和解决问题，大大限制了该技术的工程应用，特别是在机械振动和故障诊断中的应用受到极大限制。从 70 年代中期，有关学者开始致力于频谱校正理论的研究以

6、期解决离散频谱误差较大的问题。1975 年 John CBurges 等从事电学领域研究工作的学者采用插值法对加矩形窗的离散化频谱进行校正，解决了电学中的离散高次谐波参数的精确测量问题。1983 年 Thomas Grandke 提出了加 Hanning 窗的内插法，进一步提高了离散高次谐波参数的分析精度。1993 年，丁康和谢明提出了三点卷积法幅值校正法，提高了频率间隔较大的信号的离散频谱幅值精度，解决了工程实际中的一些问题。1994 年，谢明、丁康等提出和发展了比例频谱校正方法，使内差法系统地发展成为一种通用的频谱校正方法，解决了频率间隔较大的离散化频谱幅值、相位和频率的精确求解问题，并开

7、始对离散频谱的校正方法和误差分析进行了深入系统的分析和研究。1996 年，余佳兵，史铁林等提出了采用复调制细化谱分析将已产生频谱干涉的密集频率成分分离开，消除干涉，再用比例法进行校正以解决密集频率成分的离散频谱的校正问题。1997 年，谢明、丁康等分析了离散频谱中的负频率成分和多频谱成分的干涉现象，提出了离散频谱中用相位和幅值综合判定和识别单频率成分的方法，实现了单频率成分和频率间隔较大的多频率成分的自动识别和自动校正，并提出了在不采长样的基础上利用轴系旋转识别和校正两个己发生干涉的密集频率成分的自动判定和校正的方法。1998 年刘渝提出了一段信号作 N 点和 N/2 点的校正方法，利用相位3

8、信息可以得到比较精确的频率。1999 年，丁康、谢明等提出了对连续时域信号分前后两段作傅里叶变换，利用其对应离散谱线的相位差校正出谱峰处的准确频率和相位的校正方法一相位差校正法，该方法可在不知道窗谱函数表达式的情况下，直接用其相位差进行频率和相位校正。2001 年，徐培民、杨积东、闻邦椿提出了自动识别和修正离散频谱中两临近谱峰参数的方法，不仅能识别间距不到一个频率分辨率的两个密集频率成分，而且能识别峰间距为 1-6 个频率的临近谱峰参数。1.2 设计要求利用 C 语言在 CCS 环境中编写一个 FFT 程序，并能利用已设计好的频谱分析仪对常用信号进行频谱分析。1.3 设计思路简介在 TMS32

9、0C54x 系统开发环境 CCS（Code Composer Studio）下对 FFT 算法的 DSP 实现原理进行讨论。利用 C 语言设计相应的频谱分析仪，通过实验仿真，从输入信号和输出信号的时域和频域曲线可看出在 DSP 上实现的信号的频谱分析能完成预定的任务。第 2 章 系统开发平台与环境2.1 CCS 开发环境CCS 提供了配置、建立、调试、跟踪和分析程序的工具，它便于实时、嵌入式信号处理程序的编制和测试，它能够加速开发进程，提高工作效率。CCS 提供了基本的代码生成工具，它们具有一系列的调试、分析能力。CCS支持如下图 1 所示的开发周期的所有阶段。4图 12.2 SEED-DEC

10、2812 开发实验箱SEED-DECxxxx 系列嵌入式 DSP 开发板本着模块化、总线型、开放式、系列化的设计思想，采用统一的系统结构、模块结构和机械结构，以多种典型 DSP处理器构成具有标准总线和相同物理尺寸的高性能嵌入式 DSP 开发板。SEED-DEC2812 嵌入式 DSP 开发板原理框图如图 2 所示：图 2第 3 章 频谱分析仪设计过程3.1 FFT算法实现快速傅里叶变换（FFT）是离散傅里叶变换的快速算法，它是根据离散傅里叶变换的奇、偶、虚、实等特性，对离散傅里叶变换的算法进行改进获得的，它对离散傅里叶变换并没有新的发现。有限长序列x(n)及其频域表示X(k)可由以下离散傅立叶

11、变换得出设 计概 念 性 规 划 编 程 和 编 译创 建 工 程 文 件 、 编 写源 代 码 和 配 置 文 件 调 试语 法 检 查 、 探测 点 设 置 和 日志 保 存 等 分 析实 时 调 试 、 统计 和 跟 踪510x(k)=DFT(n)()NnkxW01N（8）10()IX(k)()Nnkn（9）其中 。式（8）称为离散傅立叶正变换，式（9）称为离散傅立叶逆2jnknkNWe变换，x(n)与X(k)构成了离散傅立叶变换对。根据上述公式，计算一个X(k)，需要N次复数乘法和N-1次复数加法，而计算全部X(k)( )，共需要 次复数乘法和N(N-1)次复数加法。实现01k2N一次

12、复数乘法需要四次实数乘法和两次实数加法，一次复数加法需要两次实数加法，因此直接计算全部X(k)共需要4 次实数乘法和2N(2N-1)次实数加法。2当N较大时，对实时信号处理来说，对处理器计算速度有十分苛刻的要求，于是如何减少计算离散傅里叶变换运算量的问题变得至关重要。为减少运算量，提高运算速度，就必须改进算法。计算DFT过程中需要完成的运算的系数里，存在相当多的对称性。通过研究这种对称性，可以简化计算过程中的运算，从而减少计算DFT所需的时间。如前所述，N点的DFT的复乘次数等于 。显然，把N点的DFT分解为几个较2短的DFT，可是乘法的次数大大减少。另外，旋转因子 具有明显的周期性和mNW对

13、称性，其周期为： 22()jmlNjmmlNNWee其对称性表现为：或N*mNWFFT算法就是不断地把长序列的DFT分解成几个短序列的DFT，并利用 的mNW周期性和对称性来减少DFT的运算次数。具有以下固有特性：nkNW(1) 的周期性： ()(NnkNknkW）6(2) 的对称性： nkNW()nknknNkNW（ ）(3) 的可约性： /,nn另外， 。/2(/2)1,NkNk利用 的上述特性，将x(n)或X(k)序列按一定规律分解成短序列进行运nk算，这样可以避免大量的重复运算，提高计算DFT的运算速度。算法形式有很多种，但基本上可以分为两大类，即按时间抽取(Decimation In

14、 Time，DIT)FFT算法和按频率抽取(Decimation In Frequency，DIF)FFT算法。如果序列x(n)的长度 ，其中M是整数(如果不满足此条件，可以人为地增2N补零值点来达到)，在时域上按奇偶抽取分解成短序列的DFT，使最小DFT运算单元为2点。通常将FFT运算中最小DFT运算单元称为基(radix)，因而把这种算法称为基-2时间抽取FFT(DIT-FFT)算法 4。将x(n)按n为奇偶分解成两个子序列，当n为偶数时，令n=2r；当n为奇数时，令n=2r+l；可得到12(2),()(,0.,12Nxrxrr（10）则其 DFT 可写成 1122(21)00()()(N

15、NrkrkNrXkxWx1122(21)00()rkrkNr1122/ /200()()rkrkNrxx12kNXW（11）和 均分别是N/2点序列 和 的DFT，而且r与k的取值满1()Xk2()1()xn2足0,1,N/2-1。而X(k)是一个N点的DFT，因此式（11）只计算了X(k)的前N/2的值。由DFT和 的性质可得到X(k)的后N/2的值为：nk7（12）21()()()2NkNXkWX12()kNWX式（11）和式（12）表明，只要计算出两个 N/2 点的 DFT 和 ，1()xk2经过线性组合，即可求出全部 N 点的 X(k)。由于 ， 仍为偶数，2M/M因而这样的分解可以继

16、续进行下去，直到最后的单元只需要做 2 点 DFT 为止。若 Xm(p)和 Xm(q)为输入数据， 和 为输出数据， 为旋转1()mXp1()mqnkNW因子，则对于基-2DIT-FFT 算法，蝶形运算的基本公式为1()()()kmmmNXpXpXqWq其图形表示如图3所示，称Xm(P)为上结点，Xm(q)为下结点。图3 时间抽取蝶形运算单元对于一个8点的FFT，根据上述算法可以得到一个完整的N=8的基-2DIT-FFT的运算流图，如图4所示。图4 N=8 DIT-FFT运算流图8根据上述算法原理及运算流图，可以得出基-2DIT-FFT的基本特点，特点如下。(1)级数分解：对于 。共分了M级，

17、每级包含N/2个蝶形运算单元，总共2N所需蝶形运算个数为 。2logN(2)运算量估计：每个蝶形运算需要一次复数乘法和两次复数加(减)法，N点FFT共需要 次复数乘法， 次复数加(减)法。实际上有些蝶形运2log2l算不需要做复乘。(3)原位运算：当数据输入到存储器以后，每一组蝶形运算后，结果仍然存放在这同一组存储器中的同一位置，不需要另辟存储单元，直到最后输出。(4)位码倒序：由图2可以看到，FFT输出的X(k)的次序正好是顺序排列的，即X(0)，X(1)，,X(7)，而输入X(n)是按x(0)，X(4)，,X(7)的倒序存入存储单元，即为倒序输入，正序输出。这种顺序看起来相当杂乱，然而它是

18、有规律的，即位码倒序规则。(5)旋转因子的确定：由8点FFT的三次迭代运算可以看出 的变化。在第一级kNW迭代中，只有一种类型的蝶形运算系数，即 ，参加蝶形运算的两个数据点间08隔为l；在第二级迭代中，有两种类型的蝶形运算系数，分别是 和 ，参加082蝶形运算的两个数据点间隔为2；在第三级迭代中，有四种类型的蝶形运算系数，分别是 ， ， ， ，参加蝶形运算的两个数据点间隔4。可见，每次迭08W1283代的蝶形类型比前一迭代增加一倍，间隔也增大一倍。最后一次迭代的蝶形类型最多，参加蝶形运算的两个数据点的间隔也最大，为N/2。3.2 频谱分析系统设计3.2.1硬件方案系统的硬件实现方案如图5所示：

19、9信号隔离信号滤波 AD 转换电源、复位时钟TMS320F2812USB 电平转换待检测信号上位机图5（1）信号前向处理电路设计信号的前向处理主要进行了隔离和低通滤波。隔离采用精密隔离放大器IS0124来实现，它通过封装在器件内部的2pF的隔离电容来实现输入级和输出级信号的隔离。它防止了可能损害信号的不同地电位之间的环路电流的干扰存在，阻止了可使信号受损的噪声通过信号的传播，将每个信号分离到一个干净的信号子系统。由IS0124输出的隔离后的信号带有大量的高频信号，所以在进行AD采样前要对信号进行模拟滤波。模拟滤波采用连续时间滤波器MAX274来完成。（2） AD转换设计F2812内部带有12位

20、的AD转换模块，考虑到信号采集和分析的精度，本系统外扩了14位的AD转换芯片MAXl320。其独立的采样保持(TH)电路为每一通道提供了同时采样，可以保持输入信号间的相对相位信息。3.2.2软件设计思路系统上电，进人main()函数后，首先进行系统的初始化，包括系统时钟初始化、外设中断向量表初始化、CPIO初始化、外部扩展总线XINTF初始化、SCI初始化、定时器初始化等，初始化完成后，启动定时器，使能中断，接着进入主循环，系统开始正常运转。开始初始化工作变量调用波形发生子程序产生波形(3 个正弦波)调用 FFT 子程序计算功率谱波形发生计算步长用标准 C 的 sin 函数计算当前波形值(12

21、8 点)结束FFT按照编码逆序排列输入序列返回计算结果用蝶形算法计算计算功率谱10图 63.3 源程序代码#include “DSP281x_Device.h“ / DSP281x Headerfile Include File#include “DSP281x_Examples.h“ / DSP281x Examples Include File#include “DSP281x_Device.h“ / DSP281x Headerfile Include File#include “DSP281x_Examples.h“ / DSP281x Examples Include File#in

22、clude “f2812a.h“#include“math.h“#define PI 3.1415926#define SAMPLENUMBER 128/ Prototype statements for functions found within this file.interrupt void adc_isr(void);void InitForFFT();void MakeWave();/void FFT(float dataRSAMPLENUMBER,float dataISAMPLENUMBER);/ Global variables used in this example:Ui

23、nt16 LoopCount;Uint16 ConversionCount;Uint16 Voltage11024;Uint16 Voltage21024;int i;int INPUTSAMPLENUMBER,DATASAMPLENUMBER;float fWaveRSAMPLENUMBER,fWaveISAMPLENUMBER,wSAMPLENUMBER;float sin_tabSAMPLENUMBER,cos_tabSAMPLENUMBER;main() InitSysCtrl();DINT;/InitPieCtrl(); IER = 0x0000;IFR = 0x0000;InitP

24、ieVectTable(); / Interrupts that are used in this example are re-mapped to/ ISR functions found within this file. EALLOW; / This is needed to write to EALLOW protected register11PieVectTable.ADCINT = EDIS; AdcRegs.ADCTRL1.bit.RESET = 1; / Reset the ADC moduleasm(“ RPT #10 | NOP“);AdcRegs.ADCTRL3.all

25、 = 0x00C8;AdcRegs.ADCTRL3.bit.ADCBGRFDN = 0x3; AdcRegs.ADCTRL3.bit.ADCPWDN = 1; / Power up rest of ADC/ Enable ADCINT in PIEPieCtrlRegs.PIEIER1.bit.INTx6 = 1;IER |= M_INT1; / Enable CPU Interrupt 1EINT; / Enable Global interrupt INTMERTM; / Enable Global realtime interrupt DBGMLoopCount = 0;Conversi

26、onCount = 0;/ Configure ADCAdcRegs.ADCMAXCONV.all = 0x0001; / Setup 2 convs on SEQ1AdcRegs.ADCCHSELSEQ1.bit.CONV00 = 0x0; / Setup ADCINA3 as 1st SEQ1 /conv.AdcRegs.ADCCHSELSEQ1.bit.CONV01 = 0x1; / Setup ADCINA2 as 2nd SEQ1 /conv.AdcRegs.ADCTRL2.bit.EVA_SOC_SEQ1 = 1; / Enable EVASOC to start /SEQ1Adc

27、Regs.ADCTRL2.bit.INT_ENA_SEQ1 = 1; / Enable SEQ1 interrupt /(every EOS)/ Configure EVA/ Assumes EVA Clock is already enabled in InitSysCtrl();EvaRegs.T1CMPR = 0x0080; / Setup T1 compare valueEvaRegs.T1PR = 0x10; / Setup period registerEvaRegs.GPTCONA.bit.T1TOADC = 1; / Enable EVASOC in EVAEvaRegs.T1

28、CON.all = 0x1042; / Enable timer 1 compare /(upcount mode)/ Wait for ADC interruptwhile(1)LoopCount+;12interrupt void adc_isr(void)Voltage1ConversionCount = AdcRegs.ADCRESULT0 4;Voltage2ConversionCount = AdcRegs.ADCRESULT1 4;/ If 40 conversions have been logged, start overif(ConversionCount = 1023)

29、ConversionCount = 0;else ConversionCount+;/ Reinitialize for next ADC sequenceAdcRegs.ADCTRL2.bit.RST_SEQ1 = 1; / Reset SEQ1AdcRegs.ADCST.bit.INT_SEQ1_CLR = 1; / Clear INT SEQ1 bitPieCtrlRegs.PIEACK.all = PIEACK_GROUP1; / Acknowledge interrupt /to PIEInitForFFT();MakeWave();for ( i=0;i0 ) b=b*2; i-;

30、 /* b= 2(L-1) */for ( j=0;j0 ) /* p=pow(2,7-L)*j; */p=p*2; i-;p=p*j;for ( k=j;kGraph-Time/Frequency进行如下图所示设置15图 7 观察窗口设置 1选择菜单 View-Graph-Time/Frequency进行如下图所示设置。图 8 观察窗口设置 2选择菜单 View-Graph-Time/Frequency进行如下图所示设置。16图 9 观察窗口设置 35清除显示：在以上打开的窗口中单击鼠标右键，选择弹出式菜单中“Clear Display”功能。6设置断点：在程序 FFT.c 中有注释“bre

31、ak point”的语句上设置软件断点。7.运行并观察结果选择“Debug”菜单的“Animate”项，或按 F12 键运行程序。观察“FFT”窗口中时域和频域图形。4.2 仿真图在现场环境中，通过仿真器与设备PCB板相连接，在CCS环境中可以查看采集到的波形以及经过FFT变换后的波形1718图 10 结果显示 119图 11 结果显示 220图 12 结果显示 3图 14 结果显示 421第 5 章 总结本论文学习和研究了快速傅里叶变换（FFT）的算法，把重点放在了时间抽取法基-2FFT 算法上。以及在 DSP 基础上用 FFT 变换对信号进行频谱分析。明确了 FFT 在 DSP 芯片上的实

32、现的关键。基于 DSP 的快速傅里叶变换频谱分析的研究使 FFT 能够有效的在 DSP 芯片上实现，有助于我们能够更及时的了解信息，对我们的生活生产以及科技研究有很大的帮助。自从快速傅里叶变换（FFT）出现以后，频谱分析技术便很快的发展起来，而且越来越贴近我们的生活生产，如医疗器械，无线电通信等等。但是我们对频谱分析技术的研究并未达到最高的层次，未来发展具有很广阔的空间。参考文献：1方勇.数字信号处理原理与实践M.北京：清华大学出版社.2006.2丁康，张晓飞.频谱校正理论的发展J.振动工程学报2000,13(1)：14-22. 3高西全,丁玉美. 数字信号处理M.西安：西安电子科技大学出版社

33、. 2008：75-86.4孙仲康快速傅里叶变换及其应用M. 北京：人民邮电出版社1982.5郑阿奇.MATLAB实用教程M.北京：电子工业出版社.2009：104-105,417-418.6郭仕剑，王宝顺，贺志国，杨可心.MATLAB7.X数字信号处理M。北京：中国邮电出版社.2006：24-26,105-108. 7赵桂芳等基于DSP的快速傅立叶变换的实现J黄石理工学院学报2007，23(5)：27-30.8乔瑞萍，崔涛，张芳娟.TMS32054xDSP原理及应用M.西安：西安电子科技大学出版社.2005：183-189.9杨宇，叶宇风，王洪.基于 DSP 的实时信号频谱分析模块设计J.电子测量技术.2006 年 4 月第二期：124-128.