19.6.4 调和分析

以公式或经验给出的以 $2\pi$ 为周期的周期函数 $f\left(x\right)$ ,应以三角多项式或傅里叶和逼近:

$$\begin{array}{}\text{(19.210)}& g\left(x\right)=\frac{a_0}{2}+\sum _{k=1}^n\left(a_k\cos kx+b_k\sin kx\right),\end{array}$$

其中系数 $a_0,a_k,b_k$ 是未知实数. 确定这些系数正是调和分析的课题.

19.6.4.1 三角插值公式

1. 傅里叶系数公式

因为函数系 $1,\cos kx,\sin kx\left(k=1,2,\cdots ,n\right)$ 在区间 $\left[0,2\pi \right]$ 内关于权函数 $\omega =$ 1 是正交的, 根据 (19.172) 应用连续最小二乘法得到系数公式

$$a_k=\frac{1}{\pi }{\int }_0^{2\pi }f\left(x\right)\cos kx\mathrm{d}x,b_k=\frac{1}{\pi }{\int }_0^{2\pi }f\left(x\right)\sin kx\mathrm{d}x\left(k=0,1,2,\cdots ,n\right).$$

(19.211)

由公式 (19.211) 计算得到的系数 $a_k,b_k$ 称为周期函数 $f\left(x\right)$ 的傅里叶系数 (参见第 633 页 7.4).

若 (19.211) 中的积分复杂或函数 $f\left(x\right)$ 仅在离散点上已知,则傅里叶系数只可由数值积分近似确定.

使用有 $N+1$ 个等距节点的梯形公式 (参见第 1253 页 19.3.2.2)

$$\begin{array}{}\text{(19.212)}& x_{\nu }=\nu h\left(\nu =0,1,\cdots ,N\right),h=\frac{2\pi }{N}\end{array}$$

得近似公式

$$a_k\approx {\tilde{a}}_k=\frac{2}{N}\sum _{\nu =1}^{N}f\left(x_{\nu }\right)\cos k{x}_{\nu },b_k\approx {\tilde{b}}_k=\frac{2}{N}\sum _{\nu =1}^{N}f\left(x_{\nu }\right)\sin k{x}_{\nu }\left(k=0,1,2,\cdots ,n\right).$$

(19.213)

在周期函数的情况下梯形公式变为非常简单的矩形公式. 如下事实使它有更高的精度: 若 $f\left(x\right)$ 是周期函数且 $2m+2$ 次可微,则梯形公式的误差阶为 $O\left(h^{2m+2}\right)$ .

2. 三角插值

以 ${\tilde{a}}_k,{\tilde{b}}_k$ 为近似系数的某些特殊的三角多项式有重要的性质. 这里介绍其中的两个性质.

(1) 插值 假设 $N=2n$ 成立. 系数由 (19.123) 给出的特殊的三角多项式

$$\begin{array}{}\text{(19.214)}& {\tilde{g}}_1\left(x\right)=\frac{1}{2}{\tilde{a}}_0+\sum _{k=1}^{n-1}\left({\tilde{a}}_k\cos kx+{\tilde{b}}_k\sin kx\right)+\frac{1}{2}{\tilde{a}}_n\cos nx\end{array}$$

在插值点 $x_v\left(19.212\right)$ 满足插值条件

$$\begin{array}{}\text{(19.215)}& {\tilde{g}}_1\left(x_{\nu }\right)=f\left(x_{\nu }\right)\left(\nu =1,2,\cdots ,N\right).\end{array}$$

由 $f\left(x\right)$ 的性质,成立 $f\left(x_0\right)=f\left(x_N\right)$ .

(2)平均近似假设 $N=2n$ 成立. 特殊的三角多项式

$$\begin{array}{}\text{(19.216)}& {\tilde{g}}_2\left(x\right)=\frac{1}{2}{\tilde{a}}_0+\sum _{k=1}^m\left({\tilde{a}}_k\cos kx+{\tilde{b}}_k\sin kx\right),\end{array}$$

对于 $m<n$ ,系数 (19.123) 按离散二次平均关于 $N$ 个节点 $x_{\nu }$ (19.212) 逼近函数 $f\left(x\right)$ ,即残量平方和

$$\begin{array}{}\text{(19.217)}& F=\sum _{\nu =1}^N{[f\left(x_{\nu }\right)-{\tilde{g}}_2\left(x_{\nu }\right)]}^2\end{array}$$

