快速傅里叶变换的C程序实现思路分析MATLAB仿真.docx
《快速傅里叶变换的C程序实现思路分析MATLAB仿真.docx》由会员分享，可在线阅读，更多相关《快速傅里叶变换的C程序实现思路分析MATLAB仿真.docx（14页珍藏版）》请在冰点文库上搜索。
快速傅里叶变换的C程序实现思路分析MATLAB仿真
现代通信原理
设计报告
快速傅里叶变换的C程序实现思路分析MATLAB仿真
组长：
孙建东
组员：
赵玉兰
2013年12月8号
快速傅里叶变换的C程序实现思路分析MATLAB仿真
一、设计目的
快速傅里叶变换在频谱分析、信号处理方面应用十分广泛，在实际实现中常使用DSP器件作为硬件平台，利用C程序设计实现的信号的快速傅里叶变换，并在频域中完成对信号的处理。
故如何实现用C程序实现快速傅里叶变换（FFT）便成为了十分重要的任务。
该次仿真设计以MATLAB为软件仿真平台用C语言的编程思路来编写MATLAB程序完成任意点数的FFT运算，并通过编写相应的仿真测试文件，通过调用MATLAB内部自带的FFT与调用自带的FFT进行对比，验证该FFT的C程序设计思路思路是否正确。
二、FFT的基本原理
FFT是一种DFT的高效算法，称为快速傅立叶变换（fastFouriertransform）。
FFT算法可分为按时间抽取算法和按频率抽取算法，简要介绍FFT的基本原理。
从DFT运算开始，说明FFT的基本原理。
DFT的运算为：
[1]
式中
由这种方法计算DFT对于X（K）的每个K值，需要进行4N次实数相乘和（4N-2）次相加，对于N个k值，共需N*N乘和N（4N-2）次实数相加。
改进DFT算法，减小它的运算量，利用DFT中
的周期性和对称性，使整个DFT的计算变成一系列迭代运算，可大幅度提高运算过程和运算量，这就是FFT的基本思想。
FFT基本上可分为两类，时间抽取法和频率抽取法，而一般的时间抽取法和频率抽取法只能处理长度N=2^M的情况，另外还有组合数基四FFT来处理一般长度的FFT
基于时间抽取的FFT算法
设N点序列x（n）,
将x（n）按奇偶分组，公式如下图
改写为：
一个N点DFT分解为两个N/2点的DFT，
继续分解，迭代下去，其运算量约为
[1]
其算法有如下规律
两个4点组成的8点DFT
蝶形算法
四个2点组成的8点DFT
按时间抽取的8点DFT
原位计算
当数据输入到存储器中以后，每一级运算的结果仍然储存在同一组存储器中，直到最后输出，中间无需其它存储器
序数重排
对按时间抽取FFT的原位运算结构，当运算完毕时，这种结构存储单元A
（1）、A
（2），…，A（8）中正好顺序存放着X（0），X
（1），X
（2），…，X（7），因此可直接按顺序输出，但这种原位运算的输入x（n）却不能按这种自然顺序存入存储单元中，而是按X（0），X（4），X
（2），X（6），…，X（7）的顺序存入存储单元，这种顺序看起来相当杂乱，然而它也是有规律的。
当用二进制表示这个顺序时，它正好是“码位倒置”的顺序。
蝶形类型随迭代次数成倍增加
每次迭代的蝶形类型比上一次蝶代增加一倍，数据点间隔也增大一倍
三、基于时间抽取的设计思路
通过对上述思路的进一步研究发现，FFT的每级蝶形运算都是有规律可循的
接下来我们便探究FFT的规律；
1、倒位序算法设计
在FFT的基于时间抽取的算法中，输入序列必须按倒位序进行重排
对于8点输入序列x（n）
输入原始序列下标号二进制重排后序列下标号二进制
X（0）000x（0）000
X
（1）001x（4）100
X
（2）010x
（2）010
X（3）011x（6）110
X（4）100x
（1）001
X（5）101x（5） 101
X（6）110x（3）110
X（7）111x（7）111
由上面的序列重排我们发现重排的序列与原始序列只是在脚标上按其二进制的顺序高低位进行了交换一次实现倒排；便实现倒位序的算法，只需将每个元素的组号转换成二进制按倒序排列后，再把倒序好的二进制转化成十进制作为元素的组号，并把对应原始数组存入该数组的元素中，程序实现见程序五
2、蝶形运算中旋转因子的规律
每级旋转因子及其个数
第M级个数旋转因子x（n）序列中x1（n）与x2（n）的间隔
N=8
1
2
3
N=16
1
2
3
4
……
N点
1
2
3
…
（k=0,1,2,3）
k
…
（k=0,1,…,
）
M
…
（k=0,1,2,3,…,
）
由上述的规律可得，第一级蝶形运算旋转因子只有一个，原因是重排后两点相距为1，故旋转因子只有一个，以后每级旋转因子随级数的变化按上表的规律变化
3、每级蝶形运算序列的原位运算规律
第1级
第2级
第3级
第4级
0
1
0
2
0
4
0
8
2
3
1
3
1
5
1
9
4
5
4
6
2
6
2
10
6
7
5
7
3
7
3
11
8
9
8
10
8
12
4
12
10
11
9
11
9
13
5
13
12
13
12
14
10
14
6
14
14
15
13
15
11
15
7
15
我们尝试用一个通项公式表示出这张表中所有的排列：
定义：
i:
蝶形运算的级数
a:
序列所分的组数；k：
每组中元素的个数
第1级：
i=0；a=8；k=0；
第2级：
i=1:
a=4;k=0,1;
第3级：
i=2;a=2;k=0,1,2,3;
第4级：
i=3;a=0;k=0,1,2,3,4,5,6,7;
x1（n）的通项公式为：
x1（n）=
*a+k，x2（n）的通项公式：
x2（n）=
*（a+1）+k
且a满足a=a+2；
验证：
当a=4；i=1;k=0,1;{x1（n）,x2（n）}={0,2,1,3,4,6,5,7,8,10,9,11,12,14,13,15};
对于N点M级的蝶形运算；
则i=1,2,…,M-1;k=
-1;而每级只有N/2个蝶形运算则a的最大取值就随每级k的取值而定
有了以上的FFT运算规律的分析便能轻易的实现FFT的算法的设计了。
四、用MATLAB模拟C程序设计思路编程实现
实现该FFT设计的matlab程序文件如下
程序一：
完成乘方运算
%自定义函数名称:
jiecheng（）完成乘方运算
%n是底数m是指数
functiony=jiecheng（n,m）
temp=1;%存储累乘的结果
fori=1:
m%完成m次累乘运算
temp=temp*n;%实现n^m的运算
end
y=temp;%输出乘方的结果
end
程序二：
完成复数的加减运算
%自定义函数名称:
Fujia完成复数加法
%x1，x2是参与蝶形运算的复数点,k判断复数运算是加还是减
functiony=Fujia（x1,x2,k）
ifk==0%k=0进行复数加运算
y（1,1）=x1（1,1）+x2（1,1）;%实部加法运算
y（2,1）=x1（2,1）+x2（2,1）;%虚部加法运算
elseifk==1%k=1进行复数加法运算
y（1,1）=x1（1,1）-x2（1,1）;%完成实部的加法运算
y（2,1）=x1（2,1）-x2（2,1）;%完成虚部的加法运算
end
end
程序三：
完成复数的乘法运算
%自定义Fucheng完成复数乘法（a+bj）（c+dj）=（ac-bd）+（ad+bc）j
%x是输入序列得点w是参与运算旋转因子
functiony=Fucheng（x,w）
y（1,1）=x（1,1）*w（1,1）-x（2,1）*w（2,1）;%完成实部的乘法运算
y（2,1）=x（2,1）*w（1,1）+w（2,1）*x（1,1）;%完成虚部的乘法运算
end
程序四：
完成一个蝶形运算
%自定义函数名称:
diexing完成一个蝶形单元的运算
%x1，x2参与蝶形运算的输入序列,w参与蝶形运算的旋转因子
function[x1,x2]=diexing（x1,x2,w）
%因为蝶形运算是满足原位运算所以必须将x1，x2*w的保存到暂存变量中
temp=x1;%x1暂存到temp中
temp1=Fucheng（x2,w）;%将x2与w的复数乘积保存到temp1中
x1=Fujia（temp,temp1,0）;%完成x1的原位运算
x2=Fujia（temp,temp1,1）;%完成x2的原位运算
end
程序五：
完成对一个序列按位倒序重排
%自定义函数名称:
daoxiu完成对输入序列下标号按二进制倒序重排
%x输入序列N需要重排的点数M下标号转成二进制所需要的位数也即FFT运算的级数
functiontem=daoxiu（x,N,M）
fori=0:
N-1%完成对N个点顺序重排
k=dec2bin（i,M）;%将下标号转换成二进制
forj=1:
floor（M/2）%完成倒序
t=k（1,j）;
k（1,j）=k（1,M-j+1）;%实现二进制高低位数据的交换
k（1,M-j+1）=t;
end
n=bin2dec（k）;%将将换后的二进制数转换成十进制
tem（1,n+1）=x（1,i+1）;%将原始数组的数据放到重排后的数组中
end
end
程序六：
完成N点M级的蝶形运算，即实现FFT的算法
%函数名称:
FFT_N完成N点M级蝶形运算的FFT（快速傅里叶变换）
%x:
进行FFT运算的输入序列w:
参与FFT运算的旋转因子N:
运算的点数M:
蝶形运算的级数
functionout=FFT_N（x,w,N,M）
x1=daoxiu（x（1,:
）,N,M）;%对输入序列的实部进行位倒序
x2=daoxiu（x（2,:
）,N,M）;%对输入序列的虚部进行位倒序
x0=[x1;x2];%将实部和虚部放入x0，第一行实部,第二行虚部
fori=0:
M-1%进行M级蝶形运算
k=1;a=0;
forj=1:
N/2%每级进行N/2次蝶形运算
m=（k-1）*jiecheng（2,M-i-1）;%每次参与蝶形运算所需的旋转因子
p=jiecheng（2,i）*a+k;%每次蝶形运算所需要的x1（n）
q=jiecheng（2,i）*（a+1）+k;%每次蝶形运算所需要的x2（n）
%进行的蝶形运算时原位运算
[x0（:
p）,x0（:
q）]=diexing（x0（:
p）,x0（:
q）,w（:
m+1））;%完成一次蝶形运算
ifk==jiecheng（2,i）%每级蝶形运算x1，x2的间隔
k=1;a=a+2;%切换到下一组间隔为k的输入序列x0参与蝶形迭代运算
else
k=k+1;%切换到下一组蝶形运算
end
end
end
out=x0;%将经过FFT运算的向量输出到向量out
end
五、编写仿真程序进行验证分析
程序七：
通过调用matlab内部FFT函数和自定义的FFT_N函数分析信号频谱
%基2算法的128点7级蝶形运算
%输入信号频率f0=80Hz,采样频率fs=500Hz
N=256;M=8;f0=100;f1=200;fs=1000;
x=zeros（2,N）;%开辟输入序列存储空间,数组第一行保存保存输入序列的实部,第二行保存输入序列的虚部
stepf=fs/N;n=0:
N-1;%stepf频率的步长即FFT的最小分辨率,n开辟N点的数组
t=2*pi*n/fs;n1=0:
stepf:
fs/2-stepf;%定义时域采样序列的长度,n1频谱显示的横轴
x0=sin（f0*t）+cos（f1*t）;%输入序列（实序列输入）
plot（t,x0,'g','LineWidth',2）;%显示输入序列的时域波形
title（'输入序列的时域波形x（t）'）;
figure;
x（1,:
）=x0;%将实序列放入进行FFT变换的数组的第一行
Y=abs（fft（x0））;%调用matlab的FFT函数进行频谱分析
plot（n1,abs（Y（1:
N/2））,'r','LineWidth',2）;gridon%显示matlab函数分析出的频谱
title（'调用matlab函数分析的频谱'）;
figure;
%用自定义FFT_N分析的频谱
fork=0:
N/2-1;%生成参与N点FFT运算的旋转因子
%w1=exp（-j*2*pi*k/N）;
w（1,k+1）=cos（-2*pi*k/N）;%旋转因子的实部
w（2,k+1）=sin（-2*pi*k/N）;%旋转因子的虚部
end
Y1=FFT_N（x,w,N,M）;%调用自己编写的FFT运算出频域的Y1（k）
forj=1:
N/2;
Y2（1,j）=sqrt（Y1（1,j）^2+Y1（2,j）^2）;%求个点信号的幅度
end
plot（n1,abs（Y2）,'b','LineWidth',2）;gridon%显示调用自定的FFT函数生成的频谱
title（'调用自定义的FFT_N函数分析的频谱'）;
六、仿真结果
调用matlab的FFT函数分析出的信号频谱如图1
如图1
调用自定义的FFT函数分析出的信号频谱如图2
七、结果分析与总结
本次设计顺利完成了FFT算法的C程序实际思路基于MATLAB仿真平台的设计，从编写的仿真测试程序来看，本次设计的算法很好地完成了FFT的设计任务，分析出的频谱与调用MATLAB分析出的频谱一模一样，且该仿真的FFT思路在DSP器件上的实现也可以像matlab里的FFT一样，可将其封装成子函数模块，当需要调用时只需向该子函数传递需要进行FFT变换的序列x（n）及点数N，蝶形运算的级数M，该函数就能完成FFT的变换，其返回的序列不仅包含了幅度信息，还包含了相位信息，但该设计还需要进一步探索，现在通过仿真证实DSP上是可实现的，接下来需要在DSP真正的实现该算法并考虑是否能进一步优化程序结构及存储空间的分配等。
升级会员 