DSP课程设计报告.docx
《DSP课程设计报告.docx》由会员分享，可在线阅读，更多相关《DSP课程设计报告.docx（18页珍藏版）》请在冰豆网上搜索。
DSP课程设计报告
JIANGSUUNIVERSITY
DSP原理及应用课程设计报告
——FFT的DSP实现
姓名：
专业班级：
通信1002
学号：
31006010
指导老师：
宋雪桦
设计日期：
2013.05.06~2013.05.15
一、设计目的
1、加深对DFT算法原理和基本性质的理解；
2、了解并学习使用FFT算法，以及其在TMS320C54X上的运用；
3、学习DSP中FFT的设计和编程思想；
4、练习使用CCS的探针和图形工具来观察器观察波形和频谱情况。
二、设计内容
用C语言及汇编语言进行编程，实现FFT运算，对于C语言，实现8点和16点的FFT运算，对于汇编语言，需调试出8点的FFT运算结果。
三、设计原理
快速傅里叶变换（FFT）是一种高效实现离散傅里叶变换（DFT）的快速算法，是数字信号处理中最为重要的工具之一，它在声学，语音，电信和信号处理等领域有着广泛的应用。
1、离散傅里叶变换DFT
对于长度为N的有限长序列x（n），它的离散傅里叶变换（DFT）为
（1）
式中，
，称为旋转因子或蝶形因子。
从DFT的定义可以看出，在x（n）为复数序列的情况下，对某个k值，直接按
（1）式计算X（k）只需要N次复数乘法和（N-1）次复数加法。
因此，对所有N个k值，共需要N2次复数乘法和N（N-1）次复数加法。
对于一些相当大有N值（如1024点）来说，直接计算它的DFT所需要的计算量是很大的，因此DFT运算的应用受到了很大的限制。
2、快速傅里叶变换FFT
旋转因子WN有如下的特性:
对称性：
周期性：
利用这些特性，既可以使DFT中有些项合并，减少了乘法积项，又可以将长序列的DFT分解成几个短序列的DFT。
FFT就是利用了旋转因子的对称性和周期性来减少运算量的。
FFT的算法是将长序列的DFT分解成短序列的DFT。
例如：
N为偶数时，先将N点的DFT分解为两个N/2点的DFT，使复数乘法减少一半：
再将每个N/2点的DFT分解成N/4点的DFT，使复数乘又减少一半，继续进行分解可以大大减少计算量。
最小变换的点数称为基数，对于基数为2的FFT算法，它的最小变换是2点DFT。
一般而言，FFT算法分为按时间抽取的FFT（DITFFT）和按频率抽取的FFT（DIFFFT）两大类。
DIFFFT算法是在时域内将每一级输入序列依次按奇/偶分成2个短序列进行计算，而DIFFFT算法是在频域内将每一级输入序列依次奇/偶分成2个短序列进行计算。
两者的区别是旋转因子出现的位置不同，但算法是一样的。
在DIFFFT算法中，旋转因子
出现在输入端，而在DIFFFT算法中它出现在输入端。
假定序列x（n）的点数N是2的幂，按照DIFFFT算法可将其分为偶序列和奇序列，记偶序列为
，奇序列为
，则x（n）的FFT表示为
由于
，则（3）式可表示为
式中，
和
分别为
和
的N/2的DFT。
由于对称性，
则
。
因此，N点
可分为两部分：
前半部分：
（4）
后半部分：
（5）
从式（4）和式（5）可以看出，只要求出0~N/2-1区间
和
的值，就可求出0~N-1区间
的N点值。
以同样的方式进行抽取，可以求得N/4点的DFT，重复抽取过程，就可以使N点的DFT用上组2点的DFT来计算，这样就可以大减少运算量。
在基数为2的FFT中，设N=2M，共有M级运算，每级有N/2个2点FFT蝶形运算，因此，N点FFT总共有MN/2个蝶形运算。
蝶形运算如图1所示。
图1蝶形运算
设蝶形输入为
和
，输出为
和
，则有
（6）
（7）
在基数为2的FFT中，设N=2M，共有M级运算，每级有N/2个2点FFT蝶形运算，因此，N点FFT总共有
个蝶形运算。
例如：
基数为2的FFT，当N=8时，共需要3级，12个基2DITFFT的蝶形运算。
其信号流程如图2所示。
图28点基2DIFFFT蝶形运算
从图2可以看出，输入是经过比特反转的倒位序列，称为位码倒置，其排列顺序为
。
输出是按自然顺序排列，其顺序为
。
3、FFT运算的实现
（1）实现输入数据的比特反转
输入数据的比特反转实际上就是将输入数据进行码位倒置，以便在整个运算后输出序列是一个自然序列。
在用汇编指令进行码位倒置时，使用位码倒置寻址可以大大提高程序执行速度和使用存储器的效率。
在这种寻址方式下，AR0存放的整数N的FFT点的一半，一个辅助寄存器指向一个数据存放单元。
当使用位码倒置寻址将AR0加到辅助寄存器时，地址将以位码倒置的方式产生。
（2）实现N点复数FFT
N点复数FFT算法的实现可分为三个功能块，即第一级蝶形运算、第二级蝶形运算、第三级至log2N级蝶形运算。
队以任何一个2的整数幂N=2^M，总可以通过M次分解后最后成为二点的DFT运算。
通过这样的M次分解，可以构成M级迭代计算，每级由N/2个蝶形运算完成。
（4）输出FFT结果
FFT算法的程序流程图如下图所示。
图3FFT运算的程序流程图
四、基于C语言的FFT算法
1、源程序FFT.c
#include"myapp.h"
#include"ICETEK-VC5509-EDU.h"
#include"scancode.h"
#include
#definePI3.1415926
#defineSAMPLENUMBER16
voidMakeWave（）;
intINPUT[SAMPLENUMBER];
structcompx{floatreal,imag;};
structcompxEE（structcompx,structcompx）;
structcompxxin[SAMPLENUMBER];
structcompxfWave[SAMPLENUMBER];
floaty[SAMPLENUMBER],data[SAMPLENUMBER];
structcompxEE（structcompxb1,structcompxb2）
{
structcompxb3;
b3.real=b1.real*b2.real-b1.imag*b2.imag;b3.imag=b1.real*b2.imag+b1.imag*b2.real;
return（b3）;
}
main（）
{
inti;
MakeWave（）;
for（i=0;i{
fWave[i].real=INPUT[i];
fWave[i].imag=0.0f;
y[i]=0.0f;
}
FFT（fWave）;
for（i=0;i{
data[i]=y[i];
}
while
（1）;//breakpoint
}
voidFFT（structcompx*xin）
{
intf,m,nv2,nm1,i,k,j=0,l;
structcompxv,w,t;
nv2=SAMPLENUMBER/2;
f=SAMPLENUMBER;
for（m=1;（f=f/2）!
=1;m++）{;}
nm1=SAMPLENUMBER-1;
for（i=0;i{
if（i{
t=xin[j];
xin[j]=xin[i];
xin[i]=t;
}
k=nv2;
while（k<=j）
{
j=j-k;
k=k/2;
}
j=j+k;
}
{intle,lei,ip;
for（l=1;l<=m;l++）
{
le=2<<（l-1）;
lei=le/2;
v.real=1.0;v.imag=0.0;
w.real=cos（PI/lei）;
w.imag=-sin（PI/lei）;
for（j=0;j<=lei-1;j++）
{
for（i=j;i<=SAMPLENUMBER;i=i+le）
{
ip=i+lei;
t=EE（xin[ip],v）;
xin[ip].real=xin[i].real-t.real;
xin[ip].imag=xin[i].imag-t.imag;
xin[i].real=xin[i].real+t.real;
xin[i].imag=xin[i].imag+t.imag;
}
v=EE（v,w）;
}
}
}
for（i=0;i{
y[i]=sqrt（xin[i].real*xin[i].real+xin[i].imag*xin[i].imag）;
}
}
voidMakeWave（）
{
inti;
for（i=0;i{
INPUT[i]=i+1;
}
}
2、链接命令文件ICETEK-VC5509-A.cmd
-w
-stack500
-sysstack500
-lrts55x.lib
MEMORY
{
DARAM:
o=0x100,l=0x7f00
VECT:
o=0x8000,l=0x100
DARAM2:
o=0x8100,l=0x7f00
SARAM:
o=0x10000,l=0x30000
SDRAM:
o=0x40000,l=0x3e0000
}
SECTIONS
{
.text:
{}>DARAM
.vectors:
{}>VECT
.trcinit:
{}>DARAM
.gblinit:
{}>DARAM
frt:
{}>DARAM
.cinit:
{}>DARAM
.pinit:
{}>DARAM
.sysinit:
{}>DARAM
.bss:
{}>DARAM2
.far:
{}>DARAM2
.const:
{}>DARAM2
.switch:
{}>DARAM2
.sysmem:
{}>DARAM2
.cio:
{}>DARAM2
.MEM$obj:
{}>DARAM2
.sysheap:
{}>DARAM2
.sysstack{}>DARAM2
.stack:
{}>DARAM2
}
3、FFT程序的使用方法
（1）根据N值，修改FFT.c中的中的常数，如N=8，将#defineSAMPLENUMBER16语句中的“16”修改为8。
（2）编译、汇编、链接，得到.out文件，加载。
（3）将data加入观察窗，可看到FFT运算输出结果。
4、运行结果
8点的FFT运算，且输入为1、2、3、…、8时，运算结果如图4所示，16点的FFT运算，且输入为1、2、3、…、16时，运算结果如图5所示。
图48点FFT运算结果
图516点FFT运算结果
五、基于汇编语言的FFT算法
1、汇编源程序fft.asm
.title"fft.asm"
.mmregs
.include"coeff.inc"
.include"in.inc"
.defstart
sine:
.usect"sine",512
sine1:
.usect"sine1",512
cosine:
.usect"cosine",512
cosine1:
.usect"cosine1",512
fft_data:
.usect"fft_data",1024
fft_out:
.usect"fft_out",512
STACK.usect"STACK",10
K_DATA_IDX_1.set2
K_DATA_IDX_2.set4
K_DATA_IDX_3.set8
K_TWID_TBL_SIZE.set512
K_TWID_IDX_3.set128
K_FLY_COUNT_3.set4
K_FFT_SIZE.set8
K_LOGN.set3
PA0.set0
.bssd_twid_idx,1
.bssd_data_idx,1
.bssd_grps_cnt,1
.sect"fft_prg"
********************位码倒置程序**************************
.asgAR2,REORDERED
.asgAR3,ORIGINAL_INPUT
.asgAR7,DATA_PROC_BUF
start:
SSBXFRCT
STM#STACK+10,SP
STM#sine,AR1
RPT#K_TWID_TBL_SIZE-1
MVPD#sine1,*AR1+
STMcosine,AR1
RPT#K_TWID_TBL_SIZE-1
MVPD#cosine1,*AR1+
STM#d_input,ORIGINAL_INPUT
STM#fft_data,DATA_PROC_BUF
MVMMDATA_PROC_BUF,REORDERED
STM#K_FFT_SIZE-1,BRC
RPTBDbit_rev_end-1
STM#K_FFT_SIZE,AR0
MVDD*ORIGINAL_INPUT+,*REORDERED+
MVDD*ORIGINAL_INPUT-,*REORDERED+
MAR*ORIGINAL_INPUT+0B
bit_rev_end:
*********************FFTCODE*********************************
.asgAR1,GROUP_COUNTER
.asgAR2,PX
.asgAR3,QX
.asgAR4,WR
.asgAR5,WI
.asgAR6,BUTTERFLY_COUNTER
.asgAR7,STAGE_COUNTER
*******************第一级蝶形运算stage1************************
STM#0,BK
LD#-1,ASM
STM#fft_data,PX
STM#fft_data+K_DATA_IDX_1,QX
STMK_FFT_SIZE/2-1,BRC
LD*PX,16,A
RPTBDstage1end-1
STM#K_DATA_IDX_1+1,AR0
SUB*QX,16,A,B
ADD*QX,16,A
STHA,ASM,*PX+
STB,*QX+
||LD*PX,A
SUB*QX,16,A,B
ADD*QX,16,A
STHA,ASM,*PX+0%
STB,*QX+0%2
||LD*PX,A
stage1end:
******************第二级蝶形运算stage2***************************
STM#fft_data,PX
STM#fft_data+K_DATA_IDX_2,QX
STM#K_FFT_SIZE/4-1,BRC
LD*PX,16,A
RPTBDstage2end-1
STM#K_DATA_IDX_2+1,AR0
;1stbutterfly
SUB*QX,16,A,B
ADD*QX,16,A
STHA,ASM,*PX+
STB,*QX+
||LD*PX,A
SUB*QX,16,A,B
ADD*QX,16,A
STHA,ASM,*PX+
STHB,ASM,*QX+
;2ndbutterfly
MAR*QX+
ADD*PX,*QX,A
SUB*PX,*QX-,B
STHA,ASM,*PX+
SUB*PX,*QX,A
STB,*QX
||LD*QX+,B
STA,*PX
||ADD*PX+0%,A
STA,*QX+0%
||LD*PX,A
stage2end:
********************第三级至最后一级蝶形运算************************
STM#K_TWID_TBL_SIZE,BK
ST#K_TWID_IDX_3,d_twid_idx
STM#K_TWID_IDX_3,AR0
STM#cosine,WR
STM#sine,WI
STM#K_LOGN-2-1,STAGE_COUNTER
ST#K_FFT_SIZE/8-1,d_grps_cnt
STM#K_FLY_COUNT_3-1,BUTTERFLY_COUNTER
ST#K_DATA_IDX_3,d_data_idx
stage:
STM#fft_data,PX
LDd_data_idx,A
ADD*（PX）,A
STLMA,QX
MVDKd_grps_cnt,GROUP_COUNTER
group:
MVMDBUTTERFLY_COUNTER,BRC
RPTBDbutterflyend-1
LD*WR,T
MPY*QX+,A
MAC*WI+0%,*QX-,A
ADD*PX,16,A,B
STB,*PX
||SUB*PX+,B
STB,*QX
||MPY*QX+,A
MAS*QX,*WR+0%,A
ADD*PX,16,A,B
STB,*QX+
||SUB*PX,B
LD*WR,T
STB,*PX+
||MPY*QX+,A
butterflyend:
PSHMAR0
MVDKd_data_idx,AR0
MAR*PX+0
MAR*QX+0
BANZDgroup,*GROUP_COUNTER-
POPMAR0
MAR*QX-
LDd_data_idx,A
SUB#1,A,B
STLMB,BUTTERFLY_COUNTER
STLA,1,d_data_idx
LDd_grps_cnt,A
STLA,ASM,d_grps_cnt
LDd_twid_idx,A
STLA,ASM,d_twid_idx
BANZDstage,*STAGE_COUNTER-
MVDKd_twid_idx,AR0
fft_end:
***********************计算功率谱
STM#fft_data,AR2
;STM#fft_data,AR3
STM#fft_out,AR4
STM#K_FFT_SIZE*2-1,BRC
RPTBpower_end-1
SQUR*AR2+,A
SQURA*AR2+,A
STHA,*AR4+
power_end:
STM#fft_out,AR4
RPT#K_FFT_SIZE-1
PORTW*AR4+,PA0
NOP
NOP
here:
Bhere
.end
2、链接命令文件fft.cmd
vector.obj
rfft_task.obj
-orfft_task.obj
-mrfft_task.map
-erfft_task
MEMORY
{
PAGE0:
EPROM:
org=0E000Hlen=1000H
VECS:
org=0FF80Hlen=0080H
PAGE1:
SPRAM:
org=0060Hlen=0020H
DARAM:
org=0400Hlen=0600H
RAM:
org=8000Hlen=1400H
}
SECTIONS
{
sine1:
>EPROMPAGE0
cosine1:
>EPROMPAGE0
fft_prg:
>EPROMPAGE0
.bss:
>SPRAMPAGE1
sine:
align（512）{}>DARAMPAGE1
cosine:
align（512）{}>DARAMPAGE1
d_input:
>RAMPAGE1
fft_data:
>RAMPAGE1
fft_out:
>RAMPAGE1
STACK:
>SPRAMPAGE1
.vectors:
>VECSPAGE0
}
3、FFT程序的使用方法
（1）根据N值，修改rfft_task.asm中的两个常数，如N=64.
K_FFT_SIZE.set64
K_LOGN.set6
（2）准备输入数据文件in.dat。
输入数据按实部、虚部，实部、虚部，……顺序存放。
（3）汇编、链接、仿真执行，得到输出数据文件out.dat。
（4）根据out.dat作图，就可以得到输入信号的功率谱图。
六、设计心得体会
升级会员 