基4FFT 原理及MATLAB 算法实现

图 5 运算时间测量
使用 tic，toc 命令对自定义基 4FFT 和系统 FFT 进行测量发现， 系统 FFT 远远快于自
6
定义基 4FFT。如图 5 所示，每次运行程序，第一行为系统 FFT 运算时间，第二行为自定义 基 4FFT 运算时间。 分析认为，这是由于两方面原因造成的。一是 matlab 本身是解释性的高级语言，与 vb 属于同一种类型，运行时，要先对语句命令进行解释，翻译成计算机可识别的机器语言，然 后再执行。这就大大增加了自定义基 4FFT 运算时间。而系统自带的 FFT 显然是预先编译好 的内部代码，因此执行效率非常高。二是程序本身未最大限度的优化。例如本程序中，将递 推公式中重复计算的部分用变量代换： 令 X1=x(k1);X2=Wn^m*x(k2);X3=Wn^(2*m)*x(k3);X4=Wn^(3*m)*x(k4); 递推公式化简为： temp(k1)=X1+X2+X3+X4; temp(k2)=X1-1j*X2-X3+1j*X4; temp(k3)=X1-X2+X3-X4; temp(k4)=X1+1j*X2-X3-1j*X4; 从而避免了 X2，X3，X4 表达式中的乘法在 4 个递推表达式中各被重复计算 4 次。如 未作这样的处理，显然效率很低。由于个人能力原因，本程序仍有很大优化空间。 3.进一步提高运算速度的方法 1）采用 c 语言编程，编译成机器代码，大大提高执行效率。 2）事先建立旋转因子表，以空间换取时间的方法，提高运算速度。计算过程中，通过 查表的方式，获取旋转因子值，避免实时乘法运算增加的运算时间。 四、总结 本文在分析基 2FFT，参考借鉴多篇文章的基础上，对基 4FFT 的原理进行了分析，推 导出 4 点迭代运算公式，运用 matlab 实现该算法。所编程序可以对 4L 点采样数据进行 4L 点基 4FFT （L=1~8） 。在此基础上，对比了系统自带 FFT 与自定义基 4FFT 的运算效率。 提出了进一步提高运算效率的方法。
N
4
L l
，第 g 组对应的序号为 k k0 g
N
令 k1 k， k2 k1 如下信号流图
N
4
Ll
， k3 k2
N
4
Ll
， k4 k3
N
4Ll
，则可将公式（5）进一步简化，并得到
Xl‐1（
）
Xl（
）
Xl‐1（
）
Xl（
）
Xl‐1（
）
Xl（
）
k1=k;k2=k1+group_interval_1;k3=k2+group_interval_1;k4=k3+group_interval_1; X1=x(k1);X2=Wn^m*x(k2);X3=Wn^(2*m)*x(k3);X4=Wn^(3*m)*x(k4); 式中重复计算部分用变量代换，减少运算次数 %根据递推公式计算，结果存储到临时数组 temp temp(k1)=X1+X2+X3+X4; temp(k2)=X1-1j*X2-X3+1j*X4; temp(k3)=X1-X2+X3-X4; temp(k4)=X1+1j*X2-X3-1j*X4; end end x=temp; end subplot(2,2,4),plot(abs(x)); axis([0 4096 0 2500]),title('自定义基 4FFT 计算出的频谱'); %将 temp 中临时存储的第 l 级结果赋值给 x，作为次级运算的输入
2N 4
L l 1
（5.2）
X l(k
2N 4
L l 1
) X l 1(k ) W Nm X l 1(k
N
4
L l 1
) W N2m X l 1(k
) W N3m X l 1(k
3N 4L l 1 3N 4
L l 1
)
（5.3）
X l(k
0
0 0 0
3
K 1
(4n n 0 )
，N
K 4
（3）
X(M )
n0 0
W Kmn )] W 4m 0n0 ， M m Km 0 [W Nmn0( x[4n n0 ]
n 0
3
K 1
（4）
由上式可得如下矩阵变换过程：
0 1 2 3
, 乘以 , , ,
4 5 6 7

WN x[n ] n
0
N
mn
（1）
当 N 等于 4 时的四点 DFT 运算为： 1 1 = 1 1 1 1 1 1 1 1 1 1 ∙ （2）
∙
可以看到，类似 2 点 DFT，4 点 DFT 运算也无需乘法，可以简少运算量。 对(1)式进行分解：
X (M )
x[4n n ] W NM n n
∙ ∙ ∙ ∙
… … … … … … … …
, , , ,
4 3 2 1 ∙ ∙ ∙ ∙
∙ ∙ ∙ ∙
, 行 , , ,
, , , ,
… … … …
, , , ,
∙ ∙ ∙ ∙
∙ ∙ ∙ ∙
, , , ,
∙ ∙ ∙ ∙
0
列
1 2 3 1 1 1
2 3
… … … …
2 3
1 1 1 1
可以看到，第一步先对最开始的采样点矩阵每一行进行 K 点 DFT，然后第二步每项对 应乘以旋转因子 ，n0 为行（0 到 3 行） ，m 为列（0 到 K-1 列）,最后按列做 4 点 DFT, 得到一个按行顺序排列的最终结果。 显然，第一步中的行 DFT 可以进一步分解成两步，即第一步进行 K/4 点 DFT，第二步 乘以对应旋转因子 ，即 ∙ ，n0 为行（0 到 3 行） ，m 为列（0 到 K/4-1 列）, 按列
3N 4
L l 1
) X l 1(k ) jWNm X l 1(k
N
4
L l 1
) W N2m X l 1(k
2N 4
L l 1
) jWN3m X l 1(k
)
（5.4）
k k0 g
N
4
L l
， （ k0 0, 1, 2, …,
N
4
Ll1
%分组上限 %组内数据上限
%每一级中包含的组循环，遍历每一组，计算各组中的数据
for k0=0:K %遍历每一组中的所有数据，计算次级数据 4
k=k0+g*group_interval_2+1; m=group_cont_2*k0; 数据中每一级所乘旋转因子的指数因子
基 4FFT 原理及 MATLAB 算法实现
FFT( 快速傅里叶变换 ) 算法与 DFT( 离散傅里叶变换 ) 算法比较 , 其运算量显著减少 , 用计算机实现时速度大为提高。其思想是利用 2 点 DFT 运算无需乘法的特点，以减少过程 在乘法运算上的时间开销。考虑到 4 点 DFT 也无需乘法，可以减少过程中的乘法次数，且 运算级数减少，无疑能加快 DFT 的实现。 通常的FFT幂法都是“基2分解法”，即长度为N的D FT序列由两个长度为N / 2的DFT 列的组合表示；而这两个长度为N / 2的DFT序列各自又分别由两个长度为N / 4的D F T序列 的组合表示；⋯⋯，按照这一做法对序列进行反复分解，直到每个序列的长度等于2为止。 这个分解、 组合过程如同一棵标准二叉树。 按分解的逆过程进行组合运算便得到所要求的频 谱序列X( n)，(n =0，l，⋯，N一1)整个变换过程共需要进行N/2log2N次复数乘法运算。 参考上述的分解、组合方法，对序列进行“基4分解” ，即长度为N的DFT序列由四个长 度为N/4的DFT序列组合表示。 一、基本原理 N 点 DFT 公式为： X ( m )
1）,（g=0, 1, 2, …, 4Ll 1）,（ l =1, 2, 3, ..., L ）
m 4Llk0
。 4L l 可以验证，对于第一级 4 点行 DFT，递推公式仍适用。因此可以通过递推公式（5）反 复迭代 L 次得到最终 N 点 DFT 结果。 每一级中分为 g= 4L l 组，组与组间距为
数据计算循环：for
0：
；
利用当前级数据 X， 利用递推公式计算出次级数据并存入临时 数组 temp，最后用临时数组中的次级数据覆盖 X
最终得到 N 点 DFT 结果 X
图 2 程序流程图 3
2.matlab 程序
clear; a=0:4095;x=sin(2*pi/3*a)+sin(2*pi/4*a)+sin(2*pi/5*a); subplot(2,1,1),plot(x); axis([0 4096 -3 3]),title('时域信号波形'); subplot(2,2,3),plot(abs(fft(x))); axis([0 4096 0 2500]),title('系统 FFT 计算出的频谱'); N=4096; %N 点 DFT，N 为 4 的整数次幂 L=log(N)/log(4); Wn=exp(-2j*pi/N); temp=zeros(1,N); %4 点 DFT 分解级数 %旋转因子 %定义中间临时数组 %输入三频率信号
N
4
L l 1
) W N2m X l 1(k
L l 1
) W N3m X l 1(k
3N 4
L l 1
)
（5.1）
X l(k
N
4
L l 1
) X l 1(k ) jWNm X l 1(k
N
4
L l 1
) WN2m X l 1(k