取得极小. 公式 (19.213) 正是用不同方法有效计算傅里叶系数的出发点.

19.6.4.2 快速傅里叶变换 (FFT)

1. 计算傅里叶系数的计算量

公式 (19.213) 中的求和也与离散傅里叶变换相关联, 例如在电工学、脉冲和图像处理中. 这里 $N$ 可能非常大,因为计算傅里叶系数的 $N$ 个近似值 (19.213) 需要大约 $N^2$ 个加法和乘法运算,所以必须以合理的方法计算求和.

对于 $N=2^p$ 的特殊情况,借助所谓快速傅里叶变换 (FFT),乘法运算的计算量可由 $N^2\left(=2^{2p}\right)$ 减少为 $pN\left(=p{2}^p\right)$ . 该减少量可以由右边的例子说明.(19.218)该方法如此有效地减少计算量和计算时间, 使得在重要的应用领域小型计算机也够用了.

	$N^2$	$pN$
10	$\sim {10}^6$	$\sim {10}^4$

FFT 将 $N$ 单位根,即方程 $z^N=1$ 的解的性质应用于 (19.213) 中的逐次求和.

2. 傅里叶和的复数表达

若在傅里叶和 (19.210) 中代入如下公式:

$$\begin{array}{}\text{(19.219)}& \cos kx=\frac{1}{2}\left({\mathrm{e}}^{\mathrm{i}kx}+{\mathrm{e}}^{-\mathrm{i}kx}\right),\sin kx=\frac{\mathrm{i}}{2}\left({\mathrm{e}}^{-\mathrm{i}kx}-{\mathrm{e}}^{\mathrm{i}kx}\right),\end{array}$$

则 FFT 的原理可以相当简单地表述为复数形式

$$g\left(x\right)=\frac{1}{2}{a}_0+\sum _{k=1}^n\left(a_k\cos kx+b_k\sin kx\right)$$$$\begin{array}{}\text{(19.220)}& =\frac{1}{2}{a}_0+\sum _{k=1}^n\left(\frac{a_k-\mathrm{i}{b}_k}{2}{\mathrm{e}}^{\mathrm{i}kx}+\frac{a_k+\mathrm{i}{b}_k}{2}{\mathrm{e}}^{-\mathrm{i}kx}\right).\end{array}$$

通过代换

$$\begin{array}{}\text{(19.221a)}& c_k=\frac{a_k-\mathrm{i}{b}_k}{2},\end{array}$$

由 (19.211), 有

$$\begin{array}{}\text{(19.221b)}& c_k=\frac{1}{2\pi }{\int }_0^{2\pi }f\left(x\right){\mathrm{e}}^{-\mathrm{i}kx}\mathrm{d}x,\end{array}$$

故 (19.220) 化为傅里叶和的复数表达式

$$\begin{array}{}\text{(19.222)}& g\left(x\right)=\sum _{k=-n}^n{c}_k{\mathrm{e}}^{\mathrm{i}kx}\text{,其中 }{c}_{-k}={\overline{c}}_k.\end{array}$$

若复系数 $c_k$ 已知,则要求的实傅里叶系数可以由下面简单的方法得到

$$\begin{array}{}\text{(19.223)}& a_0=2c_0,a_k=2\mathrm{Re}\left(c_k\right),b_k=-2\mathrm{Im}\left(c_k\right)\left(k=1,2,\cdots ,n\right).\end{array}$$

3. 复傅里叶系数的数值计算

为数值确定 $c_k$ ,类似于 (19.212) 和 (19.213),可对 (19.221b) 应用梯形公式,从而得到离散复傅里叶系数 ${\tilde{c}}_k$ :

$$\begin{array}{}\text{(19.224a)}& {\tilde{c}}_k=\frac{1}{N}\sum _{\nu =0}^{N-1}f\left(x_{\nu }\right){\mathrm{e}}^{-\mathrm{i}k{x}_{\nu }}=\sum _{\nu =0}^{N-1}{f}_{\nu }{\omega }_N^{k\nu }\left(k=0,1,2,\cdots ,n\right),\end{array}$$$$\begin{array}{}\text{(19.224b)}& f_{\nu }=\frac{1}{N}f\left(x_{\nu }\right),x_{\nu }=\frac{2\pi \nu }{N}\left(\nu =0,1,2,\cdots ,N-1\right),{\omega }_N={\mathrm{e}}^{-\frac{2\pi \mathrm{i}}{N}}.\end{array}$$

