【数字信号处理】FFT

FFT

2023年11月18日
#elecEngeneer

【数字信号处理】DFT

---

---

1. 快速傅里叶变换-FFT

由于乘法是最慢的操作，衡量DFT在计算机上的标准就是乘法是数量。在前面的分析中DFT可以看成是矩阵乘以向量，所以乘法数量是 $N^2$ ， $N$ 是采样点的个数，也是变换的长度，对于大多数稳态， $N$ 被选为最少是 $256$ ，来获得对信号有良好近似的频谱，计算速度也因此变得非常重要。
上世纪60年代中期开始，高效率执行DFT的算法逐渐出现。这些算法就是FFT算法，这些算法基于一个事实，即DFT存在大量冗余的计算。
将DFT写成：
$$F[k]\stackrel{\text{ DFT }}{=}\sum _{n=0}^{N-1}f[n]W_N^{kn},k=0,1,\cdots ,N-1$$
$$W_N=e^{-j\frac{2\pi }{N}}$$

[!quote]- DFT矩阵
$$\left[\begin{array}{c}F[0]\\ F[1]\\ F[2]\\ ⋮\\ F[k]\\ ⋮\\ F[N-1]\end{array}\right]=\left[\begin{array}{cccccc}1& 1& 1& 1& \cdots & 1\\ 1& W& W^2& W^3& \cdots & W^{N-1}\\ 1& W^2& W^4& W^6& \cdots & W^{2(N-1)}\\ 1& W^3& W^6& W^9& \cdots & W^{3(N-1)}\\ ⋮& & & & & \\ 1& W^k& W^{2k}& W^{3k}& \cdots & W^{k(N-1)}\\ ⋮& & & & & \\ 1& W^{N-1}& W^{2(N-1)}& W^{3(N-1)}& \cdots & W^{(N-1)(N-1)}\end{array}\right]\left[\begin{array}{c}f[0]\\ f[1]\\ f[2]\\ ⋮\\ f[n]\\ ⋮\\ f[N-1]\end{array}\right]$$

从DFT矩阵可以看出，相同的 $W_N^{kn}$ 会被计算很多次，比如 $kn=16$ ，可能有
$$\begin{array}{rl}n=2,& k=8\\ n=8,& k=2\\ n=1,& k=16\\ n=16,& k=1\\ n=4,& k=4\end{array}$$
这几种情况，而且 $W_N$ 的整数幂是一个周期为 $N$ 的函数，只有 $N$ 个确定的值。比如 $N=8$ （FFT最简单的时候 $N$ 是 $2$ 的整数幂）：
$$W_8^1=e^{-j2\pi /8}=e^{-j45°}=\frac{1-j}{\sqrt{2}}=a$$
$$\begin{array}{cccc}{a}^1=a,& a^2=-j,& a^3=-ja=-a^{\ast },& a^4=-1,\\ a^5=-a,& a^6=j,& a^7=ja=a^{\ast },& a^8=1,\\ a^9=a^{8}a=a& & & \end{array}$$
所以：
$$\begin{array}{rl}{W}_8^4=& -W_8^0\\ W_8^5=& -W_8^1\\ W_8^6=& -W_8^2\\ W_8^7=& -W_8^3\end{array}$$
实际上 $W_N$ 就是个幅值为 $1$ 的复数。 $N$ 相当于把复平面上的单位圆分成 $N$ 等分，如下图：

注意，从定义可知 $W_N^x$ 是逆时针排列的！！！
对于复数，负号相当于旋转 $180°$ ，当 $N$ 为 $2$ 的整数幂时， $W_N^x$ 存在落在坐标轴上的点。

1.1 时间抽取FFT（Decimation-in-time algorithm）

