快速傅里叶变换原理介绍及MATLAB实现

最新推荐文章于 2023-10-16 14:41:41 发布

YHCANDOU

最新推荐文章于 2023-10-16 14:41:41 发布

阅读量2k

点赞数 7

分类专栏：数字信号处理

本文链接：https://blog.csdn.net/yhcwjh/article/details/119173535

版权

数字信号处理专栏收录该内容

2 篇文章 1 订阅

订阅专栏

FFT算法主要有两种，按时间抽选的FFT的算法（DIT-FFT）和按频率抽选的FFT算法（DIF-FFT）。

时间抽取基-2 FFT算法

设 $N={{2}^{M}}$ ，有限长序列 $x\left( n \right)$ 的离散傅里叶变换(DFT)为
$X(k)=\sum\limits_{n=0}^{N-1}{x\left( n \right){{e}^{-jnk\frac{2\pi }{N}}}},0\le k\le N-1$
直接使用DFT，计算一个 $X\left( k \right)$ 值，需计算：N次复数相乘和 $\left( N-1 \right)$ 次复数相加。计算N点 $X\left( k \right)$ 值，需计算： ${{N}^{2}}$ 次复数相乘和 $N\left( N-1 \right)$ 次复数相加。可以看出，在DFT计算中，不论是乘法还是加法，运算量均与 ${{N}^{2}}$ 成正比。这给处理器的运算性能带来了极大的挑战。令
${{W}_{N}}={{e}^{-j\frac{2\pi }{N}}}$
那么
$X(k)=\sum\limits_{n=0}^{N-1}{x\left( n \right)W_{N}^{nk}},0\le k\le N-1$
将序列 $x\left( n \right)$ 按下标偶数项和奇数项分解为两组，偶数项为一组，奇数项为一组，得到两个 $N / 2$ 点的子序列，即
$\left\{ \begin{aligned} & {{x}_{1}}\left( r \right)=x\left( 2r \right) \\ & {{x}_{2}}\left( r \right)=x\left( 2r+1 \right) \\ \end{aligned} \right.r=0,1,\cdots ,N/2-1$
则
$X\left( k \right)=\sum\limits_{r=0}^{N/2-1}{x\left( 2r \right)}W_{N}^{2rk}+\sum\limits_{r=0}^{N/2-1}{x\left( 2r+1 \right)}W_{N}^{(2r+1)k}$
又因为
$W_{N}^{2rk}={{e}^{-j\frac{2\pi }{N}2rk}}={{e}^{-j\frac{2\pi }{N/2}rk}}=W_{N/2}^{rk}$
所以
$X\left( k \right)=\sum\limits_{r=0}^{N/2-1}{x\left( 2r \right)}W_{N/2}^{rk}+\sum\limits_{r=0}^{N/2-1}{x(2r+1)W_{N}^{k}W_{N/2}^{rk}}$
上式又可以写作
$\begin{aligned} & X\left( k \right)=\sum\limits_{r=0}^{N/2-1}{{{x}_{1}}}\left( r \right)W_{N/2}^{rk}+W_{N}^{k}\sum\limits_{r=0}^{N/2-1}{{{x}_{2}}}\left( r \right)W_{N/2}^{rk} \\ & \begin{matrix} \begin{matrix} {} & {} \\ \end{matrix} & ={{X}_{1}} \\ \end{matrix}\left( k \right)+W_{N}^{k}{{X}_{2}}\left( k \right) \\ \end{aligned}$
同理可得
$X\left( k+\frac{N}{2} \right)={{X}_{1}}\left( k \right)-W_{N}^{k}{{X}_{2}}\left( k \right)\begin{matrix} {} & k=0,1,...,N/2-1 \\ \end{matrix}$
那么，傅里叶变换可以写作
$\left\{ \begin{aligned} & X\left( k \right)={{X}_{1}}\left( k \right)+W_{N}^{k}{{X}_{2}}\left( k \right) \\ & X\left( k+\frac{N}{2} \right)={{X}_{1}}\left( k \right)-W_{N}^{k}{{X}_{2}}\left( k \right) \\ \end{aligned} \right.\begin{matrix} {} & k=0,1,...,N/2-1 \\ \end{matrix}$
通过M次分解最后分解成2点的DFT运算。则，时间抽取基-2 FFT总共需要的运算量为：复乘： $N/2\cdot M=N/2\cdot {{\log }_{2}}N$ ,复加： $N\cdot M=N{{\log }_{2}}N$ 。显然，要比直接DFT运算快得多。

频率抽取基-2 FFT算法

DIF-FFT算法是将输入序列 $X\left( k \right)$ 分成前后两个部分
$\begin{aligned} & X\left( k \right)=\sum\limits_{n==0}^{N-1}{x\left( n \right)W_{N}^{nk}}=\sum\limits_{n=0}^{\frac{N}{2}-1}{x\left( n \right)W_{N}^{nk}}+\sum\limits_{n=\frac{N}{2}}^{N-1}{x\left( n \right)W_{N}^{nk}} \\ & =\sum\limits_{n=0}^{N-1}{x\left( n \right)W_{N}^{nk}}+\sum\limits_{n=0}^{N-1}{x\left( n+\frac{N}{2} \right)W_{N}^{k\left( n+\frac{N}{2} \right)}} \\ & =\sum\limits_{n=0}^{\frac{N}{2}-1}{\left[ x\left( n \right)+x\left[ n+\frac{N}{2} \right]W_{N}^{Nk/2} \right]W_{N}^{nk}} \end{aligned}$

因为 $W_{N}^{N/2}=-1$ ，则 $W_{N}^{Nk/2}={{\left( -1 \right)}^{k}}=\left\{ \begin{matrix} 1 & k为偶数 \\ -1 & k为奇数 \\ \end{matrix} \right.$ ，所以可以得到
$X(k)=\sum\limits_{n=0}^{\frac{N}{2}-1}{\left[ x\left( n \right)+{{\left( -1 \right)}^{k}}x\left( n+\frac{N}{2} \right) \right]W_{N}^{nk}}$

将k按照奇数和偶数来分，得到
$\left\{ \begin{matrix} k=2r & k为偶数 \\ k=2r+1 & k为奇数 \\ \end{matrix} \right.r\text{=0,1,}\cdots \text{N/2-1}$
将X(k)分为两部分
$\left\{ \begin{aligned} & X\left( 2r+1 \right)=\sum\limits_{n=0}^{\frac{N}{2}-1}{\left[ x\left( n \right)-x\left( n+\frac{N}{2} \right) \right]W_{N}^{n}W_{N/2}^{rn}} \\ & X\left( 2r \right)=\sum\limits_{n=0}^{\frac{N}{2}-1}{\left[ x\left( n \right)+x\left( n+\frac{N}{2} \right) \right]W_{N/2}^{rn}} \\ \end{aligned} \right.$
令 ${{x}_{1}}=x\left( n \right)+x\left( n+\frac{N}{2} \right)$ ， ${{x}_{2}}\left( n \right)=\left[ x\left( n \right)-x\left( n+\frac{N}{2} \right) \right]W_{N}^{n}$ ，可以得到
$\left\{ \begin{aligned} & X\left( 2r+1 \right)=\sum\limits_{n=0}^{\frac{N}{2}-1}{{{x}_{2}}\left( n \right)W_{N/2}^{rn}} \\ & X\left( 2r \right)=\sum\limits_{n=0}^{\frac{N}{2}-1}{{{x}_{1}}\left( n \right)W_{N/2}^{rn}} \\ \end{aligned} \right.$
则经过M次分解，同样可以将DFT分解成为N/2个两点DFT。以8点DIF-FFT运算为例，DIF-FFT运算流图如下图所示
在这里插入图片描述
从上图可以看出，输出序列的排列规律不是从小到大按顺序的，而是按照倒叙规则排序的，即先将0-7转换为二进制数，然后将二进制数左右倒序，再转为十进制就可以得到新的数列，即：0，4，2，6，1，5，3，7。

利用MATLAB实现DFT

DFT的矩阵表示形式为
$\left[ \begin{matrix} X\left( 0 \right) \\ X\left( 1 \right) \\ X\left( 2 \right) \\ \begin{matrix} \cdots \\ X\left( N-1 \right) \\ \end{matrix} \\ \end{matrix} \right]=\left[ \begin{matrix} W_{N}^{0} & W_{N}^{0} & W_{N}^{0} & \begin{matrix} \cdots & W_{N}^{0} \\ \end{matrix} \\ W_{N}^{0} & W_{N}^{1} & W_{N}^{2} & \begin{matrix} \cdots & W_{N}^{N-1} \\ \end{matrix} \\ W_{N}^{0} & W_{N}^{2} & W_{N}^{4} & \begin{matrix} \cdots & W_{N}^{2\left( N-1 \right)} \\ \end{matrix} \\ \begin{matrix} \cdots \\ W_{N}^{0} \\ \end{matrix} & \begin{matrix} \cdots \\ W_{N}^{N-1} \\ \end{matrix} & \begin{matrix} \cdots \\ W_{N}^{2\left( N-1 \right)} \\ \end{matrix} & \begin{matrix} \begin{matrix} \cdots & \cdots \\ \end{matrix} \\ \begin{matrix} \cdots & W_{N}^{\left( N-1 \right)\left( N-1 \right)} \\ \end{matrix} \\ \end{matrix} \\ \end{matrix} \right]\left[ \begin{matrix} x\left( 0 \right) \\ x\left( 1 \right) \\ x\left( 2 \right) \\ \begin{matrix} \cdots \\ x\left( N-1 \right) \\ \end{matrix} \\ \end{matrix} \right]$
根据DFT公式和矩阵展开形式，利用MATLAB可以得到如下结果
在这里插入图片描述

利用MATLAB实现FFT（8点DIF-FFT为例）

由DIF-FFT的流程图可以得到，FFT运算的基本单元是蝶形运算单元。蝶形运算单元可以用简单的语句实现，然后可以用循环语句来进行各个蝶形运算单元的运算。
8点FFT的蝶形运算有3级，第1级有1组，每组4个蝶形运算单元，旋转因子是 $W_{8}^{0}$ 、 $W_{8}^{1}$ 、 $W_{8}^{2}$ 、 $W_{8}^{3}$ ；第2级有2组，每组2个蝶形运算单元，循环因子是 $W_{4}^{0}$ ， $W_{4}^{1}$ ；第3级有4组，每组1个蝶形运算单元，旋转因子是 $W_{2}^{0}$ 。总结运算规律，来设定循环语句。
第一层循环在1到3级间循环，循环变量 $m m = 1, 2, 3$ 。旋转因子下标 $N_{m}=2^{4-mm}=8,4,2$ 。
第二层循环在该级的各组间循环，每级有 $2^{mm-1}$ 组，每组第一行对应的 $x$ 值为：第1级是 $x (0)$ ，第2级是 $x (0)$ ， $x (4)$ ，第3级是 $x (0)$ ， $x (2)$ ， $x (4)$ ， $x (6)$ 。设第二层循环变量为 $p$ ，则在MATLAB中， $p=0:N_{m}:7$ 。
第三层循环在该组的各个蝶形运算单元间循环，每组有 $\frac{{{N}_{m}}}{2}$ 个蝶形运算单元，则循环变量 $k$ 从1到 $\frac{{{N}_{m}}}{2}$ ，旋转因子是 ${{e}^{\frac{-2\pi (k-1)j}{Nm}}}$ ，每次蝶形运算跨越的行数是 $\frac{{{N}_{m}}}{2}$ ，则参加蝶形运算的 $x$ 值为 $x (k + p)$ 和 $x(k+p+\frac{{{N}_{m}}}{2})$ 。
循环完成后则FFT运算完成，再将x序列按倒叙规律重新排列就可以得到 $X (k)$ 序列。

涉及到的MATLAB函数

exp(x)计算e的x次方,abs(x)计算x绝对值。

倒序运算主要用到的函数有，dec2bin(x,m)，是把十进制序列x转换为m位二进制数，bin2dec(x,m)也是类似功能。

fliplr函数是将一个矩阵左右颠倒，则程序中求倒序的语句就是：d=bin2dec(fliplr(dec2bin([0:N-1],m)))+1;y=x(d); 其中，x是N点序列，执行完语句后，y序列就是x序列的倒序排列。

根据上述流程，可以得到结果如下
在这里插入图片描述

调用MATLAB的FFT函数

MATLAB自带的FFT函数语法格式为：A=fft(X,N)。A返回信号X的FFT变换矩阵，输入信号X和输出信号A大小相同。N为FFT的点数，当X长度不足N，自动在数据末尾补零，当X长度超过N，则进行自动截取。结果如下在这里插入图片描述
通过对比可以发现，和上述实现的FFT结果是一致的。

实验一代码如下：

clear;
close all;
clc;
N=256;                      %信号点数
n=0:N-1;                    %时间采样
xn=cos(pi*n/6);           %信号
k=0:N-1;                    %频域采样
WN=exp(-1j*2*pi/N);
nk=n'*k;
WNnk=WN.^nk;
Xk=xn*WNnk;            %计算DFT
figure(1);
subplot(2,1,1);stem(xn);title('变换序列');xlim([0 50])
subplot(2,1,2);stem(abs(Xk));title('序列的DFT');

实验二代码如下：

clear;
close all;
clc;
N=8;                                                           %FFT点数
n=0:N-1;
x=[0 1 2 2 1 0 1 2];                                      %变换序列
figure(1);
subplot(2,1,1);stem(x);title('变换序列');
x1=x;
m=log2(N);                                                 %蝶形运算的级数
for mm=1:m                                               %循环蝶形运算
    Nm=2^(m-mm+1);                                 %旋转因子下标Nm=8， 4， 2
    for p=0:Nm:N-1                                       %循环该级1~2mm-1组蝶形运算
        for k=1:Nm/2
            kp=k+Nm/2+p;                               %蝶形运算的下标
            a=x(kp);
            x(kp)=(x(k+p)-a)*exp(-1j*2*pi*(k-1)/Nm);
            x(k+p)=x(k+p)+a;
        end
    end
end
d=bin2dec(fliplr(dec2bin([0:N-1],m)))+1;               %倒序排列
y=x(d);
mag=abs(y);
subplot(2,1,2);stem(mag);title('FFT结果');
x=[0 1 2 2 1 0 1 2];
yy=fft(x,N);
figure(2);
subplot(2,1,1);stem(x);title('变换序列');
subplot(2,1,2);stem(abs(yy));title('FFT结果');