关系式 (19.224a) 和 (19.224b) 称为数值 $f_{\nu }\left(\nu =0,1,2,\cdots ,N-1\right)$ 的长度为 $N$ 的离散复傅里叶变换.

指数 $w_N^{\nu }=z\left(\nu =0,1,2,\cdots ,N-1\right)$ 满足方程 $z^N=1$ ,称之为 $N$ 次单位根. 因 ${\mathrm{e}}^{-2\pi \mathrm{i}}=1$ ,故

$$\begin{array}{}\text{(19.225)}& {\omega }_N^N=1,{\omega }_N^{N+1}={\omega }_N^1,{\omega }_N^{N+2}={\omega }_N^2,\cdots .\end{array}$$

长度为 $N=2n$ 的离散复傅里叶变换可以按下列方式化为两个长度为 $\frac{N}{2}=n$ 的变换, 利用这一事实可有效计算 (19.224a):

a) 对于每个偶数指标的系数 ${\tilde{c}}_k$ ,即 $k=2l$ ,成立

$${\tilde{c}}_{2l}=\sum _{\nu =0}^{2n-1}{f}_{\nu }{\omega }_N^{2l\nu }=\sum _{\nu =0}^{n-1}\left[f_{\nu }{\omega }_N^{2l\nu }+f_{n+\nu }{\omega }_N^{2l\left(n+\nu \right)}\right]=\sum _{\nu =0}^{n-1}\left[f_{\nu }+f_{n+\nu }\right]{\omega }_N^{2l\nu },$$

(19.226)

这里用到等式 $w_N^{2l\left(n+\nu \right)}=w_N^{l2n}{w}_N^{2l\nu }=w_N^{2l\nu }$ .

代入

$$\begin{array}{}\text{(19.227)}& y_{\nu }=f_{\nu }+f_{n+\nu }\left(\nu =0,1,2,\cdots ,n-1\right),\end{array}$$

并考虑到 $w_N^2=w_n$ ,和式

$$\begin{array}{}\text{(19.228)}& {\tilde{c}}_{2l}=\sum _{\nu =0}^{n-1}{y}_{\nu }{\omega }_n^{l\nu }\left(\nu =0,1,2,\cdots ,n-1\right)\end{array}$$

便是以 $\frac{N}{2}$ 为长度的数值 $y_{\nu }\left(\nu =0,1,2,\cdots ,n-1\right)$ 的离散复傅里叶变换.

b) 对每个奇数指标的系数 ${\tilde{c}}_k$ ,即 $k=2l+1$ ,类似可得

$$\begin{array}{}\text{(19.229)}& {\tilde{c}}_{2l+1}=\sum _{\nu =0}^{2n-1}{f}_{\nu }{\omega }_N^{\left(2l+1\right)\nu }=\sum _{\nu =0}^{n-1}\left[\left(f_{\nu }-f_{n+\nu }\right){\omega }_N^{\nu }\right]{\omega }_N^{2l\nu },\end{array}$$

代入

$$\begin{array}{}\text{(19.230)}& y_{n+\nu }=\left(f_{\nu }-f_{n+\nu }\right){\omega }_N^{\nu }\left(\nu =0,1,2,\cdots ,n-1\right)\end{array}$$

并考虑到 $w_N^2=w_n$ ,则和式

$$\begin{array}{}\text{(19.231)}& {\tilde{c}}_{2l+1}=\sum _{\nu =0}^{n-1}{y}_{n+\nu }{\omega }_n^{l\nu }\left(\nu =0,1,2,\cdots ,n-1\right)\end{array}$$

是以 $\frac{N}{2}$ 为长度的数值 $y_{n+\nu }\left(\nu =0,1,2,\cdots ,n-1\right)$ 的离散复傅里叶变换.

若 $N$ 为 2 的幂次,即 $N=2^p$ ( $p$ 为自然数),则根据 a) 和 b) 的归化,即将离散复傅里叶变换化为两个一半长度的离散复傅里叶变换, 是可以继续下去的. 应用 $p$ 次归化后则称 FFT.

根据 (19.230) 每步归化要求 $N/2$ 次复乘法运算, FFT 方法的计算量为

