【数字信号处理】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,?,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& ?& 1\\ 1& W& W^2& W^3& ?& W^{N?1}\\ 1& W^2& W^4& W^6& ?& W^{2(N?1)}\\ 1& W^3& W^6& W^9& ?& W^{3(N?1)}\\ ?& & & & & \\ 1& W^k& W^{2k}& W^{3k}& ?& W^{k(N?1)}\\ ?& & & & & \\ 1& W^{N?1}& W^{2(N?1)}& W^{3(N?1)}& ?& 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^{?},& a^4=?1,\\ a^5=?a,& a^6=j,& a^7=ja=a^{?},& 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?4=20\\ F[1]=& 0+W_4^1?4=?4j\\ F[2]=& 16+W_4^2?4=12\\ F[3]=& 0+W_4^3?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?g)[k]=\sum _{n=0}^{k}f[n]g[k?n],k=0,1,?,2N?2$$
卷积的结果也是一个序列。
现在设有多项式：
$$F(x)=f[0]+f[1]x+f[2]x^2+?+f[N?1]x^{N?1}$$
$$G(x)=g[0]+g[1]x+g[2]x^2+?+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+?\\ & \\ =& \sum _{n=0}^{k}f[n]g[k?n]x^k,k=0,1,?,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?g=[7,2,7,2,0]\end{array}$$

1.2.2 多项式与DFT

对于多项式的系数序列表示：
$$F(x)=f[0]+f[1]x+f[2]x^2+?+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,?,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?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)

---

下链

---