也叫“库利-图基算法”。将DFT采样序列中 $n$ 为偶数和 $n$ 为奇数的项分开，设 $N$ 为 $2$ 的整数次幂，则分开后的偶数项和奇数项各有 $N/2$ 个：
$$\begin{array}{rl}F[k]\stackrel{\text{ DFT }}{=}& \sum _{n=0}^{N-1}f[n]W_N^{kn}\\ & \\ =& \sum _{m=0}^{N/2-1}f[2m]W_N^{k2m}+\sum _{m=0}^{N/2-1}f[2m+1]W_N^{k(2m+1)}\end{array}$$
而
$$W_N^{k2m}=e^{-j\frac{2\pi }{N}2mk}=e^{-j\frac{2\pi }{N/2}mk}=W_{N/2}^{km}$$
$$F[k]=\sum _{m=0}^{N/2-1}f[2m]W_{N/2}^{km}+W_N^k\sum _{m=0}^{N/2-1}f[2m+1]W_{N/2}^{km}$$
$$F[k]=G[k]+W_N^{k}H[k]$$
因此 $N$ 个采样点的DFT可以由两个 $N/2$ 个采样点的DFT计算得到，其中一个使用的是偶数的采样点，另一个使用的是奇数的采样点。计算量显然减少了。对于 $N=8$ ：
$$\{\begin{array}{l}偶数的采样点:f[0],f[2],f[4],f[6]\\ \\ 奇数的采样点:f[1],f[3],f[5],f[7]\end{array}$$
$$\begin{array}{rl}F[0]=& G[0]+W_8^{0}H[0]\\ F[1]=& G[1]+W_8^{1}H[1]\\ F[2]=& G[2]+W_8^{2}H[2]\\ F[3]=& G[3]+W_8^{3}H[3]\\ F[4]=& G[0]+W_8^{4}H[0]=G[0]-W_8^{0}H[0]\\ F[5]=& G[1]+W_8^{5}H[1]=G[1]-W_8^{1}H[1]\\ F[6]=& G[2]+W_8^{6}H[2]=G[2]-W_8^{2}H[2]\\ F[7]=& G[3]+W_8^{7}H[3]=G[3]-W_8^{3}H[3]\end{array}$$
FFT的数据流图如下：

由于 $N$ 是 $2$ 的整数幂， $N/2$ 仍然是偶数，因此可以继续将 $G[k]$ 分成两个变换组成， $H[k]$ 也分成两个变换组成，每个变换含有 $N/4$ 个采样点。
可见，时间抽取FFT可以将DFT递归地分解成许多小的变换。

$A$ 和 $B$ 是复数，因此一个蝴蝶式的计算需要一个复数乘法（ $W_N^{p}B$ ）和两个复数加法。
一个重要的发现，或者技巧，就是分解到最后，采样序列的索引是 位反转（bit-reversed） 的，如下表

对于 $N$ 个采样点，DFT需要 $N^2$ 个复数乘法，FFT需要大约 $\frac{N}{2}{\log }_2(N)$ 个。

[!example]-
$$f:[8,4,8,0]$$

$$\begin{array}{rl}F[0]=& 16+W_4^0\cdot 4=20\\ F[1]=& 0+W_4^1\cdot 4=-4j\\ F[2]=& 16+W_4^2\cdot 4=12\\ F[3]=& 0+W_4^3\cdot 4=4j\end{array}$$

1.2 FFT做多项式乘法（卷积）

1.2.1 多项式乘法与卷积

序列的卷积 两个序列 $f[n],n=0:(M-1)$ 和 $g[n],n=0:(N-1)$ 的卷积，表示 $g$ 滑过 $f$ 时依据这些点确定的重叠部分的面积。
如果 $M<N$ ，则设 $f[n]=0,n=M:(N-1)$ ，即补零法，依此得到两个长度为 $N$ 的等长序列。卷积的公式如下：
$$(f\ast g)[k]=\sum _{n=0}^{k}f[n]g[k-n],k=0,1,\cdots ,2N-2$$
卷积的结果也是一个序列。
现在设有多项式：
$$F(x)=f[0]+f[1]x+f[2]x^2+\cdots +f[N-1]x^{N-1}$$
$$G(x)=g[0]+g[1]x+g[2]x^2+\cdots +g[N-1]x^{N-1}$$
可以发现两个多项式相乘：
$$\begin{array}{rl}F(x)\times G(x)=& f[0]g[0]+(f[0]g[1]+f[1]g[0])x+(f[0]g[2]+f[1]g[1]+f[2]g[0])x^2+\cdots \\ & \\ =& \sum _{n=0}^{k}f[n]g[k-n]x^k,k=0,1,\cdots ,2N-2\end{array}$$
即多项式相乘后的系数序列是原来两个多项式系数序列的卷积。

[!example]-
$$\begin{array}{rl}F(x)=& 1+0x+x^2,f=[1,0,1]\\ & \\ G(x)=& 7+2x+0x^2,g=[7,2,0]\end{array}$$
$$\begin{array}{r}(F\times G)(x)=7+2x+7x^2+2x^3+0x^4,f\ast g=[7,2,7,2,0]\end{array}$$

1.2.2 多项式与DFT

对于多项式的系数序列表示：
$$F(x)=f[0]+f[1]x+f[2]x^2+\cdots +f[N-1]x^{N-1}$$
如果令
$$x=e^{-j\omega T}$$
则
$$\begin{array}{rl}F(x)=& \sum _{n=0}^{N-1}f[n]e^{-j\omega nT}\stackrel{\text{ DTFT }}{=}\bar{F}(j\omega )\end{array}$$
所以使用系数序列表示的多项式，相当于系数序列的离散时间傅里叶变换（DTFT）。
而根据多项式插值定理， $N-1$ 阶多项式可以被 $N$ 个点唯一确定，所以 $N-1$ 阶多项式除了通过系数序列表示，还可以通过多项式上的 $N$ 个点来表示。即多项式可以表示成
$$\begin{array}{r}F(x)=\sum _{n=0}^{N-1}f[n]e^{-j\frac{2\pi }{N}kn}\stackrel{\text{ DFT }}{=}\bar{F}[k]\end{array}$$
$$x=e^{-j\omega T}=e^{-j(\frac{2\pi }{NT})kT}=e^{-j\frac{2\pi }{N}k},k=0,1,\cdots ,N-1$$
所以使用 $N$ 个点来表示的多项式，相当于系数序列的DFT。结合多项式乘法来看，系数序列的卷积为系数序列傅里叶变换后的乘积，符合傅里叶变换的卷积性质。

1.2.3 多项式乘法与FFT

对使用 $N$ 个点表示的多项式，多项式的乘积会变得相对简单，因为我们可以快速确定新的多项式上的 $N$ 个点。
$$(x_i,F(x_i)\times G(x_i))或(k,F[k]\times G[k])$$
得到的是多项式乘积的点值表示，对其进行反FFT，即可得到多项式乘积的系数表示。
由于 $N$ 个点只能确定 $N-1$ 次多项式，在计算多项式乘法时，我们需要对多项式进行补 $0$ ，如果要用FFT，需要补成 $2$ 的整数次幂长度。

[!example]-
$$\begin{array}{rl}F(x)=& 1+2x+3x^2+4x^3\\ G(x)=& 2+3x+4x^2+5x^3\end{array}$$
$$\begin{array}{r}f:[1,2,3,4,0,0,0,0]\\ g:[2,3,4,5,0,0,0,0]\end{array}$$
$$F[k]\stackrel{\text{ FFT }}{=}[10,-0.4142-7.2426j,-2+2j,2.4142-1.2426j,-2,2.4142+1.2426j,-2-2j,-0.4142+7.2426j]$$
$$G[k]\stackrel{\text{ FFT }}{=}[14,0.5858-9.6569j,-2+2j,3.4142-1.6569j,-2,3.4142+1.6569j,-2-2j,0.5858+9.6569j]$$
$$f\ast g=\text{IFFT}(F\times G)=[2,7,16,30,34,31,20,0]$$
$$\therefore F(x)\times G(x)=2+7x+16x^2+30x^3+34x^4+31x^5+20x^6$$

clc;clear
f = [1 2 3 4 0 0 0 0];
F = fft(f)
g = [2 3 4 5 0 0 0 0];
G = fft(g)
T = F.*G
ifft(T)
syms x
F1 = 1+2*x+3*x^2+4*x^3;
G1 = 2+3*x+4*x^2+5*x^3;
T1 = F1*G1;
collect(T1,x)

---

下链

---