$$\begin{array}{}\text{(19.232)}& \frac{N}{2}p=\frac{N}{2}{\log }_{2}N\end{array}$$

4.FFT 的格式

对于特殊情况 $N=8=2^3$ ,根据 (19.227) 和 (19.230), FFT 的三个相应的归化步骤由如下格式 1 说明.

格式 1

	第 1 步	第 2 步	第 3 步
$f_0$	$y_0=f_0+f_4$	$y_0\text{:=}{y}_0+y_2$	$y_0\text{:=}{y}_0+y_1$ $y_1\text{:=}\left(y_0-y_1\right){\omega }_2^0$ $y_2\text{:=}{y}_2+y_3$ $y_3\text{:=}\left(y_2-y_3\right){\omega }_2^0$ $y_4\text{:=}{y}_4+y_5$ $y_5\text{:=}\left(y_4-y_5\right){\omega }_2^0$ $y_6\text{:=}{y}_6+y_7$ $y_7\text{:=}\left(y_6-y_7\right){\omega }_2^0$	$={\tilde{c}}_0$
$f_1$	$y_1=f_1+f_5$	$y_1\text{:=}{y}_1+y_3$		$={\tilde{c}}_4$
$f_2$	$y_2=f_2+f_6$	$y_2\text{:=}\left(y_0-y_2\right){\omega }_4^0$		$={\tilde{c}}_2$
$f_3$	$y_3=f_3+f_7$	$y_3\text{:=}\left(y_1-y_3\right){\omega }_4^1$		$={\tilde{c}}_6$
$f_4$	$y_4=\left(f_0-f_4\right){\omega }_8^0$	$y_4\text{:=}{y}_4+y_6$		$={\tilde{c}}_1$
$f_5$	$y_5=\left(f_1-f_5\right){\omega }_8^1$	$y_5\text{:=}{y}_5+y_7$		$={\tilde{c}}_5$
$f_6$	$y_6=\left(f_2-f_6\right){\omega }_8^2$	$y_6\text{:=}\left(y_4-y_6\right){\omega }_4^0$		$={\tilde{c}}_3$
$f_7$	$y_7=\left(f_3-f_7\right){\omega }_8^3$	$y_7\text{:=}\left(y_5-y_7\right){\omega }_4^1$		$={\tilde{c}}_7$
	$N=8,n\text{:=}4$ , ${\omega }_8={\mathrm{e}}^{-\frac{2\pi \mathrm{i}}{8}}$	$N\text{:=}4,n\text{:=}2,$ ${\omega }_4={\omega }_8^2$	$N\text{:=}2,n\text{:=}1$ , ${\omega }_2={\omega }_4^2$

可以注意到偶数及奇数指标项是如何出现的. 在格式 2 中 (19.233) 阐述了该方法的结构.

格式 2

