基于dsp数字信号处理器的快速傅里叶变换程序设计

基于dsp数字信号处理器的快速傅里叶变换程序设计内容摘要：

/*开中断 */ IER |= M_INT14。 // 使能第 14 组中断（ CPU 定时器 2（ RTOS））。 //已经定义 define M_INT14 0x2020。 IER 是 CPU 中断使能寄存器 EINT。 InitAdc()。 for(。 ){ //等待 AD 采样结束。 if (adconvover==1){ //adconvover==1 表示 AD 采样 // 结束。 见 ad 中断子程序 for(i=0。 i128。 i++){ INPUT[i]=Ad_data[i]。 } for(i=0。 iN。 i++) //i128 { fWaveR[i]=INPUT[i]。 fWaveI[i]=。 w[i]=。 } Mum=(int)(+log(N)/log(2))。 //N=2^Mum FFT(fWaveR,fWaveI)。 for(i=0。 iN。 i++)DATA[i]=w[i]。 //功率谱 } adconvover=0。 } } 一般而言，编写一个能运行在操作系统上的程序，都需要一个主函数。 主函数意味着建立一个独立进程，且该进程成为了程序的入口，对其它各快速傅里叶变换程序设计 10 函数进行调用，当然其它被调用函数也可以再去调用更多函数，这样整个程序的运行轨迹就像一种栈，有时我们称之为调用栈。 主函数既是程序的入口，又是程序的出口。 由于主函数肩负着入口和出口的重任，所以最好不要把太多的细节方面的逻辑直接放在主函数内，这样不利于维护和扩展。 主函数应该尽量简洁，具体的实现细节应该封装到被调用的子函数里面去。 在这个主程序中，可以调用后续的 FFT 初始化函数程序，时间抽取法 FFT程序，以及倒序运算函数程序，是整个程序的中转站。 旋转因子的软件实现 for(m=1。 m=Mum。 m++) { B=(int)(pow(2,m1)+)。 //B=2^(m1) for(j=0。 jB。 j++) //每级需要进行 B种蝶 形运算 {S=j*(int)(pow(2,Mumm)+)。 for(k=j。 k=N1。 k+=(int)(pow(2,m)+)) { X=Xr[k+B]*cos_tab[S]+Xi[k+B]*sin_tab[S]。 //采用循环寻址方式对正弦 Y=Xi[k+B]*cos_tab[S]Xr[k+B]*sin_tab[S]。 //表和余弦表进行寻址 Xr[k+B]=Xr[k]X。 Xi[k+B]=Xi[k]Y。 Xr[k]=Xr[k]+X。 Xi[k]=Xi[k]+Y。 } } } 旋转因子是复数，可表示为： 2/ c os ( 2 / ) sin ( 2 / )k j k NNW e k N j k N    （ 41） 由式（ 41）可以看出旋转因子的实部为余弦函数，虚部为正弦函数。 为了获得FFT 运算中需要的全部旋转因子，需要分别存储正弦表和余弦表，且每个表长度为 N ，对应于 0176。 ～ 180176。 ，同时，采用循环寻址方式对正弦表和余弦表进行寻址。 FFT 初始化 void InitForFFT() //FFT 初始化函数，建立正余弦函数表 { int i。 for(i=0。 iN。 i++) 快速傅里叶变换程序设计 11 { sin_tab[i]=sin(PI*2*i/N)。 //建立正余弦函数表 cos_tab[i]=cos(PI*2*i/N)。 } } 为了获得 FFT 运算中需要的全部旋转因子，需要分别存储正弦表和余弦表，且每个表长度为 N ，对应于 0176。 ～ 180176。 按时间抽取法的 FFT 程序 void FFT(float Xr[N],float Xi[N]) //时间抽取法 FFT 程序，要求采样 //点数 N为 2的整数幂次方 { //Xr[]， Xi[]分别为输入序列的实 //部和虚部 int S,B。 //S 为旋转因子的幂数， B为蝶形运算输入数据的距离，也即各 //旋转因子的个数 int m,j,k。 float X,Y。 finv(N,Xr,Xi)。 //倒序运算函数，对输入序列倒序 for(m=1。 m=Mum。 m++) { B=(int)(pow(2,m1)+)。 //B=2^(m1) for(j=0。 jB。 j++) //每级需要进行 B种蝶形运算 {S=j*(int)(pow(2,Mumm)+)。 for(k=j。 k=N1。 k+=(int)(pow(2,m)+)) { //结果的实部和虚部分别存储在原实部和虚部位置 X=Xr[k+B]*cos_tab[S]+Xi[k+B]*sin_tab[S]。 Y=Xi[k+B]*cos_tab[S]Xr[k+B]*sin_tab[S]。 Xr[k+B]=Xr[k]X。 Xi[k+B]=Xi[k]Y。 Xr[k]=Xr[k]+X。 Xi[k]=Xi[k]+Y。 } } } 在这个时间抽取法 FFT 程序中，要求采样点数 N 为 2 的整数幂次方，每级需要进行 B种蝶形运算，每种蝶形运算在某一级中需要进行 N/pow(2,m)次蝶形运算，结果的实部和虚部分别存储在原实部和虚部位置。 快速傅里叶变换程序设计 12 对于任何一个 2的整数幂 N=2M，总是可以通过 M次分解最后完全成为 2点的DFT 运算。 这样的 M次分解，就构成从 x(n)到 X(k)的 M 级运算过程。 从上面的 流图可看到，每一级运算都由 N/2 个蝶形运算构成。 因此每一级运算都需要 次 复乘和 N 次复加 (每个结作加、减各一次 )，这样，经过时间抽取后 M级运算总共需要的运算： 复乘 复加 N 当然，实际情况与这个数字稍有出入，因为 这几个系数实际上都不用乘法运算，因此在上面 N=8 的例子中，实际上只有两个系数 W 及 W 是需要乘法运算的。 用时间抽取法所需的计算量，不论是复乘还是复加都与 Nlog2N成正比，而直接运算时则与 N2成正比。 例 N=2048， N2=4194304， (N/2)log2N=11264， N2/[(N/2)log2N]=。 FFT显然要比直接法快得多。 功率谱的计算实现 for(m=0。 mN/2。 m++) { w[m]=sqrt(Xr[m]*Xr[m]+Xi[m]*Xi[m])。 //计算功率谱 } 为了便于观察 FFT 的运算结果，需要求出信号的频谱。 经过第三级到最后一级蝶形运算之后，已经得到式（ 42）所示： ( ) ( ) ( )A k AR k AI k j （ 42） 故功率谱可以通过式（ 43）计算得到 ： 22( ) ( )AR k AI k （ 43） 通过软件仿真，可以观察到输入信号的时间波形和频谱波形以及输出信号的功率谱波形。 倒序运算函数 void finv(int N1,float *xr,float *xi) { int m,n,N2,k。 //m 为正序数； n为倒序数； k 为各个权值； N2为最高位的权值 快速傅里叶变换程序设计 13 float T。 //临时变量 T N2=N1/2。 //最高位加 1 相当于十进制加上最高位的权 N1/2 n=N2。 //第一个倒序值 for(m=1。 m=N12。 m++) //第 0 个和最后一个不倒序 { if(mn) //为了避免再次调换，只需对 mn 的部分调换顺序 { T=xr[m]。 xr[m]=xr[n]。 xr[n]=T。 T=xi[m]。 xi[m]=xi[n]。 xi[n]=T。 } k=N2。 //最高位权值 while(n=k) { n=nk。 //次高位位 1，继续上下进位，满 2 置 0 k=(int)(k/2+)。 //向下权值依次比上级减半 } n=n+k。 //得到下一倒序值 } } 倒序运算函数 finv（ N1,Xr,Xi），对输入序列倒序， N1 为序列长度， Xr[],Xi[]分别为输入序列的实部和虚部。 倒序原理：倒序数的加 1是在最高位加 1，满 2向次高位进 1，最高位变 0，依次往下。 从当前倒序值可求下一倒序值。 定时器 2 中断子程序及 AD转换子程序 interrupt void ad(void) { IFR=0x0000。 //CPU 级中断标志寄存器 IFR=0，即无 CPU 级中断请求 =0xffff。 //PIEACK 置 1，禁止外部中断向 CPU 发起中断 //请求 if(adconvover==0){ Ad_data[convcount] =。 // 将结果寄存器中的 AD //转 换结果放到 Ad_data 数组中 convcount++。 //convcount 的初始值为 0 } if (convcount==(128)) //即 convcount=128,即 Ad_data 包含采样值的最 //后一位为 Ad_data[128] { convcount=0。 adconvover=1。 //接满标志 ，即 AD采样结束 快速傅里叶变换程序设计 14 } } interrupt void ISRTimer2(void) { =1。 //清除 INT SEQ1 标志 //位 ,SEQ1 转换结束时 INT_SEQ1 中断标志位置 1 =1。 //启动 SEQ1，触发模数转换。 } 通过“ ConfigCpuTimer(amp。 CpuTimer2, 150, 22)。 ”设定定时器 2产生中断的时间。 当进入定时器 2 中断子程序后，启动 SEQ1，触发模数转换，进入 AD 中断子程 序。 快速傅里叶变换程序设计 15 5 工作过程分析 程序调试 在 CCS 下调试程序步骤 （ 1）点击 CCS 图标启动 CCS。 （ 2） 打开 工程。 将已编好的汇编源文件 及其他相关源文件和命令文件 等文件添加到工程中编译，链接产生 文件。 （ 3）下载 文件到目标板中。 （ 4）显示示输入信号的时域波形。 首先，通过设置实验箱的函数信号发生器，使 AD 采集波形为 0通道的标准正弦波。 （ 5） 显示输入信号的频域波形。 （ 6）显示信号功率谱。 由配置文件可知，经程序计算得到的信号功率谱放在 DATA 数组中。 显示经程序计算得到的信号功率谱。 （ 7）改变输入信号，如改为方波，三角波等等。 并重复（ 4）（ 5）（ 6）三个步骤。 输入信号时域波形 当输入信号为正弦波时，其时域波形如图 所示。 当输入信号为方波时，其时域波形如图 所示。 当输入信号为三角波时，其时域波形如图 所示。