信号处理快速傅里叶变换(FFT)的学习

FFT是离散傅立叶变换的快速算法，可以将一个信号变换到频域。有些信号在时域上是很难看出什么特征的，但是如果变换到频域之后，就很容易看出特征了。这就是很多信号分析采用FFT变换的原因。另外，FFT可以将一个信号的频谱提取出来，这在频谱分析方面也是经常用的。采样频率要大于信号频率的两倍。为了方便进行FFT运算，通常N取2的整数次方。

一、DFT

1.1 原理

定义（来自百科）：离散傅里叶变换（Discrete Fourier Transform，缩写为DFT），是傅里叶变换在时域和频域上都呈离散的形式，将信号的时域采样变换为其DTFT的频域采样。在形式上，变换两端（时域和频域上）的序列是有限长的，而实际上这两组序列都应当被认为是离散周期信号的主值序列。即使对有限长的离散信号作DFT，也应当将其看作其周期延拓的变换。在实际应用中通常采用快速傅里叶变换计算DFT。

DFT的公式是：设x(n)为M点的有限长序列，即在0≤n≤M-1内有值，则定义x(n)的N点（N≥M。当N＞M时，补N-M个零值点），离散傅里叶变化定义为

其中，为旋转因子（为方便编辑后续记作W(N, nk)，特此说明），其计算公式为，具有以下性质：

①周期性：W(N, nk) = W(N, (n+rN)k) = W(N, n(k+rN))，其中r问整数。

②共轭对称性：(W(N, nk))* = W(N, -nk)。

③可约性：W(N, nk) = W(mN, mnk)，W(N, nk) = W(N/m, nk/m)，其中m为整数，N/m为整数。

④特殊值：W(N, N/2) = -1；  W(N, (k+N/2)) = -W(N, k)；  

                  W(N, (N-k)n) = W(N, (N-n)k) = W(N, -nk)。

所以，在计算旋转因子的过程中可以适当的使用特殊值来提高运算的效率。

在计算DFT时，如果数据的点数够计算的点数，则截取计算点数长的数据进行DFT运算，否则将数据点个数补0至计算点个数，然后进行计算，举例如下（举例只是为了说明问题，没有按照编程语言的书写格式）。

比如：数据点数组为 dataArray = {1， 2， 3， 4，5}，可以看出数据长度为5。

如果要求做4点DFT运算，则只需截取前4个数作为运算数组进行运算即可，即为{1， 2， 3， 4}；

如果要求做8点DFT运算，则需在原数组后补三个0，使长度为8后再进行计算，计算数组为{1， 2， 3， 4，5，0，0，0}。

1.2 流程图

二、FFT

2.1 原理

设序列x(n)点数为N=2^L，L为整数。将N=2^L，先按照n的奇偶分为两组，

其中r = 0, 1, ..., N/2-1，x(2r) = x1(r)，x(2r-1) = x2(r)；

则可将DFT化为

由上式可以看出一个N点的DFT可以分为两个N/2点的DFT，按照上式右组合成N点DFT。但是这里的x1(r)、x2(r)以及X1(k)、X2(k)都是N/2点的序列，即r,k满足r,k=0, 1, ..., N/2-1。而X(k)却有N点，上式计算的只是X(k)的前半项数的结果，因此还需要计算后半项的值。

将①②③代入（式1），就可将X(k)表达为前后两部分：

前半部分，X(k)，当k=0, 1, ..., N/2-1

后半部分，X(k)，当k=N/2, ..., N-1

2.2 流程图

①FFT运算总流程图，包括

②根据（1）中推导的内容，FFT的流程图可化为

蝶形运算中又三层循环

第一层（最外层），完成M次迭代过程，即算出A0(k), A1(k), ..., Am(k)，其中k=0, 1, ..., N；A2(k)为蝶形运算第2级的结果，如A0(k)=x(k), Am(k)=X(k)；

第二层（中间层），完成旋转因子的变化，步进为2^L；

第三层（最里层），完成相同旋转因子的蝶形运算。

三、DFT与FFT

1.DFT的运算量

复数乘法次数为N*N，复数加法次数为N(N-1)，若N>>1，则这两者都近似为N*N，它随N增大为急速增大。改进途径：①利用旋转因子性质减小计算量；②由于运算量和N*N成正比，因而可将N点DFT分解成小点数的DFT，以减少运算量（点数越小，计算量越小）；③改为用FFT计算，复数乘法次数为(N/2)*log(2,N)。

2.FFT可分为按照时间抽选的基-2算法（库利-图基算法DIT-FFT）和按频率抽选的基-2算法（桑德-图基算法DIF-FFT）。

3.FFT计算原理

FFT的计算要求点数必须为2的整数次幂，如果点数不够用0补齐。例如计算{2，3，5，8，4}的8点FFT，需要补3个0后进行计算，如果计算该数组的5点FFT，则先计算8点FFT后截取前5个值即可（不提倡）。

