课程设计(论文)_快速傅里叶变换程序设计

1 设计任务描述
1.1 设计题目
快速傅里叶变换程序设计
设计要求
1.2.1 设计目的
1）理解FFT的算法以及利用DSP实现的方法。

2）能熟练的调试程序并能观察其结果。

3）熟悉TMS320C54x系列DSP芯片的软件设计方法。

基本要求
1）研究FFT原理以及利用DSP实现的方法。

2）编写FFT程序。

3）调试程序，观察结果。

2 设计思路
2.1 FFT 算法原理
若给定由N 个信号样本{x （0），x （1），…，x （N -1）}组成的信号序列x （n ），DFT 可用式2-1给出：
1
()()N nk
N
n X k x n W
-==
∑ k =0,1，…，N -1
（2-1）
式2-1中，nk N W 称为旋转因子或蝶形因子，nk
N W =2/j nk N
e
π-。

从中可以看出：
当信号样本为复数时，计算单个()X k 需经过N 次复数乘法和N -1次复数加法运算，相当于4N 次实数乘法和2（2N -1）次实数加法。

完成全部N 点DFT 共需
2N 次复数乘法和N （N -1）复数加法运算。

可见，随着N 不断增加，整个DFT
运算量是相当庞大的，而FFT 算法通过对计算过程的深入分析，利用旋转因子
nk
N
W 具有的周期性与对称性，实现了降低运算复杂度的目的。

当序列长度N 为偶数时，信号序列x （n ）可被分解为奇、偶两个子序列，相应的N 点DFT 被分解为两个N /2点的DFT ：
()()()k
N X k G k W H k =+ k =0,1, …，N /2-1
（2-2）
(/2)()()k N X N k G k W H k +=- k =0,1, …，N /2-1
（2-3）
式（2-2）和（2-3）中，G(k)和()H k 分别表示x （n ）分解后得到的N /2点偶序列点奇序列的DFT 。

式（2-2）和式（2-3）表明，只要求出G(k)和()H k ，x （n ）前N /2点和后N /2点的DFT 就得到了，整个序列的DFT 也就得到了。

这样做的好处是计算N 点DFT 只需要约2N /2次复数乘法，总运算量约为直接DFT 运算量的一半
同理，当N /2为偶数时，每个N /2点的DFT 又可被分解成两个N /4点的DFT ，进一步减少了DFT 运算的复杂度。

依次类推，直到不能继续分解为止。

分解结束时，最小DFT 的点数称为称为基数，当N =2L （L 为正整数）时，经
过L -1次分解，N 点DFT 最终可被分解为N /2个两点的DFT ，即得到基数为2
的FFT 运算，使得DFT 所需复数乘法次数降至2(/2)log N
N 。

2.2 基2FFT 的蝶形运算流图
基2FFT 的蝶形运算过程可用图2-1所示，此时N =8，L =2log N
=3。

图2-1 8点基2FFT 运算过程
观察图2-1，根据DFT 的基2FFT 算法，可以总结出以下几条规律：
（1）N 点FFT 运算从输入端开始，逐级进行，共需经过L 级运算；在第m （m =1，2，…，M ）级中存在2L m -个相似的蝶形运算组（除输入数据不同外）；每个组内蝶形运算的个数为12m -，参与每个蝶形运算的两个输入数据相距12m -个点。

（2）中间数据的存储，可采用原位存储法。

即每次蝶形运算的结果存储在与原数据相同的内存单元内。

（3）为了保证输出数据按自然数序排列，在进行FFT 之前输入数据需要按照特定的顺序存放，通过位倒序寻址可以满足这种要求。

2.2.1 输入、输出序列的倒位序规律
由流程图2-1可以看出，当进行原位运算时，发现当运算完成后，FFT 的输出()X k 按自然顺序排列在存储单元中，即按(0)X ，(1)X ，…，(7)X 的顺序排
列；但是这是输入()x n 却不是按自然顺序存储的，而是按(0)x ，(4)x ，…，(7)x 的顺序存入存储单元，这种方式就称之为倒位序。

当用二进制表示这个顺序时，它正好是“码位倒置”的顺序。

例如，原来的自然顺序应是(1)x 的地方，现在存放(4)x ，用二进制码表示这一规律时，则是在
(001)x 处存放(100)x ，(011)x 出存放(110)x ，即将自然顺序的二进制码位倒置过
来，第一位码变成最末位码，这样倒置以后的顺序正是输入所需要的顺序。

表2-1中列出的是N =8时按码位倒置规律所得的顺序，其结果与按时间抽取算法FFT 流程图中的输入顺序是一致的。

同理，当N =256时，亦可以采用同样的方法进行位倒序操作。

表2-1 码位倒置顺序
2. 蝶距的计算
设N =2L ，则整个运算流图中包含L 级蝶形运算，每一级则有/2N 个蝶形单元。

蝶距即每个蝶形单元两个输入（出）节点的序号差。

以8N =为例，结合图2-1可知共包含3级蝶形运算，每一级蝶形单元的蝶距如下：第一级，蝶距为1，可以看作由11022-=得到：第二级，蝶距为2，可以看作由21122-=得到；第三级，蝶距为4，可以看作由31222-=得到。