$$\begin{array}{}\text{(19.233)}& {\tilde{c}}_k\Rightarrow \{\begin{array}{l}{\tilde{c}}_{2k}\Rightarrow \{\begin{array}{c}{\tilde{c}}_{2k}\Rightarrow \{\begin{array}{c}{\tilde{c}}_{8k}\\ {\tilde{c}}_{8k+4}\\ {\tilde{c}}_{8k+2}\\ {\tilde{c}}_{8k+2}\\ {\tilde{c}}_{8k+6}\end{array}\\ {\tilde{c}}_{8k+6}\end{array}\\ {\tilde{c}}_{2k+1}\Rightarrow \{\begin{array}{c}{\tilde{c}}_{8k+1}\Rightarrow \{\begin{array}{c}{\tilde{c}}_{8k+1}\\ {\tilde{c}}_{8k+5}\\ {\tilde{c}}_{8k+5}\\ {\tilde{c}}_{8k+6}\end{array}\\ {\tilde{c}}_{8k+3}\end{array}\end{array}\end{array}$$$$\left(k=0,1,\cdots ,7\right)\left(k=0,1,2,3\right)\left(k=0,1\right)\left(k=0\right).$$

若将系数 ${\tilde{c}}_k$ 代入格式 1,并在第 1 步前和第 3 步后考虑指标的二进制形式,则易知通过简单反转二进制形式指标的位数顺序, 可得要求的系数的阶. 这显示在格式 3 中.

格式 3

标志		第 1 步	第 2 步	第 3 步	标志
${\tilde{c}}_0$	000	${\tilde{c}}_0$	${\tilde{c}}_0$	${\tilde{c}}_0$	000
${\tilde{c}}_1$	00L	${\tilde{c}}_2$	${\tilde{c}}_4$	${\tilde{c}}_4$	L00
${\tilde{c}}_2$	0L0	${\tilde{c}}_4$	${\tilde{c}}_2$	${\tilde{c}}_2$	0L0
${\tilde{c}}_3$	0LL	${\tilde{c}}_6$	${\tilde{c}}_6$	${\tilde{c}}_6$	LL0
${\tilde{c}}_4$	L00	${\tilde{c}}_1$	${\tilde{c}}_1$	${\tilde{c}}_1$	00L
${\tilde{c}}_5$	L0L	${\tilde{c}}_3$	${\tilde{c}}_5$	${\tilde{c}}_5$	L0L
${\tilde{c}}_6$	LLO	${\tilde{c}}_5$	${\tilde{c}}_3$	${\tilde{c}}_3$	0LL
$\tilde{C}7$	LLL	$\tilde{C}7$	C7	$\tilde{C}7$	LLL

考虑以 $2\pi$ 为周期的函数 $f\left(x\right)=\{\begin{array}{ll}2{\pi }^2,& x=0,\\ x^2,& 0<x<2\pi ,\end{array}$ 将 FFT 用于离散傅里叶变换. 选取 $N=8,x_{\nu }=\frac{2\pi }{8},f_{\nu }=\frac{1}{8}f\left(x_{\nu }\right)\left(\nu =0,1,2,\cdots ,7\right),w_8={\mathrm{e}}^{-\frac{2\pi \mathrm{i}}{8}}=$ $0.707107\left(1-\mathrm{i}\right),w_8^2=-\mathrm{i},w_8^3=-0.707107\left(1+\mathrm{i}\right)$ . 得到格式 4 .

格式 4

	第 1 步	第 2 步	第 3 步
$f_0=2.467401$	$y_0=3.701102$	$y_0=6.785353$	$y_0=13.262281={\tilde{c}}_0$
$f_1=0.077106$	$y_1=2.004763$	$y_1=6.476928$	$y_1=0.308425={\tilde{c}}_4$
$f_2=0.308425$	$y_2=3.084251$	$y_2=0.616851$	$y_2=0.616851$ $+2.467402\mathrm{i}={\tilde{c}}_2$
$f_3=0.693957$	$y_3=4.472165$	$y_3=2.467402\mathrm{i}$	$y_3=0.616851$ $-2.467402\mathrm{i}={\tilde{c}}_6$
$f_4=1.233701$	$y_4=1.233700$	$y_4=1.233700$ +2.467401 i	$y_4=2.106058$ $+5.956833\mathrm{i}={\tilde{c}}_1$
$f_5=1.927657$	$y_5=-1.308537\left(1-\mathrm{i}\right)$	$y_5=0.872358$ +3.489432 i	$y_5=0.361342$ $-1.022031\mathrm{i}={\tilde{c}}_5$
$f_6=2.775826$	$y_6=2.467401\mathrm{i}$	$y_6=1.233700$	$y_6=0.361342$
		-2.467401 i	+1.022031 i $={\tilde{c}}_3$
$f_7=3.778208$	$y_7=2.180895\left(1+\mathrm{i}\right)$	$y_7=-0.872358$ +3.489432 i	$y_7=2.106058$ $-5.956833\mathrm{i}={\tilde{c}}_7$

通过三步归化, 根据 (19.233) 得到要求的实傅里叶系数 (见 (19.234)).

$$a_0=26.524562$$$$a_1=4.212116b_1=-11.913666$$$$\begin{array}{}\text{(19.234)}& a_2=1.233702b_2=-4.934804\end{array}$$$$a_3=0.722684b_3=-2.044062$$$$a_4=0.616850b_4=0$$

在该例中, 可注意到离散复傅里叶系数的一般性质

$${\tilde{c}}_{N-k}={\tilde{c}}_k.$$

当 $k=1,2,3$ 时,即有 ${\tilde{c}}_7={\tilde{c}}_1,{\tilde{c}}_6={\tilde{c}}_2,{\tilde{c}}_5={\tilde{c}}_3$ .