fft变换Word格式文档下载.docx
- 文档编号:21409733
- 上传时间:2023-01-30
- 格式:DOCX
- 页数:14
- 大小:302.27KB
fft变换Word格式文档下载.docx
《fft变换Word格式文档下载.docx》由会员分享,可在线阅读,更多相关《fft变换Word格式文档下载.docx(14页珍藏版)》请在冰豆网上搜索。
用数学表达式就是如下:
S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)
式中cos参数为弧度,所以-30度和90度要分别换算成弧度。
我们以256Hz的采样率对这个信号进行采样,总共采样256点。
按照我们上面的分析,Fn=(n-1)*Fs/N,我们可以知道,每两个点之间的间距就是1Hz,第n个点的频率就是n-1。
我们的信号有3个频率:
0Hz、50Hz、75Hz,应该分别在第1个点、第51个点、第76个点上出现峰值,其它各点应该接近0。
实际情况如何呢?
我们来看看FFT的结果的模值如图所示。
图1FFT结果
从图中我们可以看到,在第1点、第51点、和第76点附近有比较大的值。
我们分别将这三个点附近的数据拿上来细看:
1点:
512+0i
2点:
-2.6195E-14-1.4162E-13i
3点:
-2.8586E-14-1.1898E-13i
50点:
-6.2076E-13-2.1713E-12i
51点:
332.55-192i
52点:
-1.6707E-12-1.5241E-12i
75点:
-2.2199E-13-1.0076E-12i
76点:
3.4315E-12+192i
77点:
-3.0263E-14+7.5609E-13i
很明显,1点、51点、76点的值都比较大,它附近的点值都很小,可以认为是0,即在那些频率点上的信号幅度为0。
接着,我们来计算各点的幅度值。
分别计算这三个点的模值,
结果如下:
512
384
192
按照公式,可以计算出直流分量为:
512/N=512/256=2;
50Hz信号的幅度为:
384/(N/2)=384/(256/2)=3;
75Hz信号的幅度为192/(N/2)=192/(256/2)=1.5。
可见,从频谱分析出来的幅度是正确的。
然后再来计算相位信息。
直流信号没有相位可言,不用管它。
先计算50Hz信号的相位,atan2(-192,332.55)=-0.5236,结果是弧度,换算为角度就是180*(-0.5236)/pi=-30.0001。
再计算75Hz信号的相位,atan2(192,3.4315E-12)=1.5708弧度,换算成角度就是180*1.5708/pi=90.0002。
可见,相位也是对的。
根据FFT结果以及上面的分析计算,我们就可以写出信号的表达式了,它就是我们开始提供的信号。
总结:
假设采样频率为Fs,采样点数为N,做FFT之后,某一点n(n从1开始)表示的频率为:
Fn=(n-1)*Fs/N;
该点的模值除以N/2就是对应该频率下的信号的幅度(对于直流信号是除以N);
该点的相位即是对应该频率下的信号的相位。
相位的计算可用函数atan2(b,a)计算。
atan2(b,a)是求坐标为(a,b)点的角度值,范围从-pi到pi。
要精确到xHz,则需要采样长度为1/x秒的信号,并做FFT。
要提高频率分辨率,就需要增加采样点数,这在一些实际的应用中是不现实的,需要在较短的时间内完成分析。
解决这个问题的方法有频率细分法,比较简单的方法是采样比较短时间的信号,然后在后面补充一定数量的0,使其长度达到需要的点数,再做FFT,这在一定程度上能够提高频率分辨力。
具体的频率细分法可参考相关文献。
[附录:
本测试数据使用的matlab程序]
closeall;
%先关闭所有图片
Adc=2;
%直流分量幅度
A1=3;
%频率F1信号的幅度
A2=1.5;
%频率F2信号的幅度
F1=50;
%信号1频率(Hz)
F2=75;
%信号2频率(Hz)
Fs=256;
%采样频率(Hz)
P1=-30;
%信号1相位(度)
P2=90;
%信号相位(度)
N=256;
%采样点数
t=[0:
1/Fs:
N/Fs];
%采样时刻
%信号
S=Adc+A1*cos(2*pi*F1*t+pi*P1/180)+A2*cos(2*pi*F2*t+pi*P2/180);
%显示原始信号
plot(S);
title('
原始信号'
);
figure;
Y=fft(S,N);
%做FFT变换
Ayy=(abs(Y));
%取模
plot(Ayy(1:
N));
%显示原始的FFT模值结果
FFT模值'
Ayy=Ayy/(N/2);
%换算成实际的幅度
Ayy
(1)=Ayy
(1)/2;
F=([1:
N]-1)*Fs/N;
%换算成实际的频率值
plot(F(1:
N/2),Ayy(1:
N/2));
%显示换算后的FFT模值结果
幅度-频率曲线图'
Pyy=[1:
N/2];
fori="
1:
N/2"
Pyy(i)=phase(Y(i));
%计算相位
Pyy(i)=Pyy(i)*180/pi;
%换算为角度
end;
N/2),Pyy(1:
%显示相位图
相位-频率曲线图'
看完这个你就明白谐波分析了。
FFT算法原理和实现
2011-05-1522:
16:
45|
分类:
默认分类|举报|字号
订阅
在数字信号处理中常常需要用到离散傅立叶变换(DFT),以获取信号的频域特征。
尽管传统的DFT算法能够获取信号频域特征,但是算法计算量大,耗时长,不利于计算机实时对信号进行处理。
因此至DFT被发现以来,在很长的一段时间内都不能被应用到实际的工程项目中,直到一种快速的离散傅立叶计算方法——FFT,被发现,离散是傅立叶变换才在实际的工程中得到广泛应用。
需要强调的是,FFT并不是一种新的频域特征获取方式,而是DFT的一种快速实现算法。
本文就FFT的原理以及具体实现过程进行详尽讲解。
DFT计算公式
本文不加推导地直接给出DFT的计算公式:
其中x(n)表示输入的离散数字信号序列,WN为旋转因子,X(k)一组N点组成的频率成分的相对幅度。
一般情况下,假设x(n)来自于低通采样,采样频率为fs,那么X(k)表示了从-fs/2率开始,频率间隔为fs/N,到fs/2-fs/N截至的N个频率点的相对幅度。
因为DFT计算得到的一组离散频率幅度只实际上是在频率轴上从成周期变化的,即X(k+N)=X(k)。
因此任意取N个点均可以表示DFT的计算效果,负频率成分比较抽象,难于理解,根据X(k)的周期特性,于是我们又可以认为X(k)表示了从零频率开始,频率间隔为fs/N,到fs-fs/N截至的N个频率点的相对幅度。
N点DFT的计算量
根据
(1)式给出的DFT计算公式,我们可以知道每计算一个频率点X(k)均需要进行N次复数乘法和N-1次复数加法,计算N各点的X(k)共需要N^2次复数乘法和N*(N-1)次复数加法。
当x(n)为实数的情况下,计算N点的DFT需要2*N^2次实数乘法,2*N*(N-1)次实数加法。
旋转因子WN的特性
1.WN的对称性
2.WN的周期性
3.WN的可约性
根据以上这些性质,我们可以得到式(5)的一些列有用结果
基-2FFT算法推导
假设采样序列点数为N=2^L,L为整数,如果不满足这个条件可以人为地添加若干个0以使采样序列点数满足这一要求。
首先我们将序列x(n)按照奇偶分为两组如下:
于是根据DFT计算公式
(1)有:
至此,我们将一个N点的DFT转化为了式(7)的形式,此时k的取值为0到N-1,现在分为两段来讨论,当k为0~N/2-1的时候,因为x1(r),x2(r)为N/2点的序列,因此,此时式(7)可以写为:
而当k取值为N/2~N-1时,k用k’+N/2取代,k’取值为0~N/2-1。
对式(7)化简可得:
综合以上推导我们可以得到如下结论:
一个N点的DFT变换过程可以用两个N/2点的DFT变换过程来表示,其具体公式如式(10)所示DFT快速算法的迭代公式:
上式中X'
(k’)为偶数项分支的离散傅立叶变换,X'
'
(k’’)为奇数项分支的离散傅立叶变换。
式(10)的计算过程可以用图1的蝶形算法流图直观地表示出来。
图1时间抽取法蝶形运算流图
在图1中,输入为两个N/2点的DFT输出为一个N点的DFT结果,输入输出点数一致。
运用这种表示方法,8点的DFT可以用图2来表示:
图28点DFT的4点分解
根据公式(10),一个N点的DFT可以由两个N/2点的DFT运算构成,再结合图1的蝶形信号流图可以得到图2的8点DFT的第一次分解。
该分解可以用以下几个步骤来描述:
1.将N点的输入序列按奇偶分为2组分别为N/2点的序列
2.分别对1中的每组序列进行DFT变换得到两组点数为N/2的DFT变换值X1和X2
3.按照蝶形信号流图将2的结果组合为一个N点的DFT变换结果
根据式(10)我们可以对图2中的4点DFT进一步分解,得到图3的结果,分解步骤和前面一致。
图38点DFT的2点分解
最后对2点DFT进一步分解得到最终的8点FFT信号计算流图:
图48点DFT的全分解
从图2到图4的过程中关于旋转系数的变化规律需要说明一下。
看起来似乎向前推一级,在奇数分组部分的旋转系数因子增量似乎就要变大,其实不是这样。
事实上奇数分组部分的旋转因子指数每次增量固定为1,只是因为每向前推进一次,该分组序列的数据个数变少了,为了统一使用以原数据N为基的旋转因子就进行了变换导致的。
每一次分组奇数部分的系数WN,这里的N均为本次分组前的序列点数。
以上边的8点DFT为例,第一次分组N=8,第二次分组N为4,为了统一根据式(4)进行了变换将N变为了8,但指数相应的需要乘以2。
N点基-2FFT算法的计算量
从图4可以看到N点DFT的FFT变换可以转为log2(N)级级联的蝶形运算,每一级均包含有N/2次蝶形计算。
而每一个蝶形运算包含了1次复数乘法,2次复数加法。
因此N点FFT计算的总计算量为:
复数乘法——N/2×
log2(N)
复数加法——N×
log2(N)。
假设被采样的序列为实数序列,那么也只有第一级的计算为实数与复数的混合计算,经过一次迭代后来的计算均变为复数计算,在这一点上和直接的DFT计算不一致。
因此对于输入序列是复数还是实数对FFT算法的效率影响较小。
一次复数乘法包含了4次实数乘法,2次实数加法,一次复数加法包含了2次复数加法。
因此对于N点的FFT计算需要总共的实数乘法数量为:
2×
N×
log2(N);
总的复数加法次数为:
2xNxlog2(N)。
N点基-2FFT算法的实现方法
从图4我们可以总结出对于点数为N=2^L的DFT快速计算方法的流程:
1.对于输入数据序列进行倒位序变换。
该变换的目的是使输出能够得到X(0)~X(N-1)的顺序序列,同样以8点DFT为例,该变换将顺序输入序列x(0)~x(7)变为如图4的x(0),x(4),x
(2),x(6),x
(1),x(5),x(3),x(7)序列。
其实现方法是:
假设顺序输入序列一次村在A(0)~A(N-1)的数组元素中,首先我们将数组下标进行二进制化(例:
对于点数为8的序列只需要LOG2(8)=3位二进制序列表示,序号6就表示为110)。
二进制化以后就是将二进制序列进行倒位,倒位的过程就是将原序列从右到左书写一次构成新的序列,例如序号为6的二进制表示为110,倒位后变为了011,即使十进制的3。
第三步就是将倒位前和倒位后的序号对应的数据互换。
依然以序号6为例,其互换过程如下:
temp=A(6);
A(6)=A(3);
A(3)=A(6);
实际上考虑到执行效率,如果对于每一次输入的数据都需要这个处理过程是非常浪费时间的。
我们可以采用指向指针的指针来实现该过程,或者是采用指针数组来实现该过程。
2.蝶形运算的循环结构。
从图4中我们可以看到对于点数为N=2^L的fft运算,可以分解为L阶蝶形图级联,每一阶蝶形图内又分为M个蝶形组,每个蝶形组内包含K个蝶形。
根据这一点我们就可以构造三阶循环来实现蝶形运算。
编程过程需要注意旋转因子与蝶形阶数和蝶形分组内的蝶形个数存在关联。
3.浮点到定点转换需要注意的关键问题
上边的分析都是基于浮点运算来得到的结论,事实上大多数嵌入式系统对浮点运算支持甚微,因此在嵌入式系统中进行离散傅里叶变换一般都应该采用定点方式。
对于简单的DFT运算从浮点到定点显得非常容易。
根据式
(1),假设输入x(n)是经过AD采样的数字序列,AD位数为12位,则输入信号范围为0~4096。
为了进行定点运算我们将旋转因子实部虚部同时扩大2^12倍,取整数部分代表旋转因子。
之后,我们可以按照
(1)式计算,得到的结果与原结果成比例关系,新的结果比原结果的2^12倍。
但是,对于使用蝶形运算的fft我们不能采用这种简单的放大旋转因子转为整数计算的方式。
因为fft是一个非对称迭代过程,假设我们对旋转因子进行了放大,根据蝶形流图我们可以发现其最终的结果是,不同的输入被放大了不同的倍数,对于第一个输入x(0)永远也不会放大。
举一个更加形象的例子,还是以图4为例。
从图中可以看出右侧的X(0)可以直接用下式表示:
从上式我们可以看到不同输入项所乘的旋转因子个数(注意这里是个数,就算是wn^0,也被考虑进去了,因为在没有放大时wn^0等于1,放大后所有旋转因子指数模均不为1,因此需要考虑)。
这就导致输入不平衡,运算结果不正确。
经查阅相关资料,比较妥善的做法是,首先对所有旋转因子都放大2^Q倍,Q必须要大于等于L,以保证不同旋转因子的差异化。
旋转因子放大,为了保证其模为1,在每一次蝶形运算的乘积运算中我们需要将结果右移Q位来抵消这个放大,从而得到正确的结果。
之所以采用放大倍数必须是2的整数次幂的原因也在于此,我们之后可以通过简单的右移位运算将之前的放大抵消,而右移位又代替了除法运算,大大节省了时间。
4.计算过程中的溢出问题
最后需要注意的一个问题就是计算过程中的溢出问题。
在实际应用中,AD虽然有12位的位宽,但是采样得到的信号可能较小,例如可能在0~8之间波动,也就是说实际可能只有3位的情况。
这种情况下为了在计算过程中不丢失信息,一般都需要先将输入数据左移P位进行放大处理,数据放大可能会导致溢出,从而使计算错误,而溢出的极限情况是这样:
假设我们数据位宽为D位(不包括符号位),AD采样位数B位,数字放大倍数P位,旋转因此放大倍数Q位,FFT级联运算带来的最大累加倍数L位。
我们得到:
假设AD位宽12,数据位宽32,符号位1位,因此有效位宽31位,采样点数N,那么我们可以得到log2(N)+P+Q<
=19,假设点数128,又Q>
=L可以得到放大倍数P<
=5。
FFT代码示例
根据以上的分析,我整理了一份128点的FFT代码如下,该代码在keil中仿真运行,未发现问题。
#defineN128
#definePOWER6//该值代表了输入数据首先被放大的倍数,放大倍数为2^POWER
#defineP_COEF8//该值代表了旋转因子被放大的倍数,放大倍数为2^POWER
#if(N==4)
#defineL2//L的定义满足L=log2(N)
#elif(N==8)
#defineL3//L的定义满足L=log2(N)
#elif(N==16)
#defineL4//L的定义满足L=log2(N)
#elif(N==32)
#defineL5//L的定义满足L=log2(N)
#elif(N==64)
#defineL6//L的定义满足L=log2(N)
#elif(N==128)
#defineL7//L的定义满足L=log2(N)
#endif
//AD采样位数为12位,本可以采用s16x[N]定义数据空间,但是为了节省存储空间,fft结果也将存储在该变量空间。
由于受
//fft影响变量需要重新定义,定义的方式及具体原因如下:
//fft过程会乘以较大旋转因子,因此需要32位定义
//fft过程会产生负数,因此需要有符号定义
//fft过程会出现复数,因此需要定义二维数组,[][0]存放实部,[][1]存放虚部
s32x[N][2]={};
//定义*p[N]是为了在第一次指针初始化以后,该数组指针按照位倒序指向数据变量x
//p[i][0]代表了指向数据的实部,p[i][1]代表了指向数据的虚部
s32*p[N];
//定义旋转因子矩阵
//旋转因子矩阵存储了wn^0,wn^1,wn^2...wn^(N/2-1)这N/2个复数旋转因子
s16w[N>
>
1][2]={256,0,256,-13,255,-25,253,-38,251,-50,248,-62,245,-74,241,-86,237,-98,231,-109,226,
-121,220,-132,213,-142,206,-152,198,-162,190,-172,181,-181,172,-190,162,-198,152,
-206,142,-213,132,-220,121,-226,109,-231,98,-237,86,-241,74,-245,62,-248,50,-251,38,
-253,25,-255,13,-256,0,-256,-13,-256,-25,-255,-38,-253,-50,-251,-62,-248,-74,-245,-86,
-241,-98,-237,-109,-231,-121,-226,-132,-220,-142,-213,-152,-206,-162,-198,-172,-190,-181,
-181,-190,-172,-198,-162,-206,-152,-213,-142,-220,-132,-226,-121,-231,-109,-237,-98,-241,
-86,-245,-74,-248,-62,-251,-50,-253,-38,-255,-25,-256,-13};
voidinit_pointer(void)
{
unsignedchari=0;
unsignedcharj=0;
unsignedchark=0;
for(i=0;
i<
N;
i++)
j=0;
for(k=0;
k<
L;
k++)
{j|=(((i>
k)&
0x01)<
<
(L-k-1));
}
p[i]=&
x[j][0];
}
/**description:
基2fft主函数,该函数将借助指针数组p对全局变量数组x进行fft计算
*并将结果存储在数组x中
*globalvar:
rw->
x;
r->
p,w。
(r表示读,w表示写,rw表示读写)*/
voidfft2(void)
u8i;
//i用于表示蝶形图级联的阶数
u8j;
//表示蝶形分组起始点序列,蝶形分组跨度为2^i
u8k;
//k表示蝶形组内偶数分支序列点号
u8gp_distance=1;
//蝶形分组跨度
u8temp;
//temp用于记录当前组间距离的一半,同时也表示了同一碟形两分支间的距离
u8gp_hf=0;
//记录当前组的中间点序号
u8delta=N;
//旋转因子下标增量,本来下标初始值应该为N/2,但是。
。
s16*pw=&
(w[0][0]);
inttp_result[2];
//用于临时存放旋转因子和奇数分组的乘积
//输入信号序列放大
{x[i][0]<
=POWER;
x[i][1]<
}for(i=0;
i++){temp=gp_distance;
gp_distance<
=1;
for(j=0;
j<
j+=gp_d
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- fft 变换