当前位置：文档视界 › 基于DSP的FFT设计解读

基于DSP的FFT设计解读

自本1004班邓静 28号基于DSP的FFT滤波器设计

1 引言

随着数字技术与计算机技术的发展，数字信号处理(DSP)技术已深入到各个学科领域。近些年来，数字信号处理技术同数字计算器、大规模集成电路等，有了突飞猛进的发展。

在数字信号处理中，离散傅里叶变换(Discrete．Time Fourier Transform，DFT)是常用的变换方法，它在数字信号处理系统中扮演着重要角色。由离散傅里叶变换发现频率离散化，可以直接用来分析信号的频谱、计数滤波器的频率响应，以及实现信号通过线系统的卷积运算等，因而在信号的频谱分析方面有很大的作用。

由于DFT的运算量太大，即使是采用计算机也很难对问题进行实时处理，所以经过很多学者的不懈努力，便出现了通用的快速傅里叶变换（FFT）。快速傅里叶变换（Fast Fourier Transform，FFT）并不是与离散傅里叶变换不同的另一种变换，而是为了减少DFT计算次数的一种快速有效的算法。对FFT算法及其实现方式的研究是很有意义的。目前，FFT己广泛应用在频谱分析、匹配滤波、数字通信、图像处理、语音识别、雷达处理、遥感遥测、地质勘探和无线保密通讯等众多领域。在不同应用场合，需要不同性能要求的FFT处理器。在很多应用领域都要求FFT处理器具有高速度、高精度、大容量和实时处理的性能。因此，如何更快速、更灵活地实现FFT变得越来越重要。

数字信号处理器（DSP）是一种可编程的高性能处理器。它不仅是一种适用于数字信号处理，而且在图像处理、语音处理、通信等领域得到广泛的应用。DSP 处理器中集成有高速的乘法硬件，能快速的进行大量的乘法加法运算[1]。

本文主要介绍基于DSP用FFT变换实现对信号的频谱分析。研究离散傅里叶变换以及快速傅里叶变换的原理及算法。快速傅里叶变换和离散傅里叶变换的基本理论是一样的，它根据离散傅里叶变换的奇、偶、虚、实等特性，对离散傅里叶变换进行了改进。在计算机系统或者数字系统中广泛应用者快速傅里叶变换，这是一个巨大的进步。本文主要解决的问题就是如何对信号的频谱进行研究，使FFT更广泛的应用于科学研究。

2 快速傅里叶变换（FFT ）

2.1 FFT 算法基本原理

快速傅里叶变换（FFT ）是离散傅里叶变换的快速算法，它是根据离散傅里

叶变换的奇、偶、虚、实等特性，对离散傅里叶变换的算法进行改进获得的，它对离散傅里叶变换并没有新的发现。

有限长序列x(n)及其频域表示X(k)可由以下离散傅立叶变换得出

0x(k)=DFT[x(n)]=()N nk

n x n W -=∑ 01k N ≤≤- （1） 10

1x(n)=IDFT[X(k)]=()N nk

n X k W N --=∑ 01k N ≤≤- （2）其中2j

nk nk N

π-=。式（1）称为离散傅立叶正变换，式（2）称为离散傅立叶逆

变换，x(n)与X(k)构成了离散傅立叶变换对。

根据上述公式，计算一个X(k)，需要N 次复数乘法和N-1次复数加法，而计算全部X(k)( 01k N ≤≤-)，共需要2N 次复数乘法和N(N-1)次复数加法。实现一次复数乘法需要四次实数乘法和两次实数加法，一次复数加法需要两次实数加法，因此直接计算全部X(k)共需要42N 次实数乘法和2N(2N-1)次实数加法。当N 较大时，对实时信号处理来说，对处理器计算速度有十分苛刻的要求，于是如何减少计算离散傅里叶变换运算量的问题变得至关重要。

为减少运算量，提高运算速度，就必须改进算法。计算DFT 过程中需要完成的运算的系数里，存在相当多的对称性。通过研究这种对称性，可以简化计算过程中的运算，从而减少计算DFT 所需的时间。

如前所述，N 点的DFT 的复乘次数等于2N 。显然，把N 点的DFT 分解为几个较短的DFT ，可是乘法的次数大大减少。另外，旋转因子m N W 具有明显的周期性和对称性，其周期为：

22()j

m lN j

m m lN m N

W π

π-+-+=== 其对称性表现为：

m N m N N W W --=或*[]N m m

N N

W W -= FFT 算法就是不断地把长序列的DFT 分解成几个短序列的DFT ，并利用m

N W 的周期性和对称性来减少DFT 的运算次数。