因此得：第m 级蝶形单元的蝶距为：12m -。

2. 旋转因子的计算
由FFT 算法原理过程可知，若N =2L ，则共有L 级蝶形运算，各级蝶形运算
中旋转因子分别如下：第L 级的旋转因子为r
N W （r =0,1，…，/21N -）；第L -1
级的旋转因子为/2r
N W （r =0,1，…，2/21N -）；…；第一级的旋转因子为1
/2L r N W -（r =0,1，…，/21L N -）。

由此可见， 第m 级蝶形运算中旋转因子为/2L M
r
N W -，r =0，1，…，1/21L M N -+-。

2.3 FFT 算法的DSP 的实现方法
设FFT 运算的输入数据为实数，则2N 点实数FFT 算法的实现步骤为：第一步，把2N 实数输入序列组合成N 点的复数序列。

然后把该复数序列进行位倒序操作后存储在输入区。

第二步，进行N 的FFT 运算。

第三步，把N 点FFT 输出拆成2N 的复数序列，这2N 的复数序列对应于2N 点时实数输入序列的DFT 输出。

第四步，结果输出及功率谱计算。

3 软件流程图
4 各部分程序设计及参数计算4.1 初始化程序
#include"math.h"
#define N 128
void InitForFFT();
void MakeWave();
void finv(int N1,float *xr,float *xi);
int INPUT[N],DATA[N];
float fWaveR[N],fWaveI[N],w[N];
float sin_tab[N],cos_tab[N];
int Mum;
初始化程序是对采样点数N、各个子程序、数据实部和虚部、输入输出以及正余弦表进行定义，在后面的程序中运用时就不用再进行定义了。

