串行FFT递归算法蝶式递归计算原理求傅里叶变换Word文件下载.docx
- 文档编号:15842136
- 上传时间:2022-11-16
- 格式:DOCX
- 页数:16
- 大小:273.41KB
串行FFT递归算法蝶式递归计算原理求傅里叶变换Word文件下载.docx
《串行FFT递归算法蝶式递归计算原理求傅里叶变换Word文件下载.docx》由会员分享,可在线阅读,更多相关《串行FFT递归算法蝶式递归计算原理求傅里叶变换Word文件下载.docx(16页珍藏版)》请在冰豆网上搜索。
注意推导中反复使用:
。
图2.1公式图形
2.2设计步骤
对于以上的分析可画出如图2.2所示的离散傅里叶变换递归计算流图。
图2.3就是一个按此递归方法计算的n=8的FFT蝶式计算图。
FFT的蝶式递归计算图(有计算原理推出):
图2.2递归计算流图
特别的,n=8的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<
stdio.h>
//整个程序输入和输出利用同一个空间x[N],节约空间
math.h>
stdlib.h>
#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);
Pleaseinputthedatainx[N]:
(suchas:
56)\n"
for(i=0;
i<
size_x;
i++)/*输入序列对应的值*/
%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();
//实现对码位的倒置
log(size_x)/log
(2);
i++)//循环算出fft的结果
l=1<
<
i;
for(j=0;
j<
j+=2*l)
for(k=0;
k<
l;
k++)//算出第i级内j组蝶形单元的结果
{//算出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的实现函数
W=(complex*)malloc(sizeof(complex)*size_x);
/*申请size_x个复数W的空间(这部申请的空间有点多,实际上只要申请size_x/2个即可)*/
(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;
k=i;
j=0;
t=(log(size_x)/log
(2));
while((t--)>
0)
j=j<
1;
j|=(k&
1);
k=k>
>
if(j>
i)
temp=x[i];
x[i]=x[j];
x[j]=temp;
voidoutput()//输出结果实现函数
Theresultareasfollows\n"
%.4f"
x[i].real);
//输出实部
if(x[i].img>
=0.0001)//如果虚部的值大于0.0001,输出+jx.img的形式
+j%.4f\n"
x[i].img);
elseif(fabs(x[i].img)<
0.0001)
else
-j%.4f\n"
fabs(x[i].img));
//如果虚部的值小于-0.0001,输出-jx.img的形式
voidadd(complexa,complexb,complex*c)//复数加法实现函数
c->
real=a.real+b.real;
//复数实部相加
img=a.img+b.img;
//复数虚部相加
voidmul(complexa,complexb,complex*c)//复数乘法实现函数
real=a.real*b.real-a.img*b.img;
//获取相乘结果的实部
img=a.real*b.img+a.img*b.real;
//获取相乘结果的虚部
voidsub(complexa,complexb,complex*c)//复数减法实现函数
real=a.real-b.real;
//复数实部相减
img=a.img-b.img;
//复数虚部相减
voiddivi(complexa,complexb,complex*c)//复数除法实现函数
real=(a.real*b.real+a.img*b.img)/(b.real*b.real+b.img*b.img);
//获取相除结果的实部
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算
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 串行 FFT 递归 算法 计算 原理 傅里叶变换