四、蝴蝶运算

4.1 基本流程

蝶形运算，符号表示如下图所示：

所以，FFT的蝶形运算图表示为（8点，2^3 = 8，运算级数最大为L=3）

在蝶形运算中变化规律由W(N, p)推导，其中N为FFT计算点数，J为下角标的值

        L = 1时，W(N, p) = W(N, J) = W(2^L, J)，其中J = 0;

        L = 2时，W(N, p) = W(N, J) = W(2^L, J)，其中J = 0, 1;

        L = 3时，W(N, p) = W(N, J) = W(2^L, J)，其中J = 0, 1, 2, 3;

所以，W(N, p) = W(2^L, J)，其中J = 0, 1, ..., 2^(L-1)-1

又因为2^L = 2^M*2^(L-M) = N*2^(L-M)，这里N为2的整数次幂，即N=2^M，

W(N, p) = W(2^L, J) = W(N*2^(L-M), J) = W(N, J*2^(M-L))

所以，p = J*2^(M-L)，此处J = 0, 1, ..., 2^(L-1)-1，当J遍历结束但计算点数不够N时，J=J+2^L，后继续遍历，直到计算点数为N时不再循环。

4.2 举例：N=8点的FFT计算

当L=2时，J = 0, 1两个值，因此p = J*2^(M-L) = 0, 2两个值，即旋转因子有两个值W(8, 0)和W(8, 2)，计算中两行之间的距离B = 2^(L-1)=2^(2-1)=2，

代入J=0, 1可求得X(0)、X(0+2)和X(1)、X(1+2)，即可求出第2级蝶形运算的X(0), X(1), X(2), X(3)，也就是求出一半，此时J加步进2^L=4，即J=J+2^L=4, 5，

再代入J=4, 5可求出X(4), X(4+2)和X(5), X(5+2)，即可求出第2级蝶形运算的X(4), X(5), X(6), X(7)，已经全部求出，J循环结束。

4.3 二进制倒序说明

前面数的排列顺序是进行二进制倒序后的排序。二进制倒序是指将某数转化为二进制表示，将最高位看做最低位、次高位看做次低位...以此类推，计算后的纸进行排序。例如3的二进制表示为011b（3=0*2^2+1*2+1），二进制倒序值为6（011b，最高位看做最低位...即6=0+1*2+1*2^2）

即0到7的二进制排序是：

       0,                4,                 2,                6,                 1,                5,                3,                7

000→000    001→100    010→010    011→110    100→001    101→101    110→011    111→111

五、从振幅或快速傅立叶变换到dB

从振幅或快速傅立叶变换到dB是一个涉及信号处理和单位转换的问题。

函数f(t)为一元连续函数，其傅里叶变换定义为：

其中，i为虚数单位。欧拉公式：

5.1 相关名词

1.振幅（Amplitude）：振幅是指信号的最大偏离值，表示信号的强度或能量大小。在信号处理中，振幅常用于描述波形的高低或振动的幅度。

2.dB（分贝）：分贝是一种用于表示信号强度或功率比例的对数单位。在信号处理中，分贝常用于描述信号的增益或衰减程度。分贝的计算公式为：dB = 10 * log10(P1/P0)，其中P1为信号的功率，P0为参考功率。当保持输入信号的幅度不变，改变频率使输出信号降至最大值的0.707倍，即用频响特性来表述即为-3dB点处即为截止频率，它是用来说明频率特性指标的一个特殊频率。

利用系统函数的模来表示电路的放大倍数，由于20lgA（ω）=-3dB，解得

A（ω）=10^-0.15=0.707945784≈1/√2，又因为A（ω）=|H（jω）|，则|H（jω）|^2=1/2

在高频端和低频端各有一个截止频率，分别称为上截止频率和下截止频率。两个截止频率之间的频率范围称为通频带。

5.2 物理意义理解

假设采样频率为Fs，信号频率F，采样点数为N。那么FFT之后结果就是一个为N点的复数。每一个点就对应着一个频率点。

1.幅值：这个点的模值，就是该频率值下的幅度特性。具体跟原始信号的幅度有什么关系呢？假设原始信号的峰值为A，那么FFT的结果的每个点（除了第一个点直流分量之外）的模值就是A的N/2倍。而第一个点就是直流分量，它的模值就是直流分量的N倍。

2.相位：而每个点的相位呢，就是在该频率下的信号的相位。

