diff --git a/_config.yml b/_config.yml index ebf4b17..22dd983 100644 --- a/_config.yml +++ b/_config.yml @@ -200,7 +200,7 @@ neat_js: related_posts: enabled: true filter_threshold: 0.13 - related_count: 5 + related_count: 3 weight: title: 0.1 description: 0 diff --git a/source/_posts/2023-04-14-josh-dsp-part-1.md b/source/_posts/2023-04-14-josh-dsp-part-1.md index 681baf1..b866c1a 100644 --- a/source/_posts/2023-04-14-josh-dsp-part-1.md +++ b/source/_posts/2023-04-14-josh-dsp-part-1.md @@ -274,7 +274,11 @@ X_e\left(e^{j\omega}\right)=\frac{1}{2}\left[X\left(e^{j\omega}\right)+X^\ast\le # 9. 傅里叶变换、拉氏变换和 $z$ 变换的关系 -![傅里叶变换、拉氏变换和 $z$ 变换的关系](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-04-14-josh-dsp-part-1/2023-04-14-josh-dsp-part-1-010-RelationshipBetweenTransform.png?imageMogr2/thumbnail/!100p|watermark/2/text/QEpvc2ggR2Fv/font/YWhyb25iZC50dGY=/fontsize/14/dissolve/20/gravity/southeast/dx/5/dy/5){width=600px} +![傅里叶变换、拉氏变换和 $z$ 变换的关系](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-04-14-josh-dsp-part-1/2023-04-14-josh-dsp-part-1-010-RelationshipBetweenTransform.png?imageMogr2/thumbnail/!100p|watermark/2/text/QEpvc2ggR2Fv/font/YWhyb25iZC50dGY=/fontsize/80/dissolve/20/gravity/southeast/dx/5/dy/5){width=600px} + +# 参考文献 + +王世一. *数字信号处理, 修订版.* 北京理工大学出版社, 1997. [^1]: $f\left( x \right)$ 的以 $2l$ 为周期的傅氏级数的指数形式为 $$ diff --git a/source/_posts/2023-04-14-josh-dsp-part-2.md b/source/_posts/2023-04-14-josh-dsp-part-2.md index fc15f26..934bb93 100644 --- a/source/_posts/2023-04-14-josh-dsp-part-2.md +++ b/source/_posts/2023-04-14-josh-dsp-part-2.md @@ -245,3 +245,7 @@ $$\begin{equation} \end{equation}$$ > 注：此处用 $\odot$ 表示圆周卷积，不同教材中圆周卷积的符号表述可能不同。 + +# 参考文献 + +王世一. *数字信号处理, 修订版.* 北京理工大学出版社, 1997. diff --git a/source/_posts/2023-04-14-josh-dsp-part-3.md b/source/_posts/2023-04-14-josh-dsp-part-3.md new file mode 100644 index 0000000..adfae0a --- /dev/null +++ b/source/_posts/2023-04-14-josh-dsp-part-3.md @@ -0,0 +1,564 @@ +--- +title: Josh 的学习笔记之数字信号处理（Part 3——离散傅里叶变换 DFT） +mathjax: true +comments: true +copyright: true +toc: + enable: true + number: false + wrap: true + expand_all: true + max_depth: 3 +abbrlink: fd9f9a29 +date: 2023-04-14 22:17:40 +tags: + - 数字信号和处理 + - DSP + - 离散傅里叶变换 + - DFT + - 延长序列 + - 圆周卷积 + - 重叠相加 + - 重叠保留 + - 内插公式 + - 混叠 + - 频谱泄露 + - 栅栏效应 +categories: + - - Josh 的学习笔记 + - 数字信号处理 +--- + +# 1. 离散傅里叶变换（Discrete Fourier Transform, DFT）的定义 + +  取离散傅里叶级数的主值，得离散傅里叶变换 + +$$\begin{equation} +X\left(k\right)=\widetilde{X}\left(k\right)R_N\left(k\right)=\sum_{n=0}^{N-1}{x\left(n\right)W_N^{kn}},\ \ k=0,1,\cdots,N-1 +\end{equation}$$ + +$$\begin{equation} +x\left(n\right)=\widetilde{x}\left(n\right)R_N\left(n\right)=\frac{1}{N}\sum_{k=0}^{N-1}{X\left(k\right)W_N^{-kn}},\ \ n=0,1,\cdots,N-1 +\end{equation}$$ + +且 $\widetilde{X}\left(k\right)=X\left(\left(k\right)\right)_N,\widetilde{x}\left(n\right)=x\left(\left(n\right)\right)_N$。其中 $R_N\left(n\right)$ 表示矩形序列，表示取主值；$\left(\left(n\right)\right)_N$ 表示 $n$ 对 $N$ 取余。 + +> 注：不同教材中取主值、取余的符号表述可能不同。 + + + +# 2. DFT 的矩阵算法 + +  若令 $X=\left[\begin{matrix}X\left(0\right)&X\left(1\right)&\cdots&X\left(N-1\right)\\\end{matrix}\right]^\mathrm{T}，x=\left[\begin{matrix}x\left(0\right)&x\left(1\right)&\cdots&x\left(N-1\right)\\\end{matrix}\right]^\mathrm{T}$，则 DFT 和 IDFT 可分别表示为 + +$$\begin{equation} +X=D_Nx,\ \ x=D_N^{-1}X +\end{equation}$$ + +其中 Vandermonde 矩阵 + +$$\begin{equation} +D_N=\left[\begin{matrix}1&1&1&\cdots&1\\1&W_N^1&W_N^2&\cdots&W_N^{\left(N-1\right)}\\1&W_N^2&W_N^4&\cdots&W_N^{2\left(N-1\right)}\\\vdots&\vdots&\vdots&&\vdots\\1&W_N^{\left(N-1\right)}&W_N^{2\left(N-1\right)}&\cdots&W_N^{\left(N-1\right)^2}\\\end{matrix}\right] +\end{equation}$$ + +$$\begin{equation} +D_N^{-1}=\frac{1}{N}\left[\begin{matrix}1&1&1&\cdots&1\\1&W_N^{-1}&W_N^{-2}&\cdots&W_N^{-\left(N-1\right)}\\1&W_N^{-2}&W_N^{-4}&\cdots&W_N^{-2\left(N-1\right)}\\\vdots&\vdots&\vdots&&\vdots\\1&W_N^{-\left(N-1\right)}&W_N^{-2\left(N-1\right)}&\cdots&W_N^{-\left(N-1\right)^2}\\\end{matrix}\right] +\end{equation}$$ + +且有 $D_N^{-1}=\dfrac{1}{N}D_N^\ast$。 + +# 3. 离散傅里叶变换（DFT）的性质 + +## 3.1 线性性质 + +$$\begin{equation} +\mathrm{DFT}\left[ax_1\left(n\right)+{bx}_2\left(n\right)\right]=aX_1\left(k\right)+bX_2\left(k\right) +\end{equation}$$ + +## 3.2 用正变换计算逆变换 + +$$\begin{equation} +x\left(n\right)=\frac{1}{N}\left[\sum_{k=0}^{N-1}{X^\ast\left(k\right)W_N^{nk}}\right]^\ast,\ \ n=0,1,\cdots,N-1 +\end{equation}$$ + +证明： + +$$\begin{equation} +\frac{1}{N}\left[\sum_{k=0}^{N-1}{X^\ast\left(k\right)W_N^{nk}}\right]^\ast=\frac{1}{N}\left[\sum_{k=0}^{N-1}{X\left(k\right)W_N^{-nk}}\right]=x\left(n\right) +\end{equation}$$ + +## 3.3 对称定理 + +$$\begin{equation} +\mathrm{DFT}\left[\frac{1}{N}X\left(n\right)\right]=x\left(-k\right)=x\left(N-k\right) +\end{equation}$$ + +证明： + +$$\begin{equation} +\begin{aligned} +&x\left(-n\right)=x\left(N-n\right)=\frac{1}{N}\sum_{k=0}^{N-1}{X\left(k\right)W_N^{-\left(N-n\right)k}}=\frac{1}{N}\sum_{k=0}^{N-1}{X\left(k\right)W_N^{nk}} +\\ +\Longrightarrow &x\left(-k\right)=x\left(N-k\right)=\sum_{n=0}^{N-1}{\left[\frac{X\left(n\right)}{N}\right]W_N^{kn}}=\mathrm{DFT}\left[\frac{X\left(n\right)}{N}\right] +\end{aligned} +\end{equation}$$ + +## 3.4 反转定理 + +$$\begin{equation} +\mathrm{DFT}\left[x\left(-n\right)\right]=X\left(-k\right) +\end{equation}$$ + +证明： + +$$\begin{equation} +\mathrm{DFT}\left[x\left(-n\right)\right]=\sum_{n=0}^{N-1}{x\left(-n\right)W_N^{nk}}\xlongequal{令m=-n}\sum_{m=-\left(N-1\right)}^{0}{x\left(m\right)W_N^{-mk}}=\sum_{m=0}^{N-1}{x\left(m\right)W_N^{m\cdot\left(-k\right)}}=X\left(-k\right) +\end{equation}$$ + +## 3.5 序列的总和 + +$$\begin{equation} +\sum_{n=0}^{N-1}x\left(n\right)=\left.\sum_{n=0}^{N-1}{x\left(n\right)W_N^{kn}}\right|_{k=0}\left.=X\left(k\right)\right|_{k=0}=X\left(0\right) +\end{equation}$$ + +## 3.6 序列的起始值 + +$$\begin{equation} +x\left(0\right)=\frac{1}{N}\left.\sum_{k=0}^{N-1}{X\left(k\right)W_N^{-kn}}\right|_{n=0}=\frac{1}{N}\sum_{k=0}^{N-1}X\left(k\right) +\end{equation}$$ + +## 3.7 延长序列的 DFT + +  将 $N$ 点长序列 $x\left(n\right)\left(n=1,2,\cdots,N-1\right)$ 补零延长 $rN$ 点得到序列 $g\left(n\right)$，即 + +$$\begin{equation} +g\left(n\right)= +\begin{aligned} +\begin{cases} +x(n),&n=0,1,⋯,N-1\\ +0,&n=N,N+1,⋯rN-1 +\end{cases} +\end{aligned} +\end{equation}$$ + +则 $g\left(n\right)$ 的 DFT + +$$\begin{equation} +\begin{aligned} +G\left(k\right)&=\mathrm{DFT}\left[g\left(n\right)\right]=\sum_{n=0}^{rN-1}{g\left(n\right)e^{-j\frac{2\pi nk}{rN}}}=\sum_{n=0}^{N-1}{g\left(n\right)e^{-j\frac{2\pi n\left(\frac{k}{r}\right)}{N}}} +\\ +&=X\left(\frac{k}{r}\right),\ \ k=0,1,\cdots,rN-1\ +\end{aligned} +\end{equation}$$ + +由此，**$G\left(k\right)$ 与 $X\left(k\right )$具有相同的形状，但 $G\left(k\right)$ 的频谱间隔比 $X\left(k\right)$ 的小。即 {% label primary @通过补零，频谱变得更加细致。但是补零不能提高频谱分辨能力。%}** + +## 3.8 有限长序列的圆周特性 + +### 3.8.1 圆周移位（周期化→移位→取主值） + +$$\begin{equation} +x\left(\left(n-m\right)\right)_NR_N\left(n\right)= +\begin{aligned} +\begin{cases} +x\left(n-m\right),&m\leqslant n\leqslant N-1 +\\ +x\left(N-m+n\right),&0\leqslant n\leqslant m +\end{cases} +\end{aligned} +\end{equation}$$ + +### 3.9.2 圆周反转（周期化→反转→取主值） + +$$\begin{equation} +x\left(\left(-n\right)\right)_NR_N\left(n\right)= +\begin{aligned} +\begin{cases} +x\left(0\right),&n=0 +\\ +x\left(N-n\right),&1\leqslant n\leqslant N-1 +\end{cases} +\end{aligned} +\end{equation}$$ + +### 3.8.3 有限长序列的时间圆周移位定理（可以通过 DFS 证明） + +$$\begin{equation} +\mathrm{DFT}\left[x\left(\left(n+m\right)\right)_NR_N\left(n\right)\right]=W_N^{-mk}X\left(k\right) +\end{equation}$$ + +由此，**有限长序列的圆周移位导致频谱线性相移，而对频谱幅度无影响**。 + +### 3.8.4 限长序列的频率圆周移位定理（调制特性，可以通过 DFS 证明） + +$$\begin{equation} +\mathrm{IDFT}\left[X\left(\left(k+l\right)\right)_NR_N\left(k\right)\right]=W_N^{nl}x\left(n\right)=e^{-j\frac{2\pi}{N}nl}x\left(n\right) +\end{equation}$$ + +由此，**时域序列的调制等效于频域的圆周移位**，且可以推得 + +$$\begin{equation} +\mathrm{DFT}\left[x\left(n\right)\cos{\left(\frac{2\pi nl}{N}\right)}\right]=\frac{1}{2}\left[X\left(\left(k-l\right)\right)_N+X\left(\left(k+l\right)\right)_N\right]R_N\left(k\right) +\end{equation}$$ + +$$\begin{equation} +\mathrm{DFT}\left[x\left(n\right)\sin{\left(\frac{2\pi nl}{N}\right)}\right]=\frac{1}{2}\left[X\left(\left(k-l\right)\right)_N-X\left(\left(k+l\right)\right)_N\right]R_N\left(k\right) +\end{equation}$$ + +### 3.8.5 DFT的圆周对称性 + +#### 3.8.5.1 奇序列和偶序列的 DFT + +奇序列的 DFT + +$$\begin{equation} +x\left(n\right)=-x\left(-n\right)=-x\left(\left(-n\right)\right)_NR_N\left(n\right)\rightarrow X\left(k\right)=-X\left(-k\right)=-X\left(\left(N-k\right)\right)_NR_N\left(k\right) +\end{equation}$$ + +偶序列的 DFT + +$$\begin{equation} +x\left(n\right)=x\left(-n\right)=x\left(\left(-n\right)\right)_NR_N\left(n\right)\rightarrow X\left(k\right)=X\left(-k\right)=X\left(\left(N-k\right)\right)_NR_N\left(k\right) +\end{equation}$$ + +#### 3.8.5.2 序列的对称分解 + +  普通序列 $x\left(n\right)$ 的共轭对称分量和共轭反对称分量分别为 + +$$\begin{equation} +x_e\left(n\right)=\frac{1}{2}\left[x\left(n\right)+x^\ast\left(-n\right)\right],\ \ x_o\left(n\right)=\frac{1}{2}\left[x\left(n\right)-x^\ast\left(-n\right)\right] +\end{equation}$$ + +对有限长序列的 DFT 进行分析时，定义周期共轭对称分量和周期共轭反对称分量分别为 + +$$\begin{equation} +x_{ep}\left(n\right)=\frac{1}{2}\left[x\left(\left(n\right)\right)_N+x^\ast\left(\left(-n\right)\right)_N\right]R_N\left(n\right)=\frac{1}{2}\left[x\left(n\right)+x^\ast\left(N-n\right)\right] +\end{equation}$$ + +$$\begin{equation} +x_{op}\left(n\right)=\frac{1}{2}\left[x\left(\left(n\right)\right)_N-x^\ast\left(\left(-n\right)\right)_N\right]R_N\left(n\right)=\frac{1}{2}\left[x\left(n\right)-x^\ast\left(N-n\right)\right] +\end{equation}$$ + +二者关系为 + +$$\begin{equation} +x_{ep}\left(n\right)=\left[x_e\left(n\right)+x_e\left(n-N\right)\right]R_N\left(n\right),\ \ x_{op}\left(n\right)=\left[x_o\left(n\right)+x_o\left(n-N\right)\right]R_N\left(n\right) +\end{equation}$$ + + + +#### 3.8.5.3 共轭复序列的 DFT + +$$\begin{equation} +\mathrm{DFT}\left[x^\ast\left(n\right)\right]=X^\ast\left(\left(N-k\right)\right)_N +\end{equation}$$ + +$$\begin{equation} +\mathrm{DFT}\left[x^\ast\left(-n\right)\right]=X^\ast\left(k\right) +\end{equation}$$ + +证明： + +$$\begin{equation} +\mathrm{DFT}\left[x^\ast\left(n\right)\right]=\sum_{n=0}^{N-1}{x^\ast\left(n\right)W_N^{kn}}=\left[\sum_{n=0}^{N-1}{x\left(n\right)W_N^{-kn}}\right]^\ast=X^\ast\left(-k\right)=X^\ast\left(\left(N-k\right)\right)_N +\end{equation}$$ + +$$\begin{equation} +\begin{aligned} +\mathrm{DFT}\left[x^\ast\left(-n\right)\right]&=\sum_{n=0}^{N-1}{x^\ast\left(-n\right)W_N^{kn}}=\left[\sum_{n=0}^{N-1}{x\left(-n\right)W_N^{-kn}}\right]^\ast +\\ +&=\left[\sum_{m=-N+1}^{0}{x\left(m\right)W_N^{km}}\right]^\ast=\left[\sum_{n=0}^{N-1}{x\left(n\right)W_N^{kn}}\right]^\ast=X^\ast\left(k\right) +\end{aligned} +\end{equation}$$ + + + +#### 3.8.5.4 复数序列的 DFT + +$$\begin{equation} +x\left(n\right)=x_r\left(n\right)+jx_i\left(n\right) +\end{equation}$$ + +$$\begin{equation} +\quad\quad\quad\updownarrow\qquad\quad \updownarrow +\end{equation}$$ + +$$\begin{equation} +X\left(k\right)=X_{ep}\left(k\right)+X_{op}\left(k\right) +\end{equation}$$ + +证明： + +$$\begin{equation} +X_{ep}\left(k\right)=\mathrm{DFT}\left[x_r\left(n\right)\right]=\frac{1}{2}\mathrm{DFT}\left[x\left(n\right)+x^\ast\left(n\right)\right]=\frac{1}{2}\left[X\left(k\right)+X^\ast\left(N-k\right)\right] +\end{equation}$$ + +$$\begin{equation} +X_{op}\left(k\right)=\mathrm{DFT}\left[{jx}_i\left(n\right)\right]=\frac{1}{2}\mathrm{DFT}\left[x\left(n\right)-x^\ast\left(n\right)\right]=\frac{1}{2}\left[X\left(k\right)+X^\ast\left(N-k\right)\right] +\end{equation}$$ + +且由 [3.8.5.3 共轭复序列的 DFT](#3.8.5.3) 可推知 + +$$\begin{equation} +X_{ep}\left(k\right)=X_{ep}^\ast\left(N-k\right)\ \ \ (实部相等，虚部相反) +\end{equation}$$ + +$$\begin{equation} +X_{op}\left(k\right)=-X_{op}^\ast\left(N-k\right)(实部相反，虚部相等) +\end{equation}$$ + +#### 3.8.5.5 DFT 的奇偶虚实特性 + +| $x(n)$ | $X(k)$ | +| :----: | :----------------: | +| 偶 | 偶 | +| 奇 | 奇 | +| 实 | 实部为偶，虚部为奇 | +| 虚 | 实部为奇，虚部为偶 | +| 实偶 | 实偶 | +| 实奇 | 虚奇 | +| 虚偶 | 虚偶 | +| 虚奇 | 实奇 | + +### 3.8.6 圆周卷积定理 + +> 注：此处用 $\odot$ 表示圆周卷积，不同教材中圆周卷积的符号表述可能不同。 + +#### 3.8.6.1 时域圆周卷积定理 + +$$\begin{equation} +\begin{aligned} +x\left(n\right)=x_1\left(n\right)\odot x_2\left(n\right)&=\widetilde{y}\left(n\right)=\sum_{m=0}^{N-1}{x_1\left(m\right)\left[x_2\left(\left(n-m\right)\right)_NR_N\left(n\right)\right]} +\\ +&=\mathrm{IDFT}\left[X_1\left(k\right)\cdot X_2\left(k\right)\right] +\end{aligned} +\end{equation}$$ + +#### 3.8.6.2 频域圆周卷积定理 + +$$\begin{equation} +\begin{aligned} +X\left(k\right)=X_1\left(k\right)\odot X_2\left(k\right)&=\widetilde{X}\left(k\right)R_N\left(k\right)=\frac{1}{N}\sum_{l=0}^{N-1}{X_1\left(l\right)\left[X_2\left(\left(k-l\right)\right)_NR_N\left(k\right)\right]} +\\ +&=\mathrm{DFT}\left[x_1\left(n\right)\cdot x_2\left(n\right)\right] +\end{aligned} +\end{equation}$$ + +#### 3.8.6.3 圆周卷积和线性卷积的关系 + +  长度为 $L$ 的序列 $x_1\left(n\right)$ 和长度为 $M$ 的序列 $x_2\left(n\right)$ 的周期卷积是二者的线性卷积的以 $N$ 为周期的周期延拓，若满足 + +$$\begin{equation} +N\geqslant L+M-1 +\end{equation}$$ + +则周期卷积的主值序列即圆周卷积与线性卷积完全相同，即 $x_1\left(n\right)\odot x_2\left(n\right)=x_1\left(n\right)\ast x_2\left(n\right)$ + +#### 3.8.6.4 圆周卷积的矩阵算法 + +$$\begin{equation} +y=h\odot x=\left[\begin{matrix}h_1&0&\cdots&0&h_m&\cdots\\h_2&h_1&\cdots&0&0&\cdots\\\vdots&\vdots&&\vdots&\vdots&\cdots\\h_m&h_{m-1}&\cdots&0&0&\cdots\\0&h_m&\cdots&h_1&0&\cdots\\0&0&&\vdots&h_1&\cdots\\\vdots&\vdots&&h_{m-1}&\vdots&\cdots\\0&0&\cdots&h_m&h_{m-1}&\cdots\\\end{matrix}\right]_{N\times N}\cdot\left[\begin{matrix}x_1\\x_2\\\vdots\\x_n\\0\\\vdots\\0\\\end{matrix}\right]_{N\times1} +\end{equation}$$ + +也可先计算线性卷积以 $N$ 为周期进行周期化后取主值得到圆周卷积。注意周期化时 $N 在 $X\left(z\right)$ 的单位圆上均匀采样 $N$ 点，其 IDFT 是序列 $x\left(n\right)$ 以周期 $N$ 进行周期延拓后的主值序列； +> +> 在 $X\left(e^{j\omega}\right)$ 的一个周期内均匀采样 $N$ 点，其 IDFT 是序列 $x\left(n\right)$ 以周期 $N$ 进行周期延拓后的主值序列； +> +> 在 $X\left(e^{j\omega}\right)$ 的一个周期内均匀采样 $M$ 点，其 IDFT 是序列 $x\left(n\right)$ 本身。 + +由此，对于列长为M的有限长序列 $x\left(n\right)$，频域取样不失真的条件是取样点数 $N$ 满足 + +$$\begin{equation} +N\geqslant M +\end{equation}$$ + +并且对于无限长序列 $x\left(n\right)$，无论 $N$ 取值如何，都不可能消除混叠，只能随着取样点 $N$ 增加而接近 $x\left(n\right)$。 + +# 5. 内插公式 + +  用 DFT 表示 ZT + +$$\begin{equation} +\begin{aligned} +X\left(z\right)&=\sum_{n=0}^{N-1}{x\left(n\right)z^{-n}}=\sum_{n=0}^{N-1}{\left[\frac{1}{N}\sum_{k=0}^{N-1}{X\left(k\right)W_N^{-nk}}\right]z^{-n}}=\frac{1}{N}\sum_{k=0}^{N-1}\left[X\left(k\right)\sum_{n=0}^{N-1}\left(W_N^{-k}z^{-1}\right)^n\right] +\\ +&=\frac{1}{N}\sum_{k=0}^{N-1}{X\left(k\right)\frac{1-W_N^{-kn}z^{-N}}{1-W_N^{-k}z^{-1}}}=\frac{1}{N}\sum_{k=0}^{N-1}{X\left(k\right)\frac{1-z^{-N}}{1-W_N^{-k}z^{-1}}} +\\ +&=\frac{1-z^{-N}}{N}\sum_{k=0}^{N-1}\frac{X\left(k\right)}{1-W_N^{-k}z^{-1}}=\sum_{k=0}^{N-1}{X\left(k\right)\varphi_k\left(z\right)} +\end{aligned} +\end{equation}$$ + +其中 $X\left(z\right)$ 的内插公式 + +$$\begin{equation} +\varphi_k\left(z\right)=\frac{1}{N}\cdot\frac{1-z^{-N}}{1-W_N^{-k}z^{-1}} +\end{equation}$$ + +  用 DFT 表示 DTFT + +$$\begin{equation} +X\left(e^{j\omega}\right)=\left.X\left(z\right)\right|_{z=e^{j\omega}}=\sum_{k=0}^{N-1}{X\left(k\right)\varphi_k\left(e^{j\omega}\right)}=\sum_{k=0}^{N-1}X\left(k\right)\varphi\left(\omega-k\frac{2\pi}{N}\right) +\end{equation}$$ + +其中内插函数 + +$$\begin{equation} +\varphi\left(\omega\right)=\frac{1}{N}\frac{\sin{\left(\frac{\omega N}{2}\right)}}{\sin{\left(\frac{\omega}{2}\right)}}e^{-j\omega\left(\frac{N-1}{2}\right)} +\end{equation}$$ + +**时域取样定理说明，对频域带限的信号，可对其进行时域取样而不丢失任何信息。频域取样理论说明，时间有限的信号（有限长序列），也可对其进行频域取样（DFT）而不丢失任何信息。** + +# 6. 用 DFT 对连续时间信号逼近的问题 + +## 6.1 混叠现象 + +  时域采样会引起频谱周期化，通过 $f_s\geqslant 2f_h$ 消除。 + +## 6.2 频谱泄露现象 + +  取用有限个数据，即将信号截断的过程，等同于将信号乘以窗函数。如果窗函数是矩形窗，在时域内将信号截断相当于乘以矩形窗，而在频域中则相当于原信号的频谱与矩形窗的频谱卷积，使得原信号的频谱被展宽，造成频谱泄露。**通常通过选择频谱特性更接近于 $\delta\left(k\right)$ 的函数来减少频谱泄露。** + +## 6.3 栅栏效应 + +  用 DFT 计算频谱时，采样点为基频的整数倍，由此观察到的频谱特征时有限的，**若在两个取样点之间有剧烈的频谱变化，将无法检测出来**。通常通过在原序列的末端添加一些零值点来改变周期内的点数，从而在保持原有频谱连续形式不变的情况下，变更了谱线的位置，原来看不到的频谱分量就能移动到可见的位置上。 + +## 6.4 DFT 参数的选择 + +### 6.4.1 取样频率 $f_s$（由取样定理） + +$$\begin{equation} +f_s\geqslant 2f_h +\end{equation}$$ + +### 6.4.2 取样周期 $T$ + +$$\begin{equation} +T=\frac{1}{f_s}\leqslant\frac{1}{2f_h} +\end{equation}$$ + +### 6.4.3 频率分辨率（频率分量间的增量）$F$ + +$$\begin{equation} +F=\frac{f_s}{N} +\end{equation}$$ + +$F$ 越小，频率分辨率越高。 + +### 6.4.4 最小记录长度（周期性函数的有效周期）$t_p$ + +$$\begin{equation} +t_p=\frac{1}{F}=NT +\end{equation}$$ + +### 6.4.5 一个记录长度内的点数 $N$ + +$$\begin{equation} +N\geqslant\frac{2f_h}{F} +\end{equation}$$ + +> 注：最高频率 $f_h$ 和频率分辨率 $F$ 存在矛盾，即 +> +> $$\begin{equation} +> f_h\uparrow \xrightarrow{f_s\geqslant 2f_h}f_s\uparrow \xrightarrow{T=\frac{1}{f_s}}T\downarrow \xrightarrow{N\text{给定时}}t_p\downarrow \longrightarrow F\uparrow \left( \text{频率分辨率}\downarrow \right) +> \end{equation}$$ + +# 参考文献 + +王世一. *数字信号处理, 修订版.* 北京理工大学出版社, 1997. diff --git a/source/_posts/2023-04-14-josh-dsp-part-4.md b/source/_posts/2023-04-14-josh-dsp-part-4.md new file mode 100644 index 0000000..3194562 --- /dev/null +++ b/source/_posts/2023-04-14-josh-dsp-part-4.md @@ -0,0 +1,524 @@ +--- +title: Josh 的学习笔记之数字信号处理（Part 4——快速傅里叶变换 FFT） +mathjax: true +comments: true +copyright: true +toc: + enable: true + number: false + wrap: true + expand_all: true + max_depth: 4 +abbrlink: 7c6d68a2 +date: 2023-04-14 23:11:52 +tags: + - 数字信号和处理 + - DSP + - 快速傅里叶变换 + - FFT + - 基2按时间抽取 + - 基2按频率抽取 + - 基4 + - 分裂基FFT + - 蝶形图 + - 旋转因子 +categories: + - - Josh 的学习笔记 + - 数字信号处理 +--- + +# 1. 直接计算 DFT 的计算量 + +  考察列长为 $N$ 的序列 $x\left(n\right)$ 的 DFT + +$$\begin{equation} +\begin{aligned} +X\left(k\right)&=\sum_{n=0}^{N-1}{x\left(n\right)W_N^{nk}}=\sum_{n=0}^{N-1}\left\{\Re{\left[x\left(n\right)\right]}+j\Im{\left[x\left(n\right)\right]}\right\}\left\{\Re{\left[W_N^{nk}\right]}+j\Im{\left[W_N^{nk}\right]}\right\} +\\ +&=\sum_{n=0}^{N-1}{\left\{ \Re \left[ x\left( n \right) \right] \Re \left[ W_{N}^{nk} \right] -\Im \left[ x\left( n \right) \right] \Im \left[ W_{N}^{nk} \right] \right\} +j\left\{ \Re \!\:\left[ x\left( n \right) \right] \Im \!\:\left[ W_{N}^{nk} \right] +\Im \left[ \!\:x\left( n \right) \right] \Re \!\:\left[ W_{N}^{nk} \right] \right\}} +\end{aligned} +\end{equation}$$ + +由此，直接计算 DFT 时，在不忽略 $W_N^0=1$ 等特例时，实数乘法次数为 $4N^2$，实数加法次数为 $2N\left(2N-1\right)$，均与 $N^2$ 成正比，计算量非常庞大。 + + + +# 2. 利用旋转因子 $W_N^{nk}$ 的特性改善 DFT 的运算效率 + +| | $W_N^{nk}=e^{j\frac{2\pi}{N}nk}$ | +| :--: | :--: | +| 对称性 | $\left(W_N^{nk}\right)^\ast=W_N^{-nk}=W_N^{\left(N-n\right)k}=W_N^{n\left(N-k\right)}$ | +| 周期性 | $W_N^{nk}=W_N^{\left(N+n\right)k}=W_N^{n\left(N+k\right)}$ | +| 可约性 | $W_N^{nk}=W_{mN}^{mnk},W_N^{nk}=W_{\frac{N}{m}}^{\frac{nk}{m}}$ | +| 特殊点 | $W_N^0=W_N^N=1,W_N^{\frac{N}{2}}=-1,W_N^{k+\frac{N}{2}}=-W_N^k$ | + +利用旋转因子 $W_N^{nk}$ 的特性，可将直接计算 DFT 算式中的对称项合并，从使长序列的DFT分解为为更小点数的DFT，即 + +$$\begin{equation} +\Re \left[ x\left( n \right) \right] \Re \left[ W_{N}^{nk} \right] +\Re \left[ x\left( N-n \right) \right] \Re \left[ W_{N}^{\left( N-n \right) k} \right] =\left\{ \Re \left[ x\left( n \right) \right] +\Re \left[ x\left( N-n \right) \right] \right\} \Re \left[ W_{N}^{nk} \right] +\end{equation}$$ + +$$\begin{equation} +-\Im \left[ x\left( n \right) \right] \Im \left[ W_{N}^{nk} \right] -\Im \left[ x\left( N-n \right) \right] \Im \left[ W_{N}^{\left( N-n \right) k} \right] =-\left\{ \Im \left[ x\left( n \right) \right] -\Im \left[x\left( N-n \right) \right] \right\} \Im \left[ W_{N}^{nk} \right] +\end{equation}$$ + +# 3. 基-2 按时间抽取的 FFT 算法 + +（radix-2 Decimation in Time Fast Fourier Transform, radix-2 DIT-FFT, i.e. Cooley-Turkey Method） + + + +## 3.1 算法原理 + +  设序列 $x\left(n\right)$ 的列长 $N=2^\nu$，$\nu$ 为整数，则其 DFT + +$$\begin{equation} +X\left(k\right)=\sum_{n=0}^{N-1}{x\left(n\right)W_N^{nk}},\ \ k=0,1,\cdots,N-1 +\end{equation}$$ + +**将 $x\left(n\right)$ 按照 $n$ 的奇偶分成两个子序列** + +$$\begin{equation} +\begin{cases} +x_1\left(r\right)=x\left(2r\right)\\ +x_2\left(r\right)=x\left(2r+1\right) +\end{cases} +,\ \ r=0,1,\cdots,\frac{N}{2}-1 +\end{equation}$$ + +则序列 DFT 的**前半部分**可表示为 + +$$\begin{equation} +\begin{aligned} +X\left(k\right)&=\sum_{r=0}^{\frac{N}{2}-1}{x\left(2r\right)W_N^{2rk}}+\sum_{r=0}^{\frac{N}{2}-1}{x\left(2r+1\right)W_N^{\left(2r+1\right)k}}=\sum_{r=0}^{\frac{N}{2}-1}{x_1\left(r\right)W_{\frac{N}{2}}^{rk}}+W_N^k\sum_{r=0}^{\frac{N}{2}-1}{x_2\left(r\right)W_{\frac{N}{2}}^{rk}} +\\ +&=\mathrm{DFT}\left[x_1\left(n\right)\right]+W_N^k\mathrm{DFT}\left[x_2\left(n\right)\right] +\\ +&=X_1\left(k\right)+W_N^kX_2\left(k\right),\ \ k=0,1,\cdots,\frac{N}{2}-1 +\end{aligned} +\end{equation}$$ + +由旋转因子的周期性，即 + +$$\begin{equation} +W_{\frac{N}{2}}^{rk}=W_{\frac{N}{2}}^{r\left(k+\frac{N}{2}\right)} +\end{equation}$$ + +可将序列 DFT 的**后半部分**表示为 + +$$\begin{equation} +\begin{aligned} +X\left(k+\frac{N}{2}\right)&=X_1\left(k+\frac{N}{2}\right)+W_N^{\left(k+\frac{N}{2}\right)}X_2\left(k+\frac{N}{2}\right) +\\ +&=X_1\left(k\right)+W_N^{\left(k+\frac{N}{2}\right)}X_2\left(k\right) +\\ +&=X_1\left(k\right)-W_N^kX_2\left(k\right),\ \ k=0,1,\cdots,\frac{N}{2}-1 +\end{aligned} +\end{equation}$$ + +类似地，可以将分离后的序列继续按奇偶序号分成两个序列计算，**直到每个序列中只含有两个点**。此时开始计算 DFT。这种方法称为**基-2 按时间抽取的 FFT 算法**。 + +## 3.2 蝶形运算流图符号 + +![基-2 按时间抽取的 FFT 算法的鲽形图](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-04-14-josh-dsp-part-4/2023-04-14-josh-dsp-part-4-010-Radix2DITFFTButterfly.png?imageMogr2/thumbnail/!100p|watermark/2/text/QEpvc2ggR2Fv/font/YWhyb25iZC50dGY=/fontsize/80/dissolve/20/gravity/southeast/dx/5/dy/5){width=1000px} + +## 3.3 8 点基-2 DIT-FFT 推导 + +![8 点 radix-2 DIT-FFT 推导](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-04-14-josh-dsp-part-4/2023-04-14-josh-dsp-part-4-020-Radix2DITFFTDerivation.png?imageMogr2/thumbnail/!100p|watermark/2/text/QEpvc2ggR2Fv/font/YWhyb25iZC50dGY=/fontsize/40/dissolve/20/gravity/southeast/dx/5/dy/5){width=1000px} + +  对于 8 点的序列，进行两次分解后每部分序列只含有两个点，可以直接计算 DFT。为了得到自然顺序输出的 DFT 序列，将第二级分解以后的结果作为输入。下面推导每级分解后序列的结果表示。 + +  对于第一级分解，直接用[算法原理](#3.1)中的式子表示结果，有 + +$$\begin{align} +&X\left(k\right)=X_1\left(k\right)+W_N^kX_2\left(k\right),&k=0,1,\cdots,\frac{N}{2}\\ +&X\left(k+\frac{N}{2}\right)=X_1\left(k\right)+W_N^kX_2\left(k\right),&k=0,1,\cdots,\frac{N}{2} +\end{align}$$ + +对于第二级分解 + +$$\begin{equation} +\begin{aligned} +X_1\left(k\right)&=\sum_{l=0}^{\frac{N}{4}-1}{x_1\left(2l\right)W_{\frac{N}{2}}^{2lk}}+\sum_{l=0}^{\frac{N}{4}-1}{x_1\left(2l+1\right)W_{\frac{N}{2}}^{\left(2l+1\right)k}} +\\ +&=\sum_{l=0}^{\frac{N}{4}-1}{x_3\left(l\right)W_{\frac{N}{4}}^{lk}}+W_{\frac{N}{2}}^k\sum_{l=0}^{\frac{N}{4}-1}{x_4\left(l\right)W_{\frac{N}{4}}^{lk}}=X_3\left(k\right)+W_{\frac{N}{2}}^kX_4\left(k\right),\ \ k=0,1,\cdots,\frac{N}{4}-1 +\end{aligned} +\end{equation}$$ + +$$\begin{equation} +X_1\left(k+\frac{N}{4}\right)=X_3\left(k\right)-W_{\frac{N}{2}}^kX_4\left(k\right),\ \ k=0,1,\cdots,\frac{N}{4}-1 +\end{equation}$$ + +同理 + +$$\begin{equation} +X_2\left(k\right)=X_5\left(k\right)+W_{\frac{N}{2}}^kX_6\left(k\right),X_2\left(k+\frac{N}{2}\right)=X_5\left(k\right)-W_{\frac{N}{2}}^kX_6\left(k\right),\ \ k=0,1,\cdots,\frac{N}{4}-1 +\end{equation}$$ + +对于最后剩下的四个 $\dfrac{N}{4}=2$ 点的 DFT，即 $X_3\left(k\right)、X_4\left(k\right)、X_5\left(k\right)、X_6\left(k\right)，k=0,1$，可按定义计算，如对 $X_3\left(k\right)$ + +$$\begin{equation} +\begin{aligned} +&X_3\left(k\right)=\sum_{l=0}^{\frac{N}{4}-1}{x_3\left(l\right)W_{\frac{N}{4}}^{lk}}=\sum_{l=0}^{1}{x_3\left(l\right)W_{\frac{N}{4}}^{lk}}\\ +\Longrightarrow &X_3\left(0\right)=x_3\left(0\right)+W_2^0x_3\left(1\right)=x\left(0\right)+W_N^0x\left(4\right),X_3\left(1\right)=x_3\left(0\right)+W_2^1x_3\left(1\right)=x\left(0\right)-W_N^0x\left(4\right) +\end{aligned} +\end{equation}$$ + +将上述各式的旋转因子统一为 $W_{\frac{N}{2}}^k=W_N^{2k}$，则一个 8 点的 DFT 可以分解为四个 $\dfrac{N}{4}$ 点的 DFT ，先做 $\dfrac{N}{4}$ 点的 DFT，再用相应的两个 $\dfrac{N}{4}$ 点 DFT 的结果合成 $\dfrac{N}{2}$ 点的 DFT，从而得到 $X_1\left(k\right)、X_2\left(k\right)$，最后组合成为 $N$ 点的 DFT。上述过程的运算流图如下左图，实际绘制按照如下右图即可。 + +![8 点 radix-2 DIT-FFT 的蝶形图](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-04-14-josh-dsp-part-4/2023-04-14-josh-dsp-part-4-030-8PointRadix2DITFFTButterfly.png?imageMogr2/thumbnail/!100p|watermark/2/text/QEpvc2ggR2Fv/font/YWhyb25iZC50dGY=/fontsize/80/dissolve/20/gravity/southeast/dx/5/dy/5){width=1000px} + +## 3.4 算法特点 + +### 3.4.1 原位运算 + +  每级计算所得点数相同，且上一级数据不影响下级运算，数据可以保存在原寄存器中。 + +### 3.4.2 比特反位序 + +  若下标 $n$ 可用二进制表示为 $n=\left(b_2b_1b_0\right)_2$，则 FFT 的输入样本序号 $k$ 就可用二进制表示为 $k=\left(b_0b_1b_2\right)_2$。 + +### 3.4.3 计算复杂度 + +  对于长度为 $N=2^\nu$ 的序列，共有 $\nu$ 级蝶形运算，每级有 $\dfrac{N}{2}$ 个蝶形，每个蝶形有 $1$ 次复数加法、$2$ 次复数乘法，因此 $\nu$ 级蝶形运算共有 + +$$\begin{equation} +\begin{aligned} +&复数乘法：m_F=\frac{N}{2}\cdot\nu=\frac{N}{2}\log_2{N}\\ +&复数加法：a_F=N\cdot\nu=N\log_2{N} +\end{aligned} +\end{equation}$$ + +### 3.4.4 各类蝶形运算两个点相距的“距离”及 $W_N^k$ 的变化规律 + +| 点数 $N$ | 列数 | 蝶形运算的种类数 | 蝶形运算的系数 | 参加蝶形运算的
两数据点的间距 +| :--: | :--: | :--: | :--: | :--: | +| $\top$ | 第 $1$ 列 | $1$ 种 | $W_2^0/W_8^0=1$ | $1$ | +| $8$ | 第 $2$ 列 | $2$ 种 | $W_4^0/W_8^0,{W_4^1/W}_8^2$ | $2$ | +| $\bot$ | 第 $3$ 列 | $4$ 种 | $W_8^0,W_8^1,W_8^2,W_8^3$ | $4$ | +| $\top$ | 第 $1$ 列 | $2^0$ 种 | $W_2^0/W_N^0$ | $1$ | +| $\mid$ | 第 $2$ 列 | $2^1$ 种 | $W_4^0/W_N^0,{W_4^1/W}_N^{\frac{N}{4}}$ | $2$ | +| $\mid$ | 第 $3$ 列 | $2^2$ 种 | $W_8^0,W_8^1,W_8^2,W_8^3/$
$W_N^0,W_N^{\frac{8}{N}},W_N^{2\cdot\frac{8}{N}},W_8^{3\cdot\frac{8}{N}}$ | $4$ | +| $2^\nu$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | +| $\mid$ | 第 $\nu-1$ 列 | $2^{\nu-2}$ 种 | $W_{\frac{N}{2}}^0,W_{\frac{N}{2}}^1,W_{\frac{N}{2}}^2,\cdots,W_{\frac{N}{2}}^{\frac{N}{4}-1}$ $W_N^0,W_N^2,W_N^4,\cdots,W_N^{\frac{N}{2}-2}$ | $\dfrac{N}{4}$ | +| $\bot$ | 第 $\nu$ 列 | $2^{\nu-1}$ 种 | $W_N^0,W_N^1,W_N^2,\cdots,W_N^{\frac{N}{2}-1}$ | $\dfrac{N}{2}$ | + +  以上述 8 点 FFT 为例，各列蝶形运算的种类数和系数以及参加蝶形运算的两数据点的间距如上表上半部。可见，每列的蝶形类型比前一列增加一倍，参加蝶形运算的两个数据点的间距也增大一倍。最后一列系数用的最多，为 4 个，即 $W_8^0,W_8^1,W_8^2,W_8^3$；而前一列只用带它偶序号的那一半，即 ${W_4^0/W}_8^0,{W_4^1/W}_8^2$；第一列只有一个系数即 ${W_2^0/W}_8^0$。 + +  上述结论可推广至 $N=2^\nu$ 的一般情况。如上表下半部。从前向后，第一列只有一种类型的蝶形运算 $W_N^0$，以后的每列蝶形的类型逗逼前一列增加一倍，到第 $\nu$ 列是 $\dfrac{N}{2}$ 个蝶形类型，系数是 $W_N^0,W_N^1,\cdots,W_N^{\frac{N}{2}-1}$，共 $\dfrac{N}{2}$ 个。从后向前，系数是后一级的偶数序号的那一半。参加蝶形运算的两个数据点的间距，是最末一级最大，其值为 $\dfrac{N}{2}$，向前每推进一列，间距减少一半。 + +# 4. 基-2 按频率抽取的 FFT 算法 + +（radix-2 Decimation in Frequency Fast Fourier Transform, radix-2 DIF-FFT, i.e. Sande-Tukey Method） + +## 4.1 算法原理 + +  设序列 $x\left(n\right)的列长N=2^\nu$，$\nu$ 为整数，则其 DFT + +$$\begin{equation} +X\left(k\right)=\sum_{n=0}^{N-1}{x\left(n\right)W_N^{nk}},\ \ k=0,1,\cdots,N-1 +\end{equation}$$ + +**将 $x\left(n\right)$ 按照 $n$ 的顺序分成两个子序列（即可将 $X\left(k\right)$ 按照奇偶分组）** + +$$\begin{equation} +\begin{cases} +前半子序列x\left(n\right)\\ +后半子序列x\left(n+\dfrac{N}{2}\right) +\end{cases} +,\ \ n=0,1,\cdots,\frac{N}{2}-1 +\end{equation}$$ + +则序列的 DFT 可表示为 + +$$\begin{equation} +\begin{aligned} +X\left(k\right)&=\sum_{n=0}^{\frac{N}{2}-1}{x\left(n\right)W_N^{nk}}+\sum_{n=\frac{N}{2}}^{N-1}{x\left(n\right)W_N^{nk}}=\sum_{n=0}^{\frac{N}{2}-1}{x\left(n\right)W_N^{nk}}+\sum_{n=0}^{\frac{N}{2}-1}{x\left(n+\frac{N}{2}\right)W_N^{(n+\frac{N}{2})k}} +\\ +&=\sum_{n=0}^{\frac{N}{2}-1}{\left[x\left(n\right)+W_N^{\frac{N}{2}k}x\left(n+\frac{N}{2}\right)\right]W_N^{nk}} +\\ +&\xlongequal{W_{N}^{\frac{N}{2}k}=\left( -1 \right) ^k}\sum_{n=0}^{\frac{N}{2}-1}{\left[x\left(n\right)+\left(-1\right)^kx\left(n+\frac{N}{2}\right)\right]W_N^{nk}},\ \ k=0,1,\cdots,\frac{N}{2}-1 +\end{aligned} +\end{equation}$$ + +由 $W_N^{\frac{N}{2}k}=\left(-1\right)^k$，按照 $k$ 的奇偶可将 $X\left(k\right)$ 分成两部分，令 $k=2r$ 及 $k=2r+1,\ r=0,1,\cdots,\dfrac{N}{2}-1$，则 + +$$\begin{align} +&X\left(2r\right)=\sum_{n=0}^{\frac{N}{2}-1}{\left[x\left(n\right)+x\left(n+\frac{N}{2}\right)\right]W_N^{2rn}}=\sum_{n=0}^{\frac{N}{2}-1}{\left[x\left(n\right)+x\left(n+\frac{N}{2}\right)\right]W_{\frac{N}{2}}^{rn}} +\\ +&X\left(2r+1\right)=\sum_{n=0}^{\frac{N}{2}-1}{\left[x\left(n\right)-x\left(n+\frac{N}{2}\right)\right]W_N^{2rn}}=\sum_{n=0}^{\frac{N}{2}-1}{\left[x\left(n\right)-x\left(n+\frac{N}{2}\right)\right]W_N^nW_{\frac{N}{2}}^{rn}} +\end{align}$$ + +上两式种的前者是输入序列前一半和后一半之和的 $\dfrac{N}{2}$ 点 DFT，后者是输入序列前一半和后一半之差与 $W_N^n$ 之积的 $\dfrac{N}{2}$ 点 DFT。令 + +$$\begin{equation} +\begin{cases} +x_1\left(n\right)=x\left(n\right)+x\left(n+\dfrac{N}{2}\right)\\ +x_2\left(n\right)=\left[x\left(n\right)-x\left(n+\dfrac{N}{2}\right)\right]W_N^n +\end{cases} +,\ \ n=0,1,\cdots,\frac{N}{2}-1 +\end{equation}$$ + +则 + +$$\begin{equation} +X\left(2r\right)=\sum_{n=0}^{\frac{N}{2}-1}{x_1\left(n\right)W_{\frac{N}{2}}^{rn}},\ \ X\left(2r+1\right)=\sum_{n=0}^{\frac{N}{2}-1}{x_2\left(n\right)W_{\frac{N}{2}}^{rn}} +\end{equation}$$ + +由此，可将 $N$ 点 DFT 按频率 $k$ 的奇偶分解为两个序列的 $\dfrac{N}{2}$ 点 DFT。类似地，可以将分离后的序列继续按频率k的奇偶分成两个序列计算，**直到每个序列中只含有两个点**。此时开始计算 DFT。这种方法称为**基-2 按频率抽取的 FFT 算法**。 + +## 4.2 蝶形运算流图符号 + +![基-2 按频率抽取的 FFT 算法的蝶形图](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-04-14-josh-dsp-part-4/2023-04-14-josh-dsp-part-4-040-Radix2DIFFFTButterfly.png?imageMogr2/thumbnail/!100p|watermark/2/text/QEpvc2ggR2Fv/font/YWhyb25iZC50dGY=/fontsize/80/dissolve/20/gravity/southeast/dx/5/dy/5){width=1000px} + +## 4.3 8 点基-2 DIF-FFT 推导简述 + +![8 点 radix-2 DIF-FFT 推导](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-04-14-josh-dsp-part-4/2023-04-14-josh-dsp-part-4-050-Radix2DIFFFTDerivation.png?imageMogr2/thumbnail/!100p|watermark/2/text/QEpvc2ggR2Fv/font/YWhyb25iZC50dGY=/fontsize/40/dissolve/20/gravity/southeast/dx/5/dy/5){width=1000px} + +运算流图如下左图，实际绘制按照如下右图即可。 + +![8 点 radix-2 DIF-FFT 的蝶形图](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-04-14-josh-dsp-part-4/2023-04-14-josh-dsp-part-4-060-8PointRadix2DIFFFTButterfly.png?imageMogr2/thumbnail/!100p|watermark/2/text/QEpvc2ggR2Fv/font/YWhyb25iZC50dGY=/fontsize/80/dissolve/20/gravity/southeast/dx/5/dy/5){width=1000px} + +**对于任何流图，只要保持各节点所连的支路及其传输系数不变，则不论节点位置怎么排列，所得流图总是等效的，最后所得 DFT 结果也是正确的**，只是数据的提取和存放次序不同而已。 + +## 4.4 算法特点 + +### 4.4.1 原位运算 + +  每级计算所得点数相同，且上一级数据不影响下级运算，数据可以保存在原寄存器中。 + +### 4.4.2 比特反位序 + +  若下标 $n$ 可用二进制表示为 $n=\left(b_2b_1b_0\right)_2$，则 FFT 的输入样本序号 $k$ 就可用二进制表示为 $k=\left(b_0b_1b_2\right)_2$。 + +### 4.4.3 计算复杂度 + +  对于长度为 $N=2^\nu$ 的序列，共有 $\nu$ 级蝶形运算，每级有 $\dfrac{N}{2}$ 个蝶形，每个蝶形有 $1$ 次复数加法、$2$ 次复数乘法，因此 $\nu$ 级蝶形运算共有 + +$$\begin{equation} +\begin{aligned} +&复数乘法：m_F=\frac{N}{2}\cdot\nu=\frac{N}{2}\log_2{N}\\ +&复数加法：a_F=N\cdot\nu=N\log_2{N} +\end{aligned} +\end{equation}$$ + +# 5. 基-2 DIT-FFT 和 基-2 DIF-FFT 的比较 + +![radix-2 DIT-FFT 和 radix-2 DIF-FFT 的比较](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-04-14-josh-dsp-part-4/2023-04-14-josh-dsp-part-4-070-Radix2DIFFFTRadix2DIFFFTComparison.png?imageMogr2/thumbnail/!100p|watermark/2/text/QEpvc2ggR2Fv/font/YWhyb25iZC50dGY=/fontsize/40/dissolve/20/gravity/southeast/dx/5/dy/5){width=1000px} + +# 6. 分裂基 FFT 算法（split-radix FFT Method） + +## 6.1 算法原理（基-4 DIF-FFT同理） + +  设序列 $x\left(n\right)$ 的列长 $N=2^\nu$，$\nu$为整数，则其 DFT + +$$\begin{equation} +X\left(k\right)=\sum_{n=0}^{N-1}{x\left(n\right)W_N^{nk}},\ \ k=0,1,\cdots,N-1 +\end{equation}$$ + +将 $x\left(n\right)$ 按照 $n$ 的顺序分成四个子序列 + +$$\begin{equation} +\begin{cases} +第一子序列x\left(n\right)\\ +第二子序列x\left(n+\dfrac{N}{4}\right)\\ +第三子序列x\left(n+\dfrac{N}{2}\right)\\ +第四子序列x\left(n+\dfrac{3N}{4}\right) +\end{cases} +,\ \ n=0,1,\cdots,\frac{N}{2}-1 +\end{equation}$$ + +则序列的 DFT 可表示为 + +$$\begin{equation} +\begin{aligned} +X\left(k\right)&=\sum_{n=0}^{\frac{N}{4}-1}{x\left(n\right)W_N^{nk}}+\sum_{n=\frac{N}{4}}^{\frac{N}{2}-1}{x\left(n\right)W_N^{nk}}+\sum_{n=\frac{N}{2}}^{\frac{3N}{4}-1}{x\left(n\right)W_N^{nk}}+\sum_{n=\frac{3N}{4}}^{N-1}{x\left(n\right)W_N^{nk}} +\\ +&=\sum_{n=0}^{\frac{N}{4}-1}\left[x\left(n\right)W_N^{nk}+x\left(n+\frac{N}{4}\right)W_N^{\left(n+\frac{N}{4}\right)k}+x\left(n+\frac{N}{2}\right)W_N^{\left(n+\frac{N}{2}\right)k}+x\left(n+\frac{3N}{4}\right)W_N^{\left(n+\frac{3N}{4}\right)k}\right] +\\ +&=\sum_{n=0}^{\frac{N}{4}-1}{\left[x\left(n\right)+x\left(n+\frac{N}{4}\right)W_N^{\frac{N}{4}k}+x\left(n+\frac{N}{2}\right)W_N^{\frac{N}{2}k}+x\left(n+\frac{3N}{4}\right)W_N^{\frac{3N}{4}k}\right]W_N^{nk}} +\\ +&=\sum_{n=0}^{\frac{N}{4}-1}{\left[x\left(n\right)+x\left(n+\frac{N}{4}\right)W_4^k+x\left(n+\frac{N}{2}\right)W_4^{2k}+x\left(n+\frac{3N}{4}\right)W_4^{3k}\right]W_N^{nk}} +\end{aligned} +\end{equation}$$ + +按照 $k$ 除以 $4$ 的余数可将 $X\left(k\right)$ 分成四部分，分别令 $k=4r、k=4r+1、k=4r+2、k=4r+3$，其中 $r=0,1,\cdots,\dfrac{N}{4}-1$，则有 + +$$\begin{equation} +\begin{cases} +X\left(4r\right)=\displaystyle\sum_{n=0}^{\frac{N}{4}-1}{\left[x\left(n\right)+x\left(n+\dfrac{N}{4}\right)+x\left(n+\dfrac{N}{2}\right)+x\left(n+\dfrac{3N}{4}\right)\right]W_{\frac{N}{4}}^{rn}} +\\ +X\left(4r+1\right)=\displaystyle\sum_{n=0}^{\frac{N}{4}-1}{\left[x\left(n\right)-jx\left(n+\dfrac{N}{4}\right)-x\left(n+\dfrac{N}{2}\right)+jx\left(n+\dfrac{3N}{4}\right)\right]W_N^nW_{\frac{N}{4}}^{rn}} +\\ +X\left(4r+2\right)=\displaystyle\sum_{n=0}^{\frac{N}{4}-1}{\left[x\left(n\right)-x\left(n+\dfrac{N}{4}\right)+x\left(n+\dfrac{N}{2}\right)-x\left(n+\dfrac{3N}{4}\right)\right]W_N^{2n}W_{\frac{N}{4}}^{rn}} +\\ +X\left(4r+3\right)=\displaystyle\sum_{n=0}^{\frac{N}{4}-1}{\left[x\left(n\right)+jx\left(n+\dfrac{N}{4}\right)-x\left(n+\dfrac{N}{2}\right)-jx\left(n+\dfrac{3N}{4}\right)\right]W_N^{3n}W_{\frac{N}{4}}^{rn}} +\end{cases} +,\ \ r=0,1,\cdots,\frac{N}{4}-1 +\end{equation}$$ + +上述过程也是基-4 DIF-FFT 的推导。上式中的偶数序号项合并，得 + +$$\begin{equation} +\begin{cases} +X\left(2r\right)=\displaystyle\sum_{n=0}^{\frac{N}{2}-1}{\left[x\left(n\right)+x\left(n+\dfrac{N}{2}\right)\right]W_N^{2rn}}, &r=0,1,\cdots,\dfrac{N}{2}-1 +\\ +X\left(4r+1\right)=\displaystyle\sum_{n=0}^{\frac{N}{4}-1}{\left[\left(x\left(n\right)-x\left(n+\dfrac{N}{2}\right)\right)-j\left(x\left(n+\dfrac{N}{4}\right)-x\left(n+\dfrac{3N}{4}\right)\right)\right]W_N^nW_{\frac{N}{4}}^{rn}},&r=0,1,\cdots,\dfrac{N}{4}-1 +\\ +X\left(4r+3\right)=\displaystyle\sum_{n=0}^{\frac{N}{4}-1}{\left[\left(x\left(n\right)-x\left(n+\dfrac{N}{2}\right)\right)+j\left(x\left(n+\dfrac{N}{4}\right)-x\left(n+\dfrac{3N}{4}\right)\right)\right]W_N^{3n}W_{\frac{N}{4}}^{rn}},&r=0,1,\cdots,\dfrac{N}{4}-1 +\end{cases} +\end{equation}$$ + +令 + +$$\begin{equation} +\begin{cases} +x_2\left(n\right)=x\left(n\right)+x\left(n+\dfrac{N}{2}\right), &n=0,1,\cdots,\dfrac{N}{2}-1 +\\ +x_4^1\left(n\right)=\left[\left(x\left(n\right)-x\left(n+\dfrac{N}{2}\right)\right)-j\left(x\left(n+\dfrac{N}{4}\right)-x\left(n+\dfrac{3N}{4}\right)\right)\right]W_N^n,&n=0,1,\cdots,\dfrac{N}{4}-1 +\\ +x_4^2\left(n\right)=\left[\left(x\left(n\right)-x\left(n+\dfrac{N}{2}\right)\right)+j\left(x\left(n+\dfrac{N}{4}\right)-x\left(n+\dfrac{3N}{4}\right)\right)\right]W_N^{3n},&n=0,1,\cdots,\dfrac{N}{4}-1 +\end{cases} +\end{equation}$$ + +则有 + +$$\begin{equation} +\begin{cases} +X\left(2r\right)=\displaystyle\sum_{n=0}^{\frac{N}{2}-1}{x_2\left(n\right)W_{\frac{N}{2}}^{rn}}=DFT\left[x_2\left(n\right)\right], &r=0,1,\cdots,\dfrac{N}{2}-1 +\\ +X\left(4r+1\right)=\displaystyle\sum_{n=0}^{\frac{N}{4}-1}{x_4^1\left(n\right)W_{\frac{N}{4}}^{rn}}=\mathrm{DFT}\left[x_4^1\left(n\right)\right],&r=0,1,\cdots,\dfrac{N}{4}-1 +\\ +X\left(4r+3\right)=\displaystyle\sum_{n=0}^{\frac{N}{4}-1}{x_4^2\left(n\right)W_{\frac{N}{4}}^{rn}}=\mathrm{DFT}\left[x_4^2\left(n\right)\right],&r=0,1,\cdots,\dfrac{N}{4}-1 +\end{cases} +\end{equation}$$ + +由此，DFT 的偶数部分采用基-2 算法，奇数部分采用基-4 算法，将一个 $N$ 点 DFT 的分解为一个 $\dfrac{N}{2}$ 点 DFT 和两个 $\dfrac{N}{4}$ 点 DFT。可将分裂基-2/4 算法看成**时域分段、频域抽取**的方法。先将整个 $N$ 点序列在时域分成两段，对于偶数段，在时域进一步分成两段，相当于在频域抽取偶数点 $X\left(2r\right)$；对于奇数段，在时域分成 4 段，相当于在频域抽取奇数点 $X\left(4r+1\right)$ 和 $X\left(4r+3\right)$。类似地，可将分离后的序列继续分解，直到最后一级。这种方法称为**分裂基 FFT 算法**。 + +## 6.2 L 蝶形运算流图符号 + +![L 蝶形运算流图符号](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-04-14-josh-dsp-part-4/2023-04-14-josh-dsp-part-4-080-Radix24LButterfly.png?imageMogr2/thumbnail/!100p|watermark/2/text/QEpvc2ggR2Fv/font/YWhyb25iZC50dGY=/fontsize/80/dissolve/20/gravity/southeast/dx/5/dy/5){width=1000px} + +## 6.3 8 点 radix-2/4 FFT举例 + +![8 点 radix-2/4 FFT举例](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-04-14-josh-dsp-part-4/2023-04-14-josh-dsp-part-4-090-8PointRadix24LButterfly.png?imageMogr2/thumbnail/!100p|watermark/2/text/QEpvc2ggR2Fv/font/YWhyb25iZC50dGY=/fontsize/80/dissolve/20/gravity/southeast/dx/5/dy/5){width=600px} + +## 6.4 运算量分析 + +  对于 $N=2^\nu$ 点 DFT，设第 $i$ 级有 $l_i$ 个蝶形，$j=0,1,\cdots,\nu-1$，则 + +$$\begin{equation} +\begin{cases} +l_1=\dfrac{N}{4}\\ +l_i=\dfrac{N}{4}-\dfrac{l_{j-1}}{2}=\dfrac{N}{6}\left[1-\left(-\dfrac{1}{2}\right)^i\right],\ \ i\neq1 +\end{cases} +\Longrightarrow\sum_{i=1}^{\nu-1}l_i=\frac{N}{6}\left[\nu-\frac{2}{3}+\frac{2}{3}\left(-\frac{1}{2}\right)^\nu\right] +\end{equation}$$ + +其中 $l_i$ 来自 L 蝶形的特殊结构。总的复乘次数为 L 蝶形数的两倍，而总的复加次数与基-2 FFT 算法相同，即 + +$$\begin{equation} +\begin{aligned} +&复数乘法：C_F=2\cdot\sum_{i=1}^{\nu-1}l_i=\frac{N}{3}\log_2{N}-\frac{2N}{9}+\frac{\left(-1\right)^\nu2}{9}\ \ \propto\ \ \frac{N}{3}\log_2{N}\\ +&复数加法：a_F=N\cdot\nu=N\log_2{N} +\end{aligned} +\end{equation}$$ + +显然，**复乘次数相比基-2 FFT 算法下降了33%**。 + +> **L 蝶形数与流图中 $-j$ 的个数相等。** + +# 7. 实序列的 FFT 算法 + +## 7.1 用一个 $N$ 点 DFT 同时计算两个 $N$ 点实序列的 DFT + +  设 $x_1\left(n\right),x_2\left(n\right)$ 是相互独立的两个 $N$ 点实序列，将 $x_1\left(n\right)$ 和 $x_2\left(n\right)$ 分别当作一复序列的实部和虚部，即 + +$$\begin{equation} +x\left(n\right)=x_1\left(n\right)+jx_2\left(n\right) +\end{equation}$$ + +则由 DFT 的线性性质，$x\left(n\right)$ 的 DFT 可表示为 + +$$\begin{equation} +X\left(k\right)=X_1\left(k\right)+jX_2\left(k\right) +\end{equation}$$ + +再由的[复序列的 DFT 性质](https://josh-gao.top/posts/fd9f9a29.html#3.8.5.3) + +$$\begin{equation} +X\left(k\right)=X_{ep}\left(k\right)+X_{op}\left(k\right) +\end{equation}$$ + +则有 $X_1\left(k\right)$ 和 $X_2\left(k\right)$ + +$$\begin{equation} +X_1\left(k\right)=X_{ep}\left(k\right)=\frac{1}{2}\left[X\left(k\right)+X^\ast\left(N-k\right)\right] +\end{equation}$$ + +$$\begin{equation} +X_2\left(k\right)=-{jX}_{op}\left(k\right)=-\frac{j}{2}\left[X\left(k\right)-X^\ast\left(N-k\right)\right] +\end{equation}$$ + +蝶形图如下： +![用一个 N点 DFT 同时计算两个 N 点实序列的 DFT的蝶形图](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-04-14-josh-dsp-part-4/2023-04-14-josh-dsp-part-4-100-ComputeNPointComplxDFT.png?imageMogr2/thumbnail/!100p|watermark/2/text/QEpvc2ggR2Fv/font/YWhyb25iZC50dGY=/fontsize/40/dissolve/20/gravity/southeast/dx/5/dy/5){width=500px} + +## 7.2 用一个 $N$ 点 DFT 计算一个 $2N$ 点实序列的 DFT + +  设 $x\left(n\right)$ 是一个 $2N$ 点的实序列，将其按照奇偶分为两个独立的 $N$ 点实序列 $x_1\left(n\right),x_2\left(n\right)$，即 + +$$\begin{equation} +\begin{cases} +x_1\left(n\right)=x\left(2n\right)\\ +x_2\left(n\right)=x\left(2n+1\right) +\end{cases} +,\ \ n=0,1,\cdots,N-1 +\end{equation}$$ + +将 $x_1\left(n\right)$ 和 $x_2\left(n\right)$ 分别当作一复序列的实部和虚部，即 + +$$\begin{equation} +y\left(n\right)=x_1\left(n\right)+jx_2\left(n\right) +\end{equation}$$ + +则由 DFT 的线性性质，$y\left(n\right)$ 的 DFT 可表示为 + +$$\begin{equation} +Y\left(k\right)=X_1\left(k\right)+jX_2\left(k\right) +\end{equation}$$ + +则由上述讨论，有 $X_1\left(k\right)$ 和 $X_2\left(k\right)$ + +$$\begin{equation} +X_1\left(k\right)=Y_{ep}\left(k\right)=\frac{1}{2}\left[Y\left(k\right)+Y^\ast\left(N-k\right)\right] +\end{equation}$$ + +$$\begin{equation} +X_2\left(k\right)=-{jY}_{op}\left(k\right)=-\frac{j}{2}\left[Y\left(k\right)-Y^\ast\left(N-k\right)\right] +\end{equation}$$ + +因此 $x\left(n\right)$ 的傅里叶变换（类似基-2 DIT） + +$$\begin{equation} +\begin{aligned} +X\left(k\right)&=\sum_{n=0}^{2N-1}{x\left(n\right)W_N^{nk}}=\sum_{n=0}^{N-1}{x\left(2n\right)W_{2N}^{2nk}}+\sum_{n=0}^{N-1}{x\left(2n+1\right)W_{2N}^{\left(2n+1\right)k}} +\\ +&=\sum_{n=0}^{N-1}{x\left(2n\right)W_{2N}^{2nk}}+W_{2N}^k\sum_{n=0}^{N-1}{x\left(2n+1\right)W_{2N}^{2nk}} +\\ +&=\sum_{n=0}^{N-1}{x_1\left(n\right)W_N^{nk}}+W_{2N}^k\sum_{n=0}^{N-1}{x_2\left(n\right)W_N^{nk}} +\\ +&=X_1\left(k\right)+W_{2N}^kX_2\left(k\right),\ \ \ \ \ \ \ \ \ \ \ k=0,1,\cdots,N-1 +\end{aligned} +\end{equation}$$ + +$$\begin{equation} +X\left(k+N\right)=X_1\left(k\right)-W_{2N}^kX_2\left(k\right),\ \ k=0,1,\cdots,N-1\qquad\qquad\quad\ +\end{equation}$$ + +蝶形图如下： + +![用一个 N 点 DFT 计算一个 2N 点实序列的 DFT的蝶形图](https://josh-blog-1257563604.cos.ap-beijing.myqcloud.com/img/2023-04-14-josh-dsp-part-4/2023-04-14-josh-dsp-part-4-110-Compute2NPointRealDFT.png?imageMogr2/thumbnail/!100p|watermark/2/text/QEpvc2ggR2Fv/font/YWhyb25iZC50dGY=/fontsize/40/dissolve/20/gravity/southeast/dx/5/dy/5){width=500px} + +# 参考文献 + +王世一. *数字信号处理, 修订版.* 北京理工大学出版社, 1997.