4.2 子程序
旋转因子、蝶距运算
void FFT(float Xr[N],float Xi[N])
{
int S,B;
int m,j,k;
float X,Y;
finv(N,Xr,Xi);
for(m=1;m<=Mum;m++)
{
B=(int)(pow(2,m-1)+0.5);
for(j=0;j<B;j++)
{
S=j*(int)(pow(2,Mum-m)+0.5);
for(k=j;k<=N-1;k+=(int)(pow(2,m)+0.5))
该子程序是时间抽取法FFT程序，要求采样点数为2的整数幂次方，Xr[]和Xi[]分别为输入序列的实部和虚部。

S为旋转因子的幂数，B为蝶形图运算输入数据的距离，也即各级旋转因子的个数。

finv(N,Xr,Xi)语句是实现倒序运算函数，对输入序列倒序。

B=(int)(pow(2,m-1)+0.5)语句计算出蝶距 B=2^(m-1)，语句计算出旋转因子的幂数S=2^(Num-1)。

4.2.2 蝶形图变换
{
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;
}
该段子程序是实现蝶形图的变换，蝶形运算展开，结果的实部和虚部分别存储在原实部和虚部位置。

4.2.3 波形发生函数
void MakeWave()
{
int i;int j;int l;
for(i=0;i<N;i++)
{j=(int)i/20;
l=j%2;
if(l==0)INPUT[i]=1024;
else INPUT[i]=0;}
}
该程序是实现输入波形为方波，从0点开始20个点一个高电平，20个点一个低电平。

4.2.4 倒序
void finv(int N1,float *xr,float *xi)
{
int m,n,N2,k;
float T;
N2=N1/2;
n=N2;
for(m=1;m<=N1-2;m++)
{
if(m<n)
{
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=n-k;
k=(int)(k/2+0.5);
}
n=n+k;
}
}
该程序段是倒序运算函数finv(N1,Xr,Xi),对输入序列倒序。

N1为序列长度；Xr[],Xi[]分别为输入序列的实部和虚部。

倒序原理：倒序列的加1实在最高位加1，满2向次高位进1，最高位变0，依次往下，从当前倒序值可求得下一倒序值。

4.3 主程序
main()
{
int i;
InitForFFT();
MakeWave();
for(i=0;i<N;i++)
{
fWaveR[i]=INPUT[i];
fWaveI[i]=0.0;
w[i]=0.0;
}
Mum=(int)(0.5+log(N)/log(2));
FFT(fWaveR,fWaveI);
for(i=0;i<N;i++)DATA[i]=w[i];
while(1);
}
主程序的功能是调用各个子程序，Mum=(int)(0.5+log(N)/log(2))语句中Mum 为蝶形运算的级数，N=2^Num。

Num=7。

小结
在这次设计中，我完成了这次任务，在一周的课设中，我弄懂了快速傅里叶变换程序设计，并且调试出该程序的输入输出波形以及功率谱。

在这次课程设计中,我学到了很多上课时没有学到的知识。

这次设计的程序是关于快速傅里叶变换的，在设计程序之前得先弄明白快速傅里叶变换的理论知识，我翻阅了上学期学的数字信号处理，弄明白了该程序要实现的功能。

然后把整个程序分成好几个模块，先读懂每个模块，再把程序结合在一起，最后终于弄懂了整个程序。

在程序的调试过程中我遇到了许多困难。

首先，在编译时有两个错误是认为因素造成的，有一个是因为把指令的拼音打错了，还有一个是因为在程序的最开始多留出一行，由于DSP的汇编指令要求非常严格，所以在编写程序时应该注意书写的格式。

通过课程设计我收获很多，不仅对这个课程有了更深的理解，而且也学会了团队精神的重要性，个人的能力是有限的，团结才能有力量，我们都尽自己所能来完成这次课程设计。

通过这次课程设计,我懂得了理论与实际结合是很必要的。

只有理论知识是远远不够的，只有把所学的理论知识与实际结合起来，从理论中得出结论，才能真正为社会服务，从而提高自己的实际动手能力和独立思考能力。

重要的是我们在设计过程中发现了自己的不足，对前面学过的知识理解的不够深刻，掌握的不够牢固。

通过这次课程设计，我们相当于把前面所学过的知识又重新温故了一遍。

致谢
课程设计结束了，在课程设计的过程中我学到了很多课堂上所不能体会的知识。

我们做的是快速傅里叶变换程序设计，在这次设计当中，我学会了团体合作还有很多关于dsp的相关知识。

一周的课程设计结束了，以前在上课的时候老师讲解的东西以为明白就可以，但是在这次实践中，我明白了，那是远远不够的。

在老师、同学的帮助下，终于完成了这次实验。

查阅资料，请教同学，相互讨论，还好有老师和同学，一步一步的修改和讲解让我觉得非常的暖心。

可能在以后步入社会以后这种感觉可能不复存在。

所以我珍惜现在的每一次的设计。

在这里要非常感谢我的老师。

每天都来实验室为同学们解答各种问题，当我们出错误的时候，老师都会在身边所以我们感到非常的安心。

这也是一次跟老师更加近距离的学习机会。

在这里还要感谢我同组的同学以及帮助我们的同学。

在每天的课程设计中，都有你们的陪伴，在我有困难的时候，你们永远在我的身边。

非常谢谢你们对我的帮助。

这次课程设计对我以后的学习乃至工作都有很大的帮助，我对dsp有了更深一层次的了解。

虽然对程序设计还是有很多困难，很多要学的。

但我想只要我们珍惜每次设计的经过，从中我们肯定会有一些所得。

再次感谢我们的指导老师，还有曾经一直帮助我的同学。

参考文献
[1] 俞一彪，孙兵. 数字信号处理理论与应用.南京：东南大学出版社，2005
[2] 李行一. 数字信号处理.重庆：重庆大学出版社，2002
[3] 郭森楙，颜允圣. 数字信号处理器体系结构、实现与应用.北京：清华大学出版社，2005
[4]姜沫岐，许涵，俞鹏，段国强 . DSP原理与应用—从入门到提高. 北京: 机械工业出版社， 2007
[5] 俞卞章. 数字信号处理. 西安: 西北工业大学出版社， 2002
[6] 王安国. 数字信号处理基础.北京: 电子工业出版社， 2003
附录A1 程序清单#include"math.h"
#define N 128
void InitForFFT();
void MakeWave();
void finv(int N1,float *xr,float *xi);
int INPUT[N],DATA[N];
float fWaveR[N],fWaveI[N],w[N];
float sin_tab[N],cos_tab[N];
int Mum;
void FFT(float Xr[N],float Xi[N])
int S,B;
int m,j,k;
float X,Y;
finv(N,Xr,Xi);
for(m=1;m<=Mum;m++)
{
B=(int)(pow(2,m-1)+0.5);
for(j=0;j<B;j++)
{
S=j*(int)(pow(2,Mum-m)+0.5);
for(k=j;k<=N-1;k+=(int)(pow(2,m)+0.5))
{
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;
}
}
}
for(m=0;m<=N/2;m++)
{
w[m]=sqrt(Xr[m]*Xr[m]+Xi[m]*Xi[m]);
}
}
main()
{
int i;
InitForFFT();
MakeWave();
for(i=0;i<N;i++)
{
fWaveR[i]=INPUT[i];
fWaveI[i]=0.0;
w[i]=0.0;
}
Mum=(int)(0.5+log(N)/log(2));
FFT(fWaveR,fWaveI);
for(i=0;i<N;i++)DATA[i]=w[i];
while(1);
}
void InitForFFT()
{
int i;
for(i=0;i<N;i++)
{
sin_tab[i]=sin(PI*2*i/N);
cos_tab[i]=cos(PI*2*i/N);
}
}
void MakeWave()
{
int i;int j;int l;
for(i=0;i<N;i++)
{j=(int)i/20;
l=j%2;
if(l==0)INPUT[i]=1024;
else INPUT[i]=0;}
}
void finv(int N1,float *xr,float *xi) {
int m,n,N2,k;
float T;
N2=N1/2;
n=N2;
for(m=1;m<=N1-2;m++)
{
if(m<n)
{
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=n-k;
k=(int)(k/2+0.5);
}
n=n+k;
}
}
附录A2 程序图形图一输入信号时域波形
图二输入信号频域波形图三输出信号功率谱。