3.分辨率：第一个点表示直流分量（即0Hz），而最后一个点N的再下一个点（实际上这个点是不存在的，这里是假设的第N+1个点，也可以看做是将第一个点分做两半分，另一半移到最后）则表示采样频率Fs，这中间被N-1个点平均分成N等份，每个点的频率依次增加。例如某点n所表示的频率为：Fn=(n-1)*Fs/N。由上面的公式可以看出，Fn所能分辨到频率为为Fs/N，如果采样频率Fs为1024Hz，采样点数为1024点，则可以分辨到1Hz。1024Hz的采样率采样1024点，刚好是1秒，也就是说，采样1秒时间的信号并做FFT，则结果可以分析到1Hz，如果采样2秒时间的信号并做FFT，则结果可以分析到0.5Hz。如果要提高频率分辨力，则必须增加采样点数，也即采样时间。频率分辨率和采样时间是倒数关系。

4.表达式：假设FFT之后某点n用复数a+bi表示，那么这个复数的模就是An=根号a*a+b*b，相位就是Pn=atan2(b,a)。根据以上的结果，就可以计算出n点（n≠1，且n<=N/2）对应的信号的表达式为：An/(N/2)*cos(2*pi*Fn*t+Pn)，即2*An/N*cos(2*pi*Fn*t+Pn)。对于n=1点的信号，是直流分量，幅度即为A1/N。

由于FFT结果的对称性，通常我们只使用前半部分的结果，即小于采样频率一半的结果。

5.3 示例

假设我们有一个信号，它含有2V的直流分量，频率为50Hz、相位为-30度、幅度为3V的交流信号，以及一个频率为75Hz、相位为90度、幅度为1.5V的交流信号。用数学表达式就是如下（式中cos参数为弧度，所以-30度和90度要分别换算成弧度）：

S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)

我们以256Hz的采样率对这个信号进行采样，总共采样256点。按照我们上面的分析，Fn=(n-1)*Fs/N，我们可以知道，每两个点之间的间距就是1Hz，第n个点的频率就是n-1。我们的信号有3个频率：0Hz、50Hz、75Hz，应该分别在第1个点、第50个点、第76个点上出现峰值，其它各点应该接近0。

	幅值	相位（弧度）	相位（角度）
直流分量	512/N =512/256 =2	\	\
50Hz	384/(N/2) =384/(256/2) =3	atan2(-192, 332.55) =-0.5236	(-0.5236)/pi =-30.0001
75Hz	192/(N/2) =192/(256/2) =1.5	atan2(192, 3.4315E-12) =1.5708	180*1.5708/pi =90.0002

（图片横坐标为第几个点，即Hz；纵坐标为模）

5.4 波特图

1、伯德图的定义

伯德图由对数幅频特性和对数相频特性两条曲线构成。 伯德图一般是由二张图组合而成， 伯德图由两张图组成：①G(jω)的幅值(以分贝，dB表示)-频率(以对数标度)对数坐标图，其上画有对数幅频曲线；②G(jω)的相角-频率(以对数标度)对数坐标图，其上画有相频曲线。 

对数幅值的标准表达式为20 lg|G(jω)|，单位是分贝，相角的单位是度。由于增益用对数来表示（log(ab)=log(a)+log(b)），因此一传递函数乘以一常数，在伯德增益图只需将图形的纵向移动即可，二传递函数的相乘，在波德幅频图就变成图形的相加。幅频图纵轴0分贝以下具有正增益裕度、属稳定区，反之属不稳定区。

配合波德相频图可以估算一信号进入系统后，输出信号及原始信号的比例关系及相位。例如一个Asin(ωt) 的信号进入系统后振幅变原来的k倍，相位落后原信号Φ，则其输出信号则为(Ak)sin(ωt−Φ)，其中的k和Φ都是频率的函数。相频图纵轴-180度以上具有正相位裕度、属稳定区，反之属不稳定区。

2、伯德图的画法

画伯德图时，分三个频段进行，先画幅频特性，顺序是中频段、低频段和高频段。将三个频段的频率特性(或称频率响应)合起来就是全频段的幅频特性，然后再根据幅频特性画出相应的相频特性来。 作伯德图时，首先写出频率特性，然后按常数因子K、积分和微分因子(jω)、一阶因子(1+jωT)和二阶因子[1+2ζ(jω/ωn)+(jω)/ω]1这样四种基本因子分别画出伯德图，再总加而成。

下图是常见简易传递函数的伯德图趋势：

五、代码实现（C语言）

FFT.H

FFT.C

六、代码实现（Matlab）

基于MATLAB/Simulink的2ASK数字带通传输系统建模与仿真

相关文章

十分简明易懂的FFT（快速傅里叶变换）

彻底搞懂快速傅里叶变换FFT核心思想--分而治之

数字信号处理公式变程序（一）——DFT、FFT

实验笔记【3】| 利用matlab进行fft分析

信号处理之快速傅里叶变换(FFT)-通俗易懂

深入浅出的理解频谱泄露

快速傅里叶变换（FFT）和小波分析在信号处理上的应用