串行FFT递归算法蝶式递归计算原理求傅里叶变换.docx
《串行FFT递归算法蝶式递归计算原理求傅里叶变换.docx》由会员分享，可在线阅读，更多相关《串行FFT递归算法蝶式递归计算原理求傅里叶变换.docx（16页珍藏版）》请在冰豆网上搜索。
串行FFT递归算法蝶式递归计算原理求傅里叶变换
串行FFT递归算法（蝶式递归计算原理）求傅里叶变换
摘要
FFT，即为快速傅氏变换，是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。
它对傅氏变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。
设x（n）为N项的复数序列，由DFT变换，任一X（m）的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”（四次实数乘法和四次实数加法），那么求出N项复数序列的X（m）,即N点DFT变换大约就需要N^2次运算。
当N=1024点甚至更多的时候，需要N2=1048576次运算，在FFT中，利用WN的周期性和对称性，把一个N项序列（设N=2k,k为正整数），分为两个N/2项的子序列，每个N/2点DFT变换需要（N/2）^2次运算，再用N次运算把两个N/2点的DFT变换组合成一个N点的DFT变换。
这样变换以后，总的运算次数就变成N+2（N/2）^2=N+N^2/2。
继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。
而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的DFT变换就只需要Nlog
（2）（N）次的运算，N在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是FFT的优越性。
关键字：
FFT蝶式计算傅里叶变换
一．题目及要求
1.1题目
对给定的
，利用串行FFT递归算法（蝶式递归计算原理）计算其傅里叶变换的结果。
1.2要求
利用串行递归与蝶式递归原理，对给定的向量求解傅里叶变换的结果。
二．设计算法、算法原理
2.1算法原理与设计
蝶式递归计算原理：
令为n/2次单位元根，则有，
将b向量的偶数项和奇数项分别记
为和。
注意推导中反复使用：
。
图2.1公式图形
2.2设计步骤
对于以上的分析可画出如图2.2所示的离散傅里叶变换递归计算流图。
图2.3就是一个按此递归方法计算的n=8的FFT蝶式计算图。
FFT的蝶式递归计算图（有计算原理推出）：
图2.2递归计算流图
特别的，ｎ＝８的FFT蝶式计算图（展开的）：
图2.3蝶式计算图
按输入元素
展开，前面将输出序列之元素
按其偶下标（
）和（
）展开，导出
和
递归计算式，按此构造出了如图1所示的FFT递归计算流程图。
事实上，我们也可以将输入序列之元素
按其偶下标（
）和几下标（
）展开，则导出另一种形式的FFT递归计算式
。
三．算法描述、设计流程
3.1算法描述
SISD上的FFT分治递归算法：
输入:
a=（a0,a1,…,an-1）;输出:
B=（b0,b1,…,bn-1）
ProcedureRFFT（a,b）
begin
ifn=1thenb0=a0else
（1）RFFT（a0,a2,…,an-2,u0,u1,…,un/2-1）
（2）RFFT（a1,a3,…,an-1,v0,v1,…,vn/2-1）
（3）z=1
（4）forj=0ton-1do
（4.1）bj=ujmodn/2+zvjmodn/2
（4.2）z=zω
endfor
endif
end
注:
（1）算法时间复杂度t（n）=2t（n/2）+O（n）
t（n）=O（nlogn）
n=8的FFT蝶式计算图:
图3.1FFT蝶式计算图
n=6的FFT递归计算流程图:
图3.2FFT递归计算流程图
3.2流程图
飞
是
否
是
否
是
否
四．源程序代码及运行结果
4.1源程序代码
/************FFT***********/
#include//整个程序输入和输出利用同一个空间x[N]，节约空间
#include
#include
#defineN1000//定义输入或者输出空间的最大长度
typedefstruct
{
doublereal;
doubleimg;
}complex;//定义复数型变量的结构体
voidfft（）;//快速傅里叶变换函数声明
voidinitW（）;//计算W（0）~W（size_x-1）的值函数声明
voidchange（）;//码元位置倒置函数函数声明
voidadd（complex,complex,complex*）;/*复数加法*/
voidmul（complex,complex,complex*）;/*复数乘法*/
voidsub（complex,complex,complex*）;/*复数减法*/
voiddivi（complex,complex,complex*）;/*复数除法*/
voidoutput（）;/*输出结果*/
complexx[N],*W;/*输出序列的值*/
intsize_x=0;/*输入序列的长度，只限2的N次方*/
doublePI;//pi的值
intmain（）
{
inti;
system（"cls"）;
PI=atan
（1）*4;
printf（"Pleaseinputthesizeofx:
\n"）;/*输入序列的长度*/
scanf（"%d",&size_x）;
printf（"Pleaseinputthedatainx[N]:
（suchas:
56）\n"）;
for（i=0;iscanf（"%lf%lf",&x[i].real,&x[i].img）;
initW（）;//计算W（0）~W（size_x-1）的值
fft（）;//利用fft快速算法进行DFT变化
output（）;//顺序输出size_x个fft的结果
return0;
}
/*进行基-2FFT运算,蝶形算法。
这个算法的思路就是，先把计算过程分为log（size_x）/log
（2）-1级（用i控制级数）；
然后把每一级蝶形单元分组（用j控制组的第一个元素起始下标）；
最后算出某一级某一组每一个蝶形单元（用k控制个数，共l个）。
*/
voidfft（）
{
inti=0,j=0,k=0,l=0;
complexup,down,product;
change（）;//实现对码位的倒置
for（i=0;i（2）;i++）//循环算出fft的结果
{
l=1<for（j=0;j{
for（k=0;k{//算出j组中第k个蝶形单元
mul（x[j+k+l],W[（size_x/2/l）*k],&product）;/*size/2/l是该级W的相邻上标差，l是该级该组取的W总个数*/
add（x[j+k],product,&up）;
sub（x[j+k],product,&down）;
x[j+k]=up;//up为蝶形单元右上方的值
x[j+k+l]=down;//down为蝶形单元右下方的值
}
}
}
}
voidinitW（）//计算W的实现函数
{
inti;
W=（complex*）malloc（sizeof（complex）*size_x）;
/*申请size_x个复数W的空间（这部申请的空间有点多，实际上只要申请size_x/2个即可）*/
for（i=0;i<（size_x/2）;i++）
/*预先计算出size_x/2个W的值，存放,由于蝶形算法只需要前size_x/2个值即可*/
{
W[i].real=cos（2*PI/size_x*i）;//计算W的实部
W[i].img=-1*sin（2*PI/size_x*i）;//计算W的虚部
}
}
voidchange（）//输入的码组码位倒置实现函数
{
complextemp;
unsignedshorti=0,j=0,k=0;
doublet;
for（i=0;i{
k=i;
j=0;
t=（log（size_x）/log
（2））;
while（（t--）>0）
{
j=j<<1;
j|=（k&1）;
k=k>>1;
}
if（j>i）
{
temp=x[i];
x[i]=x[j];
x[j]=temp;
}
}
}
voidoutput（）//输出结果实现函数
{
inti;
printf（"Theresultareasfollows\n"）;
for（i=0;i{
printf（"%.4f",x[i].real）;//输出实部
if（x[i].img>=0.0001）//如果虚部的值大于0.0001，输出+jx.img的形式
printf（"+j%.4f\n",x[i].img）;
elseif（fabs（x[i].img）<0.0001）
printf（"\n"）;
else
printf（"-j%.4f\n",fabs（x[i].img））;
//如果虚部的值小于-0.0001，输出-jx.img的形式
}
}
voidadd（complexa,complexb,complex*c）//复数加法实现函数
{
c->real=a.real+b.real;//复数实部相加
c->img=a.img+b.img;//复数虚部相加
}
voidmul（complexa,complexb,complex*c）//复数乘法实现函数
{
c->real=a.real*b.real-a.img*b.img;//获取相乘结果的实部
c->img=a.real*b.img+a.img*b.real;//获取相乘结果的虚部
}
voidsub（complexa,complexb,complex*c）//复数减法实现函数
{
c->real=a.real-b.real;//复数实部相减
c->img=a.img-b.img;//复数虚部相减
}
voiddivi（complexa,complexb,complex*c）//复数除法实现函数
{
c->real=（a.real*b.real+a.img*b.img）/（b.real*b.real+b.img*b.img）;
//获取相除结果的实部
c->img=（a.img*b.real-a.real*b.img）/（b.real*b.real+b.img*b.img）;
//获取相除结果的虚部
}
4.2运行结果
（1）处理器p=8：
图4.1当
时串行FFT输出结果
（2）处理器p=8：
当
时输出结果与计算结果相符如图4.2所示
图4.2运行图
五．算法分析、优缺点
5.1算
升级会员 