1、6.3 频谱分析6.3.1 FFT 算法傅立叶变换是实现信号从时域变换到频域的工具，其定义为：(6.3.1)()()jtXxed其实质是把信号 分解为一系列的正弦或余弦。根据这一思想，傅立叶逆变换也成立，t频域信号 也可通过逆傅立叶变换回到时域形式：()(6.3.2)1()2jtxtXed对于采样信号 ，采用离散傅立叶变换(DFT) ，将其变换到频域。xn(6.3.3)10)()(NnkWk式中， 。根据式(6.3.3) ，要完成一次 DFT，需要进行 N2 次复数乘法和 N(N-1)次kNjke2复数加法。如果采样点数 N 较大的话，DFT 的计算量相当大。为此，J.W.Cooley 和T.
2、W.Tukey 于 1965 年提出实现 DFT 的快速算法，即 FFT。其实质是将信号序列分级分组，通过短序列的 DFT，实现整个序列的 DFT，这样可将 DFT 的计算量减小至一到二个数量级。式(6.3.3) 中的 具有对称性和周期性，如式(6.3.4) 和式( 6.3.5)所示：knNW(6.3.4)knN(6.3.5)(快速傅立叶变换 FFT 正是运用 的对称和周期性，将长序列变短，推导如下：knN设采样信号序列 ， ，采样点为 N，假设 ，M 为正整)(nx1,20 2数，因此可将信号序列 按奇偶分成两组，令偶数组 ，及奇数组 ，rn21rn，式(6.3.3)可变换为：12/,0Nr
3、 12/0)12(/0()()(NrkrNrkWxxkX(6.3.6)rkNNrkNrrkWxWx2/12/012/02/ )()( 令 ,12/02/)()(rrkNkA12/0)12(/()(rkrB即 、 分别是偶数组、奇数组的傅立叶变换，结合式(6.3.4)和式(6.3.5) ，式(6.3.6)变成：(6.3.7)()(kBWAkXN 12/,10N(6.3.8)2/k这就是 FFT 的核心公司。以 为例，这个过程如图 6.3.1 所示。8A(0)(1)A(2)(3)B(0)(1)B(2)(4) X(0)(1)X(2)(3)X(4)(5)X(6)(7)18W283808 11图 6.3
4、.1 FFT 算法分组实现原理图依此类推， 分组后得到奇数组和偶数组，可继续按其位置序号进行类似的分组，)(nx直到得到最简单的蝶形单元，每组只有两个元素，即不能再分组，仍以 为例，实现8N过程如图 6.3.2 所示。 A(0) 1(1)(2)A(3)B(0)(1)B(2)(3) X(0)(1)(2)X(3)(4)(5)X(6)(7)18W283808 1C(0)(1)D(0)(1)E(0)(1)F(0)(1)x(0)(4)x(2)(6)x(1)(5)x(3)(7) 08W08W2828 111111图 6.3.2 8 点 FFT 算法实现原理图FFT 的具体步骤为：首先,按照“码位倒置”的规
5、律对信号进行重新排序。因为 FFT 算法需要对信号序列不断进行奇偶分组，信号重排正好符合“码位倒置”的规律。如 N8，序号 用二进制数可n用 表示，第一次分组时， 为偶数组， 为奇数组；第二次分组则按)(012n0n10n划分，同样第三次分组则按 划分。N 8 时只需要三次分组，最终的排列顺序正好为2原来码位的倒置形式即 ，整个信号序列重排的过程如图 6.3.3 所示。图中，I 代表)(10顺码，J 代表倒码，如果 I=J，则不需调换；如果 I 不等于 J，则需要调换；若 I 小于 J 时进行调换，而 I 大于 J 时已经完成调换，不需再调换了。开始结束输入 NN V 2 = N / 2M N
6、 = N - 1I = J ?T = A ( J )A ( J ) = A ( I )A ( I ) = TK = N V 2K = J ?J = J + KI = I + 1I M N 1 ?J = J - KK = K / 2YNYNYN图 6.3.3 序列重排流程图信号序列经过序列重排之后，通过“分级”和“分组”，蝶形运算，完成 FFT。按照FFT 的实现原理，整个运算分成 M 级( )，每级里分成 组， 代表级数，2logNi2/每组组成一个蝶形单元，根据式(6.3.7)和(6.3.8)完成蝶形运算。如图 6.3.2，图中 N=8，分成 3 级，第一级中分成 4 组，每组组成最简单的蝶
7、形运算单元；在第二级中，有 2 组；第三级则为 1 组。“级”和“组”组成快速算法中的两个循环，FFT 算法实现流程图如图 6.3.4所示。图中，A(I)为采样数组重排之后的序列， M 为 FFT 算法要求进行的级数，LE 为每组的长度。开始结束输入 N , M , A ( I )L = 1E2L E 1 = L E / 2U = 1 1jLWeJ = 1I = JI P = I + L E 1T = A ( I P ) * UA ( I P ) = A ( I ) - TA ( I ) = A ( I ) + TI = I + L EI N ?U + U * WJ = J + 1J L E
8、1 ?L = L + 1L M ?码位倒置YYYNNN图 6.3.4 FFT 算法实现流程图6.3.2 功率谱计算在时域内，对于任意实信号 ，其平均功率为：)(tx(6.3.9)2/1lim()Txtd式中， 为信号的瞬时功率，在 T 有限情况下，式(6.3.9)的极限存在，其对应的谱为2t功率谱。若 收敛的话，则 代表信号 的总能量。由巴什瓦定理，2()xtd 2()xtd )(tx信号的功率和频谱之间为：(6.3.10)dXdtx22)(1)(式中， 为信号 幅值谱的平方，通常称为功率谱或能量谱，对于具有连续频谱2tx的信号来说，称为功率谱密度（能量谱密度），表示频率尺度上单位间隔的功率（
9、能量）。式(6.3.10)表明信号的功率（能量）可通过其功率谱密度在频率轴上的积分求得。对于离散信号序列来说，功率谱的经典估计方法主要有两种：相关估计法和周期图谱法。在和 完善 FFT 算法之前，相关估计法应用较多。但在此之后，因为 FFT 快ColeyTurk速算法大大减少傅立叶变换所需的计算量，采用 FFT 的周期图谱法 ，比相关估计法实现起来更为简单。对于采样点为 N 的采样序列，周期图谱法功率谱估计公式为：(6.3.11)2)(1)(XSx根据式(6.3.11)，只要通过 FFT 得到信号的幅值谱，就可估计出信号的功率谱。6.3.3 频谱校正以 FFT 为基础的经典功率谱估计方法，应用
10、非常广泛。但是，也存在一些缺陷：(1)数字信号处理的对象通常是具有一定点数的离散信号序列，实质上对真实信号加上了一个矩形窗，这会造成频谱分析的吉普森现象，即信号的功率谱向矩形窗频谱的旁瓣泄漏，造成信号功率谱估计的偏差；(2)经典谱估计的频率分辨率不高，它与信号采样的时间长度成反比，采样的时间越长，分辨率就越高，但这与数字信号处理系统实时运行的要求相矛盾。频谱泄漏是因为加于信号的矩形窗造成，实际应用可通过选择合适的窗来平滑信号的功率谱曲线，改善这种频谱泄漏的情况，如：汉宁窗(Hanning)、哈明窗(Hamming) 等。不过在应用功率谱通过能量最大原则估计信号频率的算法中，频谱泄漏对计算结果影
11、响不是很大，影响较大的是非整周期采样所造成的频率计算误差，最大误差为频率分辨率 的 0.5f倍。由于 ，应用中可通过减小 或提高 来降低这个最大误差，但这不是sNf/sfsN根本的解决方法而且受数字信号处理系统要求的限制，因为：减小 就会降低系统的测量sf范围，提高 会增加系统的计算量。因此，需要对信号频谱进行校正，来提高频率估计s的精度。频谱校正有多种方法，如：频谱细化、三点卷积幅值校正、重心校正和峰值搜寻等，通过这些校正方法都可估计出较准确的频率、幅值和相位。对频率估计来说，影响最大的是频率分辨率，对信号幅值相位要求不高，所以只需对频率进行校正，可以选用计算简单、物理意义明确的重心校正法。
12、重心法频谱校正原理是：挑出信号能量最大的谱线及其相邻的次高谱线，通过重心法规求取它们的重心谱线，认为这条谱线就是主瓣中心谱线。校正的公式为：(6.3.12)()1()(1)nYnYn令 ()由频率的一般形式 ，可得Nfns/(6.3.13)()/sfnf式中， 、 代表最高谱线和次高谱线的序号， 代表幅值， 、 则就是频谱校正后1Ynf的主瓣重心的序号以及对应的频率。6.3.4 源程序及注释/预处理程序void pp_sgnl() int i,j; *删掉 J*double temp_sum;float temp_aver;/去平均值开始temp_sum=0;for(i=0;imax)max=tempfloat;k=i; /开始频谱校正if(k!=0)|(k!=(NN/2-1)a0=x0k-x0k-1;a1=x0k-x0k+1;Nf=k;if(a0a1)Nf=Nf+x0k+1/(x0k+x0k+1); /频率校正max=(x0k+1*(k+1)+x0k*k)/Nf; /幅值校正tempfre=Nf*fs/NN;/校正后的频率 amp_now=max*Amp_unit; /校正后的幅值，其中，Amp_unit 是 ADC 输入的最大电压除以 1024。
