FFT的计算机实现Word下载.docx
《FFT的计算机实现Word下载.docx》由会员分享，可在线阅读，更多相关《FFT的计算机实现Word下载.docx（13页珍藏版）》请在冰豆网上搜索。
一直分为奇序列和偶序列来求其值，直到奇偶序列长度为1为止。
1.2蝶形图
对于FFT运算，我们通常用蝶形图来表示
比如
我们表示为
下面以N=8为例，运用FFT算法原理计算DFT，如图1.22、图1.23、图1.24所示。
1.3FFT算法特点
（1）原位运算。
从图2.9的FFT运算流图中我们可以看出，这种运算是很有规律可循的。
其中的每级（每列）运算都是由N/2个蝶形运算所构成，如式（1.31）所示。
（式1.31）
其中，m表示第m列迭代，k,j表示数据所在的行数，该式所对应的蝶形运算结构如图1.31所示。
由图1.31的流程图我们可以看出，对于任一个蝶形中的任何两个节点k和j，经过运算后所得的结果作为下一列（级）k,j两个节点的节点变量，同其他节点变量无关，这样就可以采用原位运算方法。
即某一列的N个数据送入存储器经蝶形运算后，结果为另一级（列）数据，而这些数据结果可以以蝶形为单位仍然存储在这同一组存储器中，直到最后输出，中间无需另加存储器。
也就是说，蝶形运算的两个输出值仍放回到蝶形的两个输入所在的存储器中（即将原先的输入值覆盖掉）。
前一级蝶形运算完成后才进行下一级蝶形原位运算，因而整个运算结
构由于采用这种原位运算而大大减少了存储单元的个数，有效地降低了成本。
（2）倒位序规律。
由图2.9我们还可看出，由于采用了原位计算方法，FFT的输出X（k）是按k值的自然顺序排列在存储单元中的，即排列序列为X（0）,X
（1）,…,X（7），而输入的时间序列都不是按照自然顺序排列，而是按x（0）,x（4）,x
（2）,…,x（7），输入到存储单元中。
这种排列初看起来像是“混乱无序”的，但实质是有序的。
如果我们用二进制数来表示0～7个数，设为n（n2n1n0）2，然后再来观察n（n0n1n2）2，如表1.1所示。
十进制数
二进制数
十进制数
000
1
001
100
4
2
010
3
011
110
6
5
101
7
111
由表我们可以看出，只要将
转换成
，即将二进制的最高位与最低位相交换、次高位与次低位相交换……所得的n就是以自然顺序排列的，故通常n的这种排列规律我们称之为倒位序。
这两条性质给我们编程带来了极大的方便。
2．C语言实现FFT算法
2.1算法基本原理
如上所述，FFT算法给编程带来了极大大方便，对比于DFT不仅运算次数大大降低，就连所需的内存空间在运算中也不会增加。
主要算法有两部分，一个是倒二进制排序，还有最重要的蝶形算法。
倒二进制算法
对数组下标进行倒二进制运算，把输入数组按倒二进制排序即可，算法很简单，具体可参见代码
蝶形算法
核心即为如下公式，只需进行几次循环分级进行蝶形算法即可，具体可参见代码。
2.2程序流程图
其中的N表示运算的总级数，即将FFT时域序列分出奇偶序列，直至每个序列点数为1所需的次数。
若有1024个点，则N为10。
3.FFT的Matlab实验
3.1改变序列的长度N对计算结果的影响
x=0.5sin（2*pi*15*t）+2*sin（2*pi*40*t）。
采样频率fs=100Hz，分别绘制N=128、1024点幅频图。
在Matlab中输入如下代码：
>
clf;
fs=100;
N=128;
%采样频率和数据点数
n=0:
N-1;
t=n/fs;
%时间序列
x=0.5*sin（2*pi*15*t）+2*sin（2*pi*40*t）;
%信号
y=fft（x,N）;
%对信号进行快速Fourier变换
mag=abs（y）;
%求得Fourier变换后的振幅
f=n*fs/N;
%频率序列
subplot（2,2,1）,plot（f,mag）;
%绘出随频率变化的振幅
xlabel（'
频率/Hz'
）;
ylabel（'
振幅'
title（'
N=128'
gridon;
subplot（2,2,2）,plot（f（1:
N/2）,mag（1:
N/2））;
%绘出Nyquist频率之前随频率变化的振幅
%对信号采样数据为1024点的处理
N=1024;
%求取Fourier变换的振幅
subplot（2,2,3）,plot（f,mag）;
N=1024'
subplot（2,2,4）
plot（f（1:
x=0.5*sin（2*pi*15*t）+2*sin（2*pi*40*t）。
采样频率fs=100Hz，当N=128、1024点幅频图如下
fs=100Hz，Nyquist频率为fs/2=50Hz。
整个频谱图是以Nyquist频率为对称轴的。
并且可以明显识别出信号中含有两种频率成分：
15Hz和40Hz。
由此可以知道FFT变换数据的对称性。
因此用FFT对信号做谱分析，只需考察0~Nyquist频率范围内的福频特性。
若没有给出采样频率和采样间隔，则分析通常对归一化频率0~1进行。
另外，振幅的大小与所用采样点数有关，采用128点和1024点的相同频率的振幅是有不同的表现值，但在同一幅图中，40Hz与15Hz振动幅值之比均为4：
1，与真实振幅0.5：
2是一致的。
为了与真实振幅对应，需要将变换后结果乘以2除以N。
由上图可得，点数取得越多，由FFT所得的频谱越精确。
4.总结
傅里叶变换在物理学、数论、组合数学、信号处理、概率论、统计学、密码学、声学、光学、海洋学、结构动力学等领域都有着广泛的应用（例如在信号处理中，傅里叶变换的典型用途是将信号分解成幅值分量和频率分量）。
而FFT则实现了傅里叶变换在计算机领域的广泛应用，学习FFT的原理和实现，对后续课程如电力系统的分析等是很有必要的。
本次课程设计，我选择了FFT在计算机上的实现这一个题目，采用C语言编程实现FFT算法，同时用Matlab对C程序进行验证和分析。
在完成课设过程中，不仅弄清楚了FFT的原理，还学到了Matlab相关的知识，C语言编程能力也得到了较大的提高。
附录程序清单
1.1复数的算法
#include"
complex.h"
CComplex:
:
CComplex（void）
{
Realpart=0;
ImagPart=0;
}
CComplex（doubleReal,doubleImag）
Realpart=Real;
ImagPart=Imag;
CComplex（doublereal）
Realpart=real;
~CComplex（void）
/////////////////////////////////////////////////////////////////////////////////
/*重载数学运算符*/
//+重载
CComplexCComplex:
operator+（CComplexOtherComplex）
CComplextemp;
temp.Realpart=Realpart+OtherComplex.Realpart;
temp.ImagPart=ImagPart+OtherComplex.ImagPart;
returntemp;
//-重载
operator-（CComplexOtherComplex）
temp.Realpart=Realpart-OtherComplex.Realpart;
temp.ImagPart=ImagPart-OtherComplex.ImagPart;
//*重载
operator*（CComplexOtherComplex）
{
temp.Realpart=Realpart*OtherComplex.Realpart-ImagPart*OtherComplex.ImagPart;
temp.ImagPart=Realpart*OtherComplex.ImagPart+ImagPart*OtherComplex.Realpart;
1.2FFT运算
#defineFFTPOINT8
//////////////////////////////////////////////////////////////////////////////////////
/*输入时域取样数据*/
voidInPut（double*pData）
doubleIpt;
for（inti=0;
i<
FFTPOINT;
i++）
{
cin>
Ipt;
*（pData+i）=Ipt;
}
/*输出变换后的频域数据，复数形式*/
voidOutput（CComplex*pData）
CComplexOut;
cout<
<
"
FFT变换后，频域数据为"
endl;
for（inti=0;
i++）
{Out=*（pData+i）;
if（（Out.Realpart）>
=0）
{
printf（"
%.4fj+%.4f\n"
Out.ImagPart,Out.Realpart）;
}
else
%.4fj%.4f\n"
///////////////////////////////////////////////////////////////////////////////////////
/*将输入的数据按倒二进制排序*/
voidFFTArry（double*pData1,CComplex*pData2）
inti=0,m=2,t=0;
for（N=0;
m<
=FFTPOINT;
N++）
m*=2;
for（i=0;
inttemp=i;
t=0;
for（intj=N-1;
j>
=0;
j--）
{
if（temp/（int）（pow（2,j）））
t+=pow（2,N-j-1）;
temp=temp%（int）（pow（2,j））;
}
pData2[t].Realpart=*（pData1+i）;
////////////////////////////////////////////////////////////////////////////////////////////////////
externintN=0;
voidmain（）
voidFFT（CComplex*pData）;
doubleFFTDataIn[FFTPOINT]={0.0};
CComplexFFTDataOut[FFTPOINT];
请输入原数据,最大个数为：
FFTPOINT<
InPut（FFTDataIn）;
//输入数据
FFTArry（FFTDataIn,FFTDataOut）;
//将输入数据按到二进制排序
FFT（FFTDataOut）;
//进行FFT变换
Output（FFTDataOut）;
//输出
voidFFT（CComplex*pData）
inti,r,steplen,k,j;
for（i=0;
N;
i++）/*实现N级运算*/
steplen=（int）（pow（2,i））;
/*计算步长*/
for（r=0;
r<
pow（2,i）;
r++）
CComplexW（cos（2*PI*r/pow（2,i+1））,-sin（2*PI*r/pow（2,i+1）））;
/*旋转因子的值*/
for（j=0;
j<
（FFTPOINT/pow（2,i+1））;
j++）/*r相同的值的组的计算*/
{
CComplextemp1,temp2;
k=（int）（r+j*pow（2,i+1））;
/*每一组成员的下标*/
temp1=pData[k]+W*pData[k+steplen];
temp2=pData[k]-W*pData[k+steplen];
pData[k]=temp1;
pData[k+steplen]=temp2;
}
//////////////////////////////////////////////////////////////////////////////////////////////////
升级会员 