W 具有以下固有特性： (1)nk N W 的周期性：()(N nk n N k n k N N N W W W ++==）

(2)nk N W 的对称性：()nk nk n n N k N N N W W W --==（） (3)nk N W 的可约性：/,n n N N n N Nn

W W W W == 另外，/2(/2)1,N k N k N N N

W W W +=-=-。利用nk N W 的上述特性，将x(n)或X(k)序列按一定规律分解成短序列进行运算，

这样可以避免大量的重复运算，提高计算DFT 的运算速度。算法形式有很多种，但基本上可以分为两大类，即按时间抽取(Decimation In Time ，DIT)FFT 算法和按频率抽取(Decimation In Frequency ，DIF)FFT 算法。

2.2 基-2FFT 算法

如果序列x(n)的长度2M N =，其中M 是整数(如果不满足此条件，可以人为地增补零值点来达到)，在时域上按奇偶抽取分解成短序列的DFT ，使最小DFT 运算单元为2点。通常将FFT 运算中最小DFT 运算单元称为基(radix)，因而把这种算法称为基-2时间抽取FFT(DIT-FFT)算法[4]

。

将x(n)按n 为奇偶分解成两个子序列，当n 为偶数时，令n=2r ；当n 为奇数时，令n=2r+l ；可得到

12(2)(),(21)(),0,...,

x r x r x r x r r =+==- （3）则其DFT 可写成

112

2(21)0

()(2)(21)N N rk r k

N N

r r X k x r W x r W --+===++∑∑ 11222(21)1

()()N

N rk r k

r r x r W

x r W --+===

+∑∑

112

/20

(2)(21)N N

rk rk

N N r r x r W

x r W --===

++∑∑ 12()()k

N X k W X k =+ （4）

1()X k 和2()X k 均分别是N/2点序列1()x n 和2()x n 的DFT ，而且r 与k 的取值满足0,1,…,N/2-1。而X(k)是一个N 点的DFT ，因此式（4）只计算了X(k)的前N/2

的值。由DFT 和nk

N W 的性质可得到X(k)的后N/2的值为：

212()()()222

k N N N N X k X k W X k ++=+++12()()k

N X k W X k =+ （5）

式（4）和式（5）表明，只要计算出两个N/2点的DFT 1()x k 和2()x k ，经过线性组合，即可求出全部N 点的X(k)。由于2M N =，1/22M N -=仍为偶数，因而这样的分解可以继续进行下去，直到最后的单元只需要做2点DFT 为止。

若Xm(p)和Xm(q)为输入数据，1()m X p +和1()m X q +为输出数据，nk

N W 为旋转

因子，则对于基-2DIT-FFT 算法，蝶形运算的基本公式为

11()()()()()(){

m m m N

m m m N X p X p X q W X q X p X q W ++=+=-

其图形表示如图1所示，称Xm(P)为上结点，Xm(q)为下结点。

图1 时间抽取蝶形运算单元

对于一个8点的FFT ，根据上述算法可以得到一个完整的N=8的基-2DIT-FFT 的运算流图，如图2所示。

图2 N=8 DIT-FFT 运算流图

根据上述算法原理及运算流图，可以得出基-2DIT-FFT 的基本特点，特点如下。

(1)级数分解：对于2M N =。共分了M 级，每级包含N/2个蝶形运算单元，总共所需蝶形运算个数为

2log 22

N N

M N ?=?。 (2)运算量估计：每个蝶形运算需要一次复数乘法和两次复数加(减)法，N 点FFT 共需要

2log 2

N ?次复数乘法，2log N N ?次复数加(减)法。实际上有些蝶形运算不需要做复乘。

(3)原位运算：当数据输入到存储器以后，每一组蝶形运算后，结果仍然存放在这同一组存储器中的同一位置，不需要另辟存储单元，直到最后输出。 (4)位码倒序：由图2可以看到，FFT 输出的X(k)的次序正好是顺序排列的，即X(0)，X(1)，…,X(7)，而输入X(n)是按x(0)，X(4)，…,X(7)的倒序存入存储单元，即为倒序输入，正序输出。这种顺序看起来相当杂乱，然而它是有规律的，即位码倒序规则。

(5)旋转因子的确定：由8点FFT 的三次迭代运算可以看出k N W 的变化。在第一级迭代中，只有一种类型的蝶形运算系数，即08W ，参加蝶形运算的两个数据点间隔为l ；在第二级迭代中，有两种类型的蝶形运算系数，分别是08W 和28W ，参加蝶

形运算的两个数据点间隔为2；在第三级迭代中，有四种类型的蝶形运算系数，分别是08W ，18W ，28W ，38W ，参加蝶形运算的两个数据点间隔4。可见，每次迭代的蝶形类型比前一迭代增加一倍，间隔也增大一倍。最后一次迭代的蝶形类型最多，参加蝶形运算的两个数据点的间隔也最大，为N/2。

3 FFT算法的DSP的实现

3.1 基于DSP实现FFT变换频谱分析

3.1.1DSP芯片和编程工具CCS 2.0的简介

(1)TMS320C5402简介

TMS320C5402是TI公司为了实现低功耗、高性能而专门设计的定点DSP 芯片。它有如下的特点：具有运算速度快，指令周期可以达10ns以内；优化的CPU结构，内部有1个40位的算术逻辑单元，2个40位的累加器，2个40位加法器，1个17×17的乘法器和40位的桶形移位器。有4条内部总线和2个地址产生器。先进的DSP结构可以高效快速实现数字信号处理中的各种算法的运算。它不仅具有标准的串行口和时分复用(TDM)串行口，还提供了自动缓冲串行DBSP和与外部处理器通信的HPI的主机接口。HPI可以与外部标准的微处理器直接接口。

(2)CCS2.0简介

DSP编程工具CCS是继“一体化的DSP解决方案"后，TI公司为了巩固其在DSP业界的地位而在开发工具方面的一次重拳出击。CCS集成了开发环境，使得DSP代码开发过程从编程、编译到调试代码的性能测试都集成在一个环境下进行，而且各项功能都有了一定程度的提升，简化了开发过程，该工具主要集成了以下几个软件工具：

(1)DSP代码产生工具(包括C编译器、汇编优化器、汇编器和连接器)。CCS 不仅支持高级语言C编程、汇编语言编程，还支持高级语言C汇编语言混合模式编程，降低了代码开发难度；

(2)软件模拟器(SIMULATOR)。模拟整个硬件的开发过程，使得系统的实现更加可靠；

(3)实时基础软件。DSP/BIOS和主机目标机之间的实时数据交换软件RTDX，它们所提供的实时分析功能为目标系统提供了一个实时窗口，不仅可以直接实时显示原始数据，还可以对原始数据进行处理。在传统的主机调试器必须通过在应用程序中插入断点，中断应用程序运行才能与目标系统交换数据，这种方法不仅麻烦，而且所得到的数据只是应用程序在高速运行中的一个侧面，为故障诊断和

系统性能评测等带来了许多不便。利用RTDX技术，就可以在不中断应用程序的前提下完成主机与目标机之间的实时数据交换，另外RTDX完成主机与目标机数据交换所使用的是DSP内部的仿真逻辑和JTAG接口，它不占用DSP系统的总线、串口等I/O资源，所以可以在应用程序背景下运行对DSP系统的影响很小[7-8]。

3.2.2利用DSP中的FFT函数进行频谱分析

（1）启动CCS，在CCS中建立一个C源文件和一个命令文件，并将这两个文件添加到工程，再编译并装载程序：

阅读Dsp原理及应用中fft 用dsp实现的有关程序。

双击，启动CCS的仿真平台的配着选项。选择C5502 Simulator。

（3）启动ccs2后建立工程文件FFT.pjt

（4）建立源文件FFT.c与链接文件FFT.cmd

（5）将这两个文件加到FFT.pjt这个工程中。

（6）创建out文件

（7）加载out文件

（8）加载数据

（9）改变参数

用View/Graph/Time/Frequency打开一个图形观察窗口；设置该图形窗口变量及参数，采用双踪观察在启动地址分别为0x3000H和0x3080H，长度为128的单元中数值的变化，数值类型为16位有符号整型变量，这两段存储单元中分别

存放的是经A/D转换后的语音信号和对该信号进行FFT变换的结果。

单击“Animate”（或按F10）运行程序；调整观察窗口并观察输入信号波形及其FFT变换结果；单击“Halt”暂停程序运行，可以得到语音信号的时域波形和对该信号进行FFT变换谱分析的静态图像

频谱分析结果（1）

频谱分析结果（2）

两图分别为输入语音信号频率大小不同情况下的结果；其中中上面的波形为语音信号的时域波形，下面的波形为对该信号进行FFT变换后的谱分析结果。

由此我们可以得出：数字信号处理（DSP）能够对信号进行实时分析，以便我们对各种信息能够更及时的了解，这也是它的优越性所在，使得他在我们的生

活生产中有着更广泛的应用。

结论

本论文学习和研究了离散傅里叶变换（DFT）和快速傅里叶变换（FFT）的算法，把重点放在了时间抽取法基-2FFT算法上。以及在DSP基础上用FFT变换对信号进行频谱分析。明确了FFT在DSP芯片上的实现的关键。基于DSP的快速傅里叶变换频谱分析的研究使FFT能够有效的在DSP芯片上实现，有助于我们能够更及时的了解信息，对我们的生活生产以及科技研究有很大的帮助。

自从快速傅里叶变换（FFT）出现以后，频谱分析技术便很快的发展起来，而且越来越贴近我们的生活生产，如医疗器械，无线电通信等等。但是我们对频谱分析技术的研究并未达到最高的层次，未来发展具有很广阔的空间。

参考文献

[1] 赵红怡.DSP 技术与应用实例[M] .电子工业出版社，2008.

[2] 汪安民.TMS320C54xx DSP 实用技术[M] .清华大学出版社，2007.

[3] 戴明帖,周健江.TMS320C54x DSP 结构原理及应用[M].北京航空航天大学出版社，2007.

[4] 清源科技.TMS320C54x DSP 应用程序设计教程[M] .机械工业出版社，2004.

[5] 潘松，黄继业.EDA 技术实用教程[M].科学出版社，2002.

[6] 万山明，《TMS320F281xDSP原理与应用实例》，北京航空航天大学出版社，2007.7.

[7] 戴明桢等编著.TMS320C54X DSP 结构原理及应用. 北京：航空航天大学出版社，第2版，2007.

[8] 彭启琮编著.DSP技术的发展与应用.北京：高等教育出版社，2002.

[9] 乔瑞平，崔涛，张芳娟.TMS320C54X DSP原理及应用[M]. 西安：西安电子科技大学出版社，2010 .

[10] 刘益成.DSP应用程序设计与开发[M].北京：北京航空航天大学出版社，2002.

附录

Cmd源文件代码：

-f 0

-w

-stack 500

-sysstack 500

-l rts55.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:{}>DARAM2

.far:{}>DARAM2

.const:{}>DARAM2

.switch:{}>DARAM2

.sysmem:{}>DARAM2

.cio:{}>DARAM2

.MEM$obj:{}>DARAM2

.sysheap:{}>DARAM2

.sysstack:{}>DARAM2

.stack:{}>DARAM2

.input:{}>DARAM2

.fftcode:{}>DARAM2

}

C文件源码：

#include "math.h"

#define sample_1 256

#define signal_1_f 60

#define signal_2_f 200

#define signal_sample_f 512

#define pi 3.1415926

int input[sample_1];

float fwaver[sample_1],fwavei[sample_1],w[sample_1]; float sin_tab[sample_1];

float cos_tab[sample_1];

void init_fft_tab();

void input_data();

void fft(float datar[sample_1],float datai[sample_1]); void main()

{

int i;

init_fft_tab();

input_data();

for (i=0;i

{

fwaver[i]=input[i];

fwavei[i]=0.0f;

w[i]=0.0f;

}

fft(fwaver,fwavei);

while(1);

}

void init_fft_tab()

{

float wt1;

float wt2;

int i;

for (i=0;i

{

wt1=2*pi*i*signal_1_f;

wt1=wt1/signal_sample_f;

wt2=2*pi*i*signal_2_f;

wt2=wt2/signal_sample_f;

input[i]=(cos(wt1)+cos(wt2))/2*32768;

}

void input_data()

{

int i;

for(i=0;i

{

sin_tab[i]=sin(2*pi*i/sample_1);

cos_tab[i]=cos(2*pi*i/sample_1);

}

void fft(float datar[sample_1],float datai[sample_1])

{

int x0,x1,x2,x3,x4,x5,x6,x7,xx;

int i,j,k,b,p,L;

float TR,TI,temp;

for(i=0;i

{

x0=x1=x2=x3=x4=x5=x6=0;

x0=i&0x01;x1=(i/2)&0x01;x2=(i/4)&0x01;x3=(i/8)&0x01;

x4=(i/16)&0x01;x5=(i/32)&0x01;x6=(i/64)&0x01;x7=(i/128)&0x01;

xx=x0*128+x1*64+x2*32+x3*16+x4*8+x5*4+x6*2+x7;

datai[xx]=datar[i];

}

for(i=0;i

{

datar[i]=datai[i];datai[i]=0;

}

for(L=1;L<=8;L++)

{

b=1;i=L-1;

while(i>0)

{

b=b*2;i--;

}

for(j=0;j<=b-1;j++)

{

p=1;i=8-L;

while(i>0)

{

p=p*2;i--;

}

p=p*j;

for(k=j;k<256;k=k+2*b)

{

TR=datar[k];TI=datai[k];temp=datar[k+b];

datar[k]=datar[k]+datar[k+b]*cos_tab[p]+datai[k+b]*sin_tab[p];

datai[k]=datai[k]-datar[k+b]*sin_tab[p]+datai[k+b]*cos_tab[p];

datar[k+b]=TR-datar[k+b]*cos_tab[p]-datai[k+b]*sin_tab[p];

datai[k+b]=TI+temp*sin_tab[p]-datai[k+b]*cos_tab[p];

}

for(i=0;i

{

w[i]=sqrt(datar[i]*datar[i]+datai[i]*datai[i]);

}