Fourier series, Fourier transform, DFT, FFT

Table of Contents

1 傅里叶级数和变换简介

法国数字家傅里叶(Jean-Baptiste Joseph Fourier)被世人铭记的最大贡献记载在他于1807年发表的传记中和1822年出版“Théorie analytique de la chaleur(热分析理论)”的一书中。傅里叶在这个领域的贡献是, 他指出任何“周期函数”都可以表示为不同频率的正弦和(或)余弦之和的形式,每个正弦项和(或)余弦项乘以不同的系数。现在,该和称为傅里叶级数(Fourier series)。 这个概念出现后,在当时遭到了怀疑,因为它不直观。

甚至,非周期函数(但该曲线下的面积是有限的)也可以用“正弦和(或)余弦乘以加权函数”的积分来表示。这种情况下的公式就是傅里叶变换(Fourier transform),其作用在多数理论和应用学科中甚至远大于傅里叶级数。

参考:信号与系统(第二版,郑君里)

1.1 数学基础知识回顾——复数

后文会用到复数,这里回顾一下关于复数的基本知识。

复数 \(C\) 的定义如下:
\[C = R + jI\]
其中, \(R\) 和 \(I\) 是实数, \(R\) 表示复数的实部, \(I\) 表示复数的虚部。 \(j\) 是虚数单位,即 \(j= \sqrt{-1}\) (注:数字上常用 \(i\) 表示虚数单位,这里用 \(j\) 是按照电工学中通常的习惯)。
一个复数 \(C\) 的共轭(Complex conjugate)表示为 \(C^{*}\) ,其定义为:
\[C^{*} = R - jI\]
复数从几何的角度可以被看成是平面(称为复平面)上的一个点,其横坐标是实轴( \(R\) 的值),其纵坐标是虚轴( \(I\) 的值)。也就是说,复数 \(R + jI\) 是复平面直角坐标系统中的点 \((R,I)\) 。
有时,在极坐标下表示复数很有用:
\[C = |C| (\cos \theta + j \sin \theta)\]
其中, \(|C|=\sqrt{R^2 + I^2}\) 是复平面的原点到点 \((R,I)\) 的向量的长度, \(\theta\) 是该向量与实轴的夹角。
利用欧拉公式(Euler’s formula):
\[e^{j\theta} = \cos \theta + j \sin \theta\]
极坐标下复数又可以表示(称为指数形式的表示)为:
\[C = |C| e^{j \theta}\]
例如,复数 \(1 + j2\) 的极坐标表示是 \(\sqrt{5} e^{j \theta}\) ,其中 \(\theta=64.4^{\circ}\) 。

前面的公式还可以用于复函数。例如,变量 \(u\) 的复函数 \(F(u)\) 可表示为 \(F(u) = R(u) + jI(u)\) ,其中 \(R(u)\) 和 \(I(u)\) 分别是实分量函数和虚分量函数。复共轭函数是 \(F^{*}(u) = R(u) - jI(u)\) 。

2 傅里叶级数

2.1 傅里叶级数基本定义

法国数字家傅里叶发现, 具有周期 \(T\) 的连续变量 \(t\) 的周期函数 \(f(t)\) 可以被描述为乘以适当系数的正弦和余弦和。 这个和就是傅里叶级数,它具有如下形式:

\begin{align*} f(t) &= a_0 + a_1 \cos(\frac{2 \pi}{T} t) + b_1 \sin(\frac{2 \pi}{T} t) + a_2 \cos(2 \frac{2 \pi}{T} t) + b_2 \sin(2 \frac{2 \pi}{T} t) + \cdots + a_n \cos(n \frac{2 \pi}{T} t) + b_n \sin(n \frac{2 \pi}{T} t) \\ &= a_0 + \sum_{n=1}^{\infty} \left[a_n \cos(n \frac{2 \pi}{T} t) + b_n \sin (n \frac{2 \pi}{T} t) \right] , \qquad (n = 1,2,\cdots) \tag{傅里叶级数的三角形式} \\ \end{align*}

其中,各个系数可按以下各式计算。
直流分量:
\[a_0 = \frac{1}{T} \int_{t_0}^{t_0 + T}f(t) \, \mathrm{d}t\]
余弦分量的幅度:
\[a_n = \frac{2}{T} \int_{t_0}^{t_0 + T}f(t) \cos(n \frac{2 \pi}{T} t) \, \mathrm{d}t , \qquad (n = 1,2,\cdots)\]
正弦分量的幅度:
\[b_n = \frac{2}{T} \int_{t_0}^{t_0 + T}f(t) \sin(n \frac{2 \pi}{T} t) \, \mathrm{d}t , \qquad (n = 1,2,\cdots)\]
为方便起见,通常积分区间 \(t_0 \sim t_0 + T\) 取为 \(0 \sim T\) 或者 \(-\frac{T}{2} \sim \frac{T}{2}\) 。

说明:并非任意周期信号都能进行傅里叶级数展开,需要满足一定的条件才能进行傅里叶级数展开。不过,在实践中,这些条件一般都会满足,这点不介绍这些条件。

2.1.1 实例:求周期锯齿脉冲信号的傅里叶级数展开

求周期函数 \(s(x)\) (如图 1 所示,它是周期锯齿脉冲信号)的傅里叶级数表示。

fourier_sawtooth_wave.png

Figure 1: Sawtooth wave, \(s(x) = x / \pi\) on the interval \((-\pi ,\pi ]\)

函数 \(s(x)\) 的周期为 \(T=2 \pi\) ,傅里叶级数中直流分量,余弦分量的幅度,正弦分量的幅度分别计算如下。
\[\begin{aligned} a_0 &= \frac{1}{2 \pi} \int_{-\pi}^{\pi} s(x) \, \mathrm{d}x = 0 \\ a_n &= \frac{2}{2 \pi} \int_{- \pi}^{\pi} s(x) \cos(n x) \, \mathrm{d}x = 0 , \qquad (n = 1,2,\cdots) \\ b_n &= \frac{2}{2 \pi} \int_{- \pi}^{\pi} s(x) \sin(n x) \, \mathrm{d}x = \frac{2 (-1)^{n+1}}{\pi n}, \qquad (n = 1,2,\cdots) \\ \end{aligned}\]
把上面系数代入到傅里叶级数展开式中,可得:
\[\begin{aligned} s(x) &= a_0 + \sum_{n=1}^{\infty} \left[a_n \cos(n \frac{2 \pi}{T} x) + b_n \sin (n \frac{2 \pi}{T} x) \right] \\ &= \frac{2}{\pi} \sum_{n=1}^{\infty} \frac{(-1)^{n+1}}{n} \sin(nx) \\ &= \frac{2}{\pi} \left[ \sin(x) - \frac{1}{2} \sin(2x) + \frac{1}{3} \sin(3x) - \frac{1}{4} \sin(x) + \cdots \right] \end{aligned}\]

由上式知,傅里叶级数展开后,要无限多项才能完全逼近原函数。实际应用中,经常采用有限项级数来代替无限项级数。如,取前1到5个正弦项,如图 2 所示。

fourier_sawtooth_wave_partial_fourier_series.gif

Figure 2: 取前1到5个正弦项来近似Sawtooth wave

2.1.1.1 函数对称性与傅里叶系数的关系

通过前面例子,可知周期锯齿脉冲信号(如图 1 所示,它是奇函数)仅含正弦分量。这不是巧合,利用“正弦函数为奇函数,余弦函数为偶函数”以及奇偶函数的性质,我们容易推导出更一般的结论:
奇函数的傅里叶级数展开后,直流分量和余弦分量为零,只含正弦分量。
偶函数的傅里叶级数展开后,正弦分量为零,只含直流分量和余弦分量。

2.2 傅里叶级数的另一种常见三角形式(幅度谱,相位谱)

如果将前面介绍的傅里叶级数中“同频率项合并”,傅里叶级数还可以写为另一种形式:

\begin{equation} f(t) = c_0 + \sum_{n=1}^{\infty} c_n \cos(n \frac{2 \pi}{T} t + \varphi_n) \tag{傅里叶级数的另一种常见三角形式} \end{equation}

或者写为:
\[f(t) = d_0 + \sum_{n=1}^{\infty} d_n \sin(n \frac{2 \pi}{T} t + \theta_n)\]
各个系数为( \(a_0, a_n, b_n, n = 1,2,\cdots\) 的计算公式前面已经介绍过):
\[\begin{aligned} & c_0 = d_0 = a_0 \\ & c_n = d_n = \sqrt{a_n^2 + b_n^2} , \qquad (n = 1,2,\cdots) \\ & \tan \varphi_n = - \frac{b_n}{a_n} \\ & \tan \theta_n = \frac{a_n}{b_n} \\ \end{aligned}\]

通常把频率为 \(1/T\) (此时角频率为 \(\frac{2 \pi}{T}\) , \(n=1\) )的分量称为“基波”,把频率为 \(2/T\) (角频率为 \(\frac{4 \pi}{T}\) )、 \(3/T\) (角频率为 \(\frac{6 \pi}{T}\) )等的分量称为“二次谐波”,“三次谐波”等。

显然,各分量幅度( \(c_n\) 或 \(d_n\) )及相位( \(\varphi_n\) 或 \(\theta_n\) )确定下来后,周期函数的傅里叶级数展开就确定了。

我们把幅度 \(c_n\) 和角频率 \(\frac{2 \pi n}{T}\) 之间的关系图称为幅度频谱(简称幅度谱),把相位 \(\varphi_n\) 和角频率 \(\frac{2 \pi n}{T}\) 之间的关系图称为相位频谱(简称相位谱)。

fourier_magnitude_and_phase.jpg

Figure 3: 某周期函数的傅里叶“幅度谱”和“相位谱”

不过,很多时候,相位谱不是很重要,常常直接忽略它。

2.3 傅里叶级数的指数形式

设周期函数 \(f(t)\) 的周期为 \(T\) ,记其角频率 \(\omega_1 = 2 \pi / T\) ,下面将从傅里叶级数的基本定义推导出傅里叶级数的指数形式。
由欧拉公式可知 \(\cos(n \omega_1 t) = \frac{1}{2} (e^{j n \omega_1 t} + e^{- j n \omega_1 t}), \, \sin(n \omega_1 t) = \frac{1}{2j} (e^{j n \omega_1 t} - e^{- j n \omega_1 t})\) ,把这两个式子代入到前面介绍的傅里叶级数基本定义中,有:
\[\begin{aligned} f(t) &= a_0 + \sum_{n=1}^{\infty} \left[a_n \cos(n \omega_1 t) + b_n \sin (n \omega_1 t) \right] \\ &= a_0 + \sum_{n=1}^{\infty} \left( a_n \frac{e^{j n \omega_1 t} + e^{- j n \omega_1 t}}{2} + b_n \frac{e^{j n \omega_1 t} - e^{- j n \omega_1 t}}{2j} \right) \\ &= a_0 + \sum_{n=1}^{\infty} \left( \frac{1}{2} (a_n - j b_n) e^{j n \omega_1 t} + \frac{1}{2} (a_n + j b_n) e^{- j n \omega_1 t} \right) \\ \end{aligned}\]

令 \(F(n \omega_1) = \frac{1}{2} (a_n - j b_n)\) ,由于 \(a_n\) 是 \(n\) 的偶函数, \(b_n\) 是 \(n\) 的奇函数,可知 \(F(- n \omega_1) = \frac{1}{2} (a_n + j b_n)\) ,从而上式可接着变形为:
\[f(t) = a_0 + \sum_{n=1}^{\infty} \left( F(n \omega_1) e^{j n \omega_1 t} + F (- n \omega_1) e^{- j n \omega_1 t} \right)\]
令 \(F(0) = a_0\) ,考虑到: \(\sum_{n=1}^{\infty} F (- n \omega_1) e^{- j n \omega_1 t} = \sum_{n=-1}^{-\infty} F (n \omega_1) e^{j n \omega_1 t}\) ,从而上式可接着变形为:
\[f(t) = \sum_{n=- \infty}^{\infty}F(n \omega_1) e^{j \omega_1 t}\]
一般我们把 \(F(n \omega_1)\) 简写为 \(F_n\) ,上式即为 周期函数 \(f(t)\) (设周期为 \(T\) )的傅里叶级数的指数形式 :

\begin{equation} f(t) = \sum_{n=- \infty}^{\infty}F_n e^{j \omega_1 t} \tag{傅里叶级数的指数形式} \end{equation}

其中,
\[\begin{aligned} \omega_1 & = \frac{2 \pi}{T} \\ F_n & = \frac{1}{T} \int_{t_0}^{t_0 + T}f(t) e^{-j n \omega_1 t} \, \mathrm{d}t, \qquad (n = 0, \pm 1, \pm 2, \cdots) \\ \end{aligned}\]
通常积分区间 \(t_0 \sim t_0 + T\) 取为 \(0 \sim T\) 或者 \(-\frac{T}{2} \sim \frac{T}{2}\) 。

显然,如果给出 \(F_n\) ,则 \(f(t)\) 唯一确定。在傅里叶级数的指数形式中傅里叶系数 \(F_n\) 不再是实数,而是复数,有模有辐角。

注: \(F_n\) 的推导过程如下: \(F_n = \frac{1}{2} (a_n - j b_n) = \frac{1}{T} \int_{t_0}^{t_0 + T} f(t) \cos (n \omega_1 t) \, \mathrm{d}t -j \frac{1}{T} \int_{t_0}^{t_0 + T} f(t) \sin (n \omega_1 t) \, \mathrm{d}t = \frac{1}{T} \int_{t_0}^{t_0 + T}f(t) e^{-j n \omega_1 t} \, \mathrm{d}t, \quad (n = 0, \pm 1, \pm 2, \cdots)\)

傅里叶级数的指数形式中系数 \(F_n\) 和傅里叶级数的三角形式中的幅度 \(c_n\) 和相位 \(\varphi_n\) 的关系如下:
\[\begin{aligned} & F_0 = c_0 \\ & F_n = |F_n| e^{j \varphi_n} , \quad (n = 1, 2, \cdots) \\ & F_{-n} = |F_{-n}| e^{-j \varphi_n} , \quad (n = 1, 2, \cdots) \\ & |F_n| = |F_{-n}| = \frac{1}{2} c_n , \quad (n = 1, 2, \cdots) \\ \end{aligned}\]

2.3.1 复数频谱(复数幅度谱,复数相位谱)

我们可以画出傅里叶级数指数形式表示的信号频谱。因为 \(F_n\) 一般是复函数,所有称这种频谱为复数频谱。根据 \(F_n = |F_n| e^{j \varphi_n}\) ,可以画出“复数幅度谱”为“复数相位谱”,如图 4 所示。

fourier_complex_spectrum.jpg

Figure 4: 周期函数的复数频谱示例

由于, \(|F_n| = |F_{-n}| = \frac{1}{2} c_n\) ,所以复数幅度谱是相对纵轴左右对称的。在三角函数的频谱(如图 3 )中每条谱线代表一个分量的幅度;在指数形式的频谱(如图 4 )中每个分量的幅度一分为二,在正、负频率相对应的位置上各为一半,所以,只有把正、负频率上对应的这两条谱线加起来才代表一个分量的幅度。

需要说明的是, 在复数频谱中出现的负频率是由于将 \(\sin(n \omega_1 t)\) 和 \(\cos(n \omega_1 t)\) 写成指数形式时,从数字的观点自然分成 \(e^{jn \omega_1 t}\) 以及 \(e^{-jn \omega_1 t}\) 两项,因而引入了 \(-jn \omega_1 t\) 项,所以,负频率的出现完全是数学运算的结果,并没有任何物理意义。

在 \(F_n\) 为实数的特殊情况下,显然相位不是0就是 \(\pi\) ,可以用 \(F_n\) 的正负分别表示相位 \(\varphi_n\) 为0和 \(\pi\) ,这时经常把幅度谱和相位谱画在同一张图上(而不是两张图)。后文将介绍这样的例子。

2.4 实例:求周期矩形脉冲信号的傅里叶频谱

设周期矩形脉冲信号(又称为方波信号) \(f(t)\) 的脉冲宽度为 \(\tau\) ,脉冲幅度为 \(E\) ,重复周期为 \(T\) (显然角频率 \(\omega_1 = 2 \pi/T\) ),如图 5 所示。求其傅里叶频谱。

fourier_ones.jpg

Figure 5: 周期矩形脉冲信号

直流分量 \(a_0\) ,余弦分量的幅度 \(a_n\) ,正弦分量的幅度 \(b_n\) 分别计算如下。
\[\begin{aligned} a_0 &= \frac{1}{T} \int_{-T/2}^{T/2} f(t) \, \mathrm{d}t = \frac{1}{T} \int_{-T/2}^{T/2} E \, \mathrm{d}t = \frac{E \tau}{T} \\ a_n &= \frac{2}{T} \int_{-T/2}^{T/2} f(t) \cos(n \omega_1 x) \, \mathrm{d}t = \frac{2E}{n \pi} \sin \left( \frac{n \pi \tau}{T} \right) , \qquad (n = 1,2,\cdots) \\ b_n &= \frac{2}{T} \int_{-T/2}^{T/2} f(t) \sin(n \omega_1 x) \, \mathrm{d}t = 0, \qquad (n = 1,2,\cdots) \\ \end{aligned}\]
从而可求得 \(c_0, c_n, \varphi_n \, (n=1,2,\cdots)\) 以及 \(F_n , \, (n = 0, \pm 1, \pm 2, \cdots)\) 。

其指数形式傅里叶级数幅度和相位如图 6 所示。

fourier_ones_spectrum.jpg

Figure 6: 方波信号指数形式傅里叶级数幅度谱和相位谱

由于, \(F_n\) 在这个例子中恰为实数(一般情况 \(F_n\) 为复数),此时相位不是0就是 \(\pi\) ,可以用 \(F_n\) 的正负分别表示相位 \(\varphi_n\) 为0和 \(\pi\) ,这时我们可以把幅度谱和相位谱画在同一张图上,如图 7 所示。

fourier_ones_spectrum_merge.jpg

Figure 7: 方波信号指数形式傅里叶级数幅度谱和相位谱(合画图)

3 傅里叶变换

前面介绍了连续“周期信号”的傅里叶级数,并得到了它的离散频谱。对于连续“非周期信号”,没有傅里叶级数,但可以推导出“傅里叶变换”。

由傅里叶级数的定义可知,对于周期为 \(T\) 的周期信号,离散频谱各个分量出现在 \(\omega_1, 2 \omega_1, 3 \omega_1, \cdots, (\omega_1 = 2 \pi /T)\) 处。显然, 当周期 \(T\) 越大,离散频谱的各个分量会越靠近(因为 \(\omega_1\) 会越小);当周期 \(T \to \infty\) 时(这时就是非周期信号),频谱会变为连续的。 以方波信号为例,其变化趋势如图 8 所示。

fourier_various_T.jpg

Figure 8: 当周期信号的周期 \(T\) 越大时,离散频谱的各个分量会越靠近

3.1 傅里叶变换定义

由前面的分析知,通过极限的方法,我们可以从周期信号的傅里叶级数导出非周期信号频谱的表示式,这称为傅里叶变换(Fourier transform)。在这里,不严格推导,直接给出非周期信号傅里叶变换的公式。

非周期信号 \(f(t)\) 的傅里叶变换为:

\begin{equation} F(\omega) = \mathscr{F}[f(t)] = \int_{- \infty}^{\infty} f(t) e^{-j \omega t} \, \mathrm{d}t \tag{傅里叶变换} \end{equation}

对应的傅里叶逆变换为:

\begin{equation} f(t) = \mathscr{F}^{-1}[F(\omega)] = \frac{1}{2 \pi} \int_{- \infty}^{\infty} F(\omega) e^{-j \omega t} \, \mathrm{d}t \tag{傅里叶逆变换} \end{equation}

式中, \(F(\omega)\) 是 \(f(t)\) 的频谱函数,它一般是复函数,可以写作:
\[F(\omega) = |F(\omega)| e^{j \varphi(\omega)}\]
其中 \(|F(\omega)|\) 是 \(F(\omega)\) 的模,它代表信号中各频率分量的相对大小, \(\varphi(\omega)\) 是 \(F(\omega)\) 的相位函数,它表示信号中各频率分量之间的相位关系。类似地, \(|F(\omega)| \sim \omega\) 与 \(\varphi(\omega) \sim \omega\) 曲线分别称为非周期信号的幅度频谱与相位频谱。

3.2 实例:求单边指数信号的傅里叶变换

已知,单边指数信号的表示式为:
\[f(t) = \begin{cases} e^{-at} & (t \ge 0) \\ 0 & (t<0) \end{cases}\]
其中 \(a\) 为正实数。

单边指数信号的傅里叶变换推导如下:
\[\begin{aligned} F(\omega) &= \int_{- \infty}^{\infty} f(t) e^{-j \omega t} \, \mathrm{d}t \\ &= \int_{0}^{\infty} e^{-at} e^{-j \omega t} \, \mathrm{d}t \\ &= \int_{0}^{\infty} e^{(a + j \omega) t} \, \mathrm{d}t \\ &= \frac{1}{a + j \omega} \\ \end{aligned}\]
从而:
\[\begin{aligned} & |F(\omega)| = \frac{1}{\sqrt{a^2 + \omega^2}} \\ & \varphi(\omega) = - \arctan \left(\frac{\omega}{a} \right) \\ \end{aligned}\]
单边指数信号的波形及频谱如图 9 所示。

fourier_transform_example.jpg

Figure 9: 单边指数信号的波形及频谱

3.3 卷积定理(Convolution theorem)

对于任意两个信号 \(f_1(t)\) 和 \(f_2(t)\) ,两者做卷积运算定义为:
\[f_1(t) \ast f_2(t) = \int_{- \infty}^{\infty} f_1(\tau) f_2(t - \tau)\, \mathrm{d} \tau\]

时域卷积定理:
设函数 \(f_1(t)\) 和 \(f_2(t)\) 的傅里叶变换为 \(F_1(\omega)\) 和 \(F_2(\omega)\) ,则:

\begin{equation} \mathscr{F}[f_1(t) \ast f_2(t)] = F_1(\omega) \cdot F_2(\omega) \tag{时域卷积定理} \end{equation}

证明如下:
\[\begin{aligned} \mathscr{F}[f_1(t) \ast f_2(t)] & = \int_{- \infty}^{\infty} \left[ \int_{- \infty}^{\infty} f_1(\tau) f_2(t - \tau)\, \mathrm{d} \tau \right] e^{-j \omega t} \, \mathrm{d}t \\ &= \int_{- \infty}^{\infty} f_1(\tau) \left[ \int_{- \infty}^{\infty} f_2(t - \tau) e^{-j \omega t} \, \mathrm{d} t \right] \, \mathrm{d} \tau \\ &= \int_{- \infty}^{\infty} f_1(\tau) F_2(\omega) e^{-j \omega t} \, \mathrm{d} \tau \\ &= F_2(\omega) \int_{- \infty}^{\infty} f_1(\tau) e^{-j \omega t} \, \mathrm{d} \tau \\ &= F_2(\omega) \cdot F_1(\omega) \\ \end{aligned}\]
它说明: 时域中两信号的卷积等效于在频域中频谱相乘。

频域卷积定理:

\begin{equation} \mathscr{F}[f_1(t) \cdot f_2(t)] = \frac{1}{2 \pi} F_1(\omega) \ast F_2(\omega) \tag{频域卷积定理} \end{equation}

频域卷积定理的证明过程和时域卷积定理的证明过程类似,这里不再重复。

4 离散傅里叶变换

设有限长序列 \(x(n)\) 长度为 \(N\) (在 \(0 \le n \le N-1\) 范围内),它的离散傅里叶变换(Discrete Fourier Transform, DFT) \(X(k)\) 仍然是一个长度为 \(N\) (在 \(0 \le n \le N-1\) 范围内)的频域有限长序列。
离散傅里叶变换为:

\begin{equation} X(k) = DFT[x(n)] = \sum_{n=0}^{N-1} x(n) W^{nk} \qquad (0 \le k \le N-1) \tag{离散傅里叶变换, DFT} \end{equation}

离散傅里叶逆变换为:

\begin{equation} x(n) = IDFT[X(k)] = \frac{1}{N} \sum_{k=0}^{N-1} X(k) W^{-nk} \qquad (0 \le n \le N-1) \tag{离散傅里叶逆变换, IDFT} \end{equation}

其中, \(W_N = e^{-j(2 \pi/N)}\) ,如果所讨论的问题中不涉及 \(N\) 的变动,可省略下标,简写为 \(W = e^{-j (2 \pi/ N)}\) 。

4.1 矩阵形式的离散傅里叶变换

前面介绍的离散傅里叶变换及逆变换的公式可以写为矩阵形式。
离散傅里叶变换的矩阵形式:
\[\begin{bmatrix} X(0) \\ X(1) \\ \vdots \\ X(N-1) \end{bmatrix} = \begin{bmatrix} W^0 & W^0 & W^0 & \cdots & W^0 \\ W^0 & W^{1 \times 1} & W^{2 \times 1} & \cdots & W^{(N-1) \times 1} \\ \vdots & \vdots & \vdots & & \vdots \\ W^0 & W^{1 \times (N-1)} & W^{2 \times (N-1)} & \cdots & W^{(N-1) \times (N-1)} \end{bmatrix} \begin{bmatrix} x(0) \\ x(1) \\ \vdots \\ x(N-1) \\ \end{bmatrix}\]

离散傅里叶逆变换的矩阵形式:
\[\begin{bmatrix} x(0) \\ x(1) \\ \vdots \\ x(n-1) \end{bmatrix} = \frac{1}{N} \begin{bmatrix} W^0 & W^0 & W^0 & \cdots & W^0 \\ W^0 & W^{-1 \times 1} & W^{-1 \times 2} & \cdots & W^{-1 \times (N-1)} \\ \vdots & \vdots & \vdots & & \vdots \\ W^0 & W^{-(N-1) \times 1} & W^{- (N-1) \times 2} & \cdots & W^{-(N-1) \times -(N-1)} \end{bmatrix} \begin{bmatrix} X(0) \\ X(1) \\ \vdots \\ X(N-1) \\ \end{bmatrix}\]

上面两个公式简写为:

\begin{equation} \boldsymbol{X}(k) = \boldsymbol{W}^{nk} \boldsymbol{x}(n) \tag{离散傅里叶变换的矩阵形式} \end{equation}


\begin{equation} \boldsymbol{x}(n) = \frac{1}{N} \boldsymbol{W}^{-nk} \boldsymbol{X}(k) \tag{离散傅里叶逆变换的矩阵形式} \end{equation}

此处, \(\boldsymbol{X}(k)\) 与 \(\boldsymbol{x}(n)\) 分别为 \(N\) 列的列矩阵,而 \(\boldsymbol{W}^{nk}\) 与 \(\boldsymbol{W}^{-nk}\) 分别为 \(N \times N\) 方阵(这两个方阵都是对称矩阵)。

4.2 实例:求有限序列的离散傅里叶变换


例:求 \(x(n) = \{1,1,1,1 \}\) (如图 10 (a)所示)的离散傅里叶变换,再由离散傅里叶逆变换求得 \(x(n)\) ,验证结果的正确性。
由 \(N=4\) 得到 \(W^1 = e^{-j (2 \pi/4)} = -j\) ,直接应用离散傅里叶变换公式,有:
\[\begin{aligned} \begin{bmatrix} X(0) \\ X(1) \\ X(2)\\ X(3) \end{bmatrix} &= \begin{bmatrix} W^0 & W^0 & W^0 & W^0 \\ W^0 & W^1 & W^2 & W^3 \\ W^0 & W^2 & W^4 & W^6 \\ W^0 & W^3 & W^6 & W^9 \end{bmatrix} \begin{bmatrix} x(0) \\ x(1) \\ x(2) \\ x(3) \end{bmatrix} \\ &= \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & -j & -1 & j \\ 1 & -1 & 1 & -1 \\ 1 & j & -1 & -j \end{bmatrix} \begin{bmatrix} 1 \\ 1 \\ 1 \\ 1 \end{bmatrix} = \begin{bmatrix} 4 \\ 0 \\ 0 \\ 0 \end{bmatrix} \\ \end{aligned}\]

离散傅里叶逆变换为:
\[\begin{bmatrix} x(0) \\ x(1) \\ x(2)\\ x(3) \end{bmatrix} = \frac{1}{4} \begin{bmatrix} W^0 & W^0 & W^0 & W^0 \\ W^0 & W^{-1} & W^{-2} & W^{-3} \\ W^0 & W^{-2} & W^{-4} & W^{-6} \\ W^0 & W^{-3} & W^{-6} & W^{-9} \end{bmatrix} \begin{bmatrix} 4 \\ 0 \\ 0 \\ 0 \end{bmatrix} = \begin{bmatrix} 1 \\ 1 \\ 1 \\ 1 \end{bmatrix}\]
显然,通过离散傅里叶逆变换可求得原始信号。

fourier_DFT_example.jpg

Figure 10: 离散傅里叶变换实例

说明:在这个例子中,离散傅里叶变换得到的序列恰好都为实数。不过,离散傅里叶变换得到的序列一般为复数,比如 \(x(n) = \{1,2,3,4 \}\) 的离散傅里叶变换是一个复数序列,计算如下:
\[\begin{aligned} \begin{bmatrix} X(0) \\ X(1) \\ X(2)\\ X(3) \end{bmatrix} &= \begin{bmatrix} W^0 & W^0 & W^0 & W^0 \\ W^0 & W^1 & W^2 & W^3 \\ W^0 & W^2 & W^4 & W^6 \\ W^0 & W^3 & W^6 & W^9 \end{bmatrix} \begin{bmatrix} x(0) \\ x(1) \\ x(2) \\ x(3) \end{bmatrix} \\ &= \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & -j & -1 & j \\ 1 & -1 & 1 & -1 \\ 1 & j & -1 & -j \end{bmatrix} \begin{bmatrix} 1 \\ 2 \\ 3 \\ 4 \end{bmatrix} = \begin{bmatrix} 10 \\ -2+j2 \\ -2 \\ -2-j2 \end{bmatrix} \\ \end{aligned}\]

4.3 快速傅里叶变换(FFT)

快速傅里叶变换(Fast Fourier transform, FFT)是一种计算离散傅里叶变换的快速算法。

4.3.1 DFT的计算量

我们先来了解一下DFT的计算量。
由DFT定义知,将 \(x(n)\) 与 \(W^{nk}\) 两两相乘再取和即可得到 \(X(k)\) ,即每计算一个 \(X(k)\) 值需要进行 \(N\) 次复数相乘(一个复数乘法包括4个实数乘法和2个实数加法),和 \(N-1\) 次复数相加。对于 \(N\) 个 \(X(k)\) 点,应重复 \(N\) 次上述运算。因此, 要完成全部DFT运算共需 \(N^2\) 次复数乘法和 \(N(N-1)\) 次复数加法。
例如,当 \(N=4\) 时,DFT的矩阵表示式为:
\[\begin{bmatrix} X(0) \\ X(1) \\ X(2)\\ X(3) \end{bmatrix} = \begin{bmatrix} W^0 & W^0 & W^0 & W^0 \\ W^0 & W^1 & W^2 & W^3 \\ W^0 & W^2 & W^4 & W^6 \\ W^0 & W^3 & W^6 & W^9 \end{bmatrix} \begin{bmatrix} x(0) \\ x(1) \\ x(2) \\ x(3) \end{bmatrix}\]
显然,求每个 \(X(k)\) 值,需要 \(N=4\) 次复数乘法和 \(N-1=3\) 次复数加法,要完成全部DFT运算共需 \(N^2=16\) 次复数乘法和 \(N(N-1)\) 次复数加法。
随着 \(N\) 值增大,运算工作量将迅速增长,例如 \(N=10\) 时需要100次复数乘法,而当 \(N=1024\) 时,就需要1048576(一百多万次)次复数乘法运算。所以,当 \(N\) 较大时,DFT的计算量太大,无法对信号进行实时处理。因此,我们需要一种更快的算法来计算DFT。

4.3.2 Cooley–Tukey FFT algorithm

从FFT算法诞生至今,有很多改进和派生的离散傅里叶变换快速算法。最基本的FFT是库利-图基FFT算法(Cooley–Tukey FFT algorithm)。

库利-图基FFT算法描述(这里不证明它的正确性,直接给出结论)如下:
设序列 \(x(n)\) 的点数 \(N\) 为2的整数次幕(即 \(N=2^L\) , \(L\) 为整数,如果不满足这个条件,可以人为地加上若干0值点,使之达到这一要求)。将 \(x(n), (n=0,1,\cdots,N-1)\) 先按 \(n\) 的偶数子序列和奇数子序列分成以下两组:
\[\begin{aligned} x_1(r) &= x(2r), \qquad (r=0,1,\cdots, \frac{N}{2} - 1) \\ x_2(r) &= x(2r+1), \qquad (r=0,1,\cdots, \frac{N}{2} - 1) \end{aligned}\]
记 \(X_1(k)\) 和 \(X_2(k)\) 分别是 \(x_1(r)\) 和 \(x_2(r)\) 的 \(N/2\) 点DFT, \(W_N^k=e^{−j(2 \pi /N)k}\) ,则 \(x(n)\) 的离散傅里叶变换的前 \(N/2\) 个值和后 \(N/2\) 个值可分别由下面式子(FFT核心公式)求得:
\[\begin{aligned} & \color{red}{X(k) = X_1(k) + W_N^k X_2(k), \qquad (k=0,1,\cdots, \frac{N}{2}-1)} \\ & \color{red}{X(k+\frac{N}{2}) = X_1(k) - W_N^k X_2(k), \qquad (k=0,1,\cdots, \frac{N}{2}-1)} \end{aligned}\]
这样,只要求得0到 \(\frac{N}{2} -1\) 区间的 \(X_1(k)\) 和 \(X_2(k)\) 值,就可以求出0到 \(N-1\) 区间内的所有 \(X(k)\) 值,显然,这是一种分治策略。对于子序列 \(X_1(k)\) 和 \(X_2(k)\) 的求解可以使用同样的思路。图 11 表示了FFT算法中利用两个子序列的DFT求得整个序列的DFT的过程。

fourier_FFT.jpg

Figure 11: FFT算法,蝶形信号流图(支路上没有系数时表示该支路传输系数为1)

参考:
《信号与系统(郑君里,第二版,下册)》9.6节 快速傅里叶变换
《数字信号处理教程(第三版,程佩青)》第四章 快速傅里叶变换

4.3.3 实例:用FFT算法计算DFT

用FFT算法求解 \(x(n) = \{1,2,3,4 \}\) 的DFT。

这个例子中 \(N=4\) ,设其偶数子序列 \(x_1(n) = \{1, 3\}\) 的DFT为 \(X_1(k)\) ,奇数子序列 \(x_2(n) = \{2, 4\}\) 的DFT为 \(X_2(k)\) ,直接利用FFT算法公式有:
\[\begin{aligned} X(0) = X_1(0) + W_4^0 X_2(0) \\ X(1) = X_1(1) + W_4^1 X_2(1) \\ X(2) = X_1(0) - W_4^0 X_2(0) \\ X(3) = X_1(1) - W_4^1 X_2(1) \end{aligned}\]

\(X_1(k)\) 和 \(X_2(k)\) 如何求解呢?也可以用FFT算法,但序列已经很短了,直接用DFT原始定义公式即可:
\[\begin{aligned} & X_1(0) = x(0) + W_2^0 x(2) = 1 + 3 = 4 \\ & X_1(1) = x(0) - W_2^0 x(2) = 1 - 3 = -2 \\ & X_2(0) = x(1) + W_2^0 x(3) = 2 + 4 = 6 \\ & X_2(1) = x(1) - W_2^0 x(3) = 2 - 4 = -2 \end{aligned}\]

从而, \(N=4\) 时FFT算法如图 12 所示。

fourier_FFT_example.jpg

Figure 12: \(N=4\) 时FFT算法

容易算得 \(W_4^0 = 1, W_4^1 = e^{-j (2 \pi/4)} = -j\) ,从而:
\[\begin{aligned} & X(0) = X_1(0) + W_4^0 X_2(0) = 4 + 6 = 10 \\ & X(1) = X_1(1) + W_4^1 X_2(1) = -2 + (-j) (-2) = -2 + 2j \\ & X(2) = X_1(0) - W_4^0 X_2(0) = 4 - 6 = -2 \\ & X(3) = X_1(1) - W_4^1 X_2(1) = -2 - (-j) (-2) = -2 -2j \end{aligned}\]

这个结果和不使用FFT算法,直接使用DFT定义公式得到的结果一样(参见 4.2 )。

4.3.4 FFT的计算量

FFT把一个 \(N\) 点的DFT分解为两个 \(N/2\) 点DFT,如果直接计算 \(N/2\) 点DFT,则由前面分析知一个 \(N/2\) 点DFT需要 \((\frac{N}{2})^2 = \frac{N^2}{4}\) 次复数乘法, \(\frac{N}{2}(\frac{N}{2} -1)\) 次复数加法,两个 \(N/2\) 点DFT,则共需要 \(2 \times (\frac{N}{2})^2 = \frac{N^2}{2}\) 将复数乘法和 \(N(\frac{N}{2} - 1)\) 次复数加法。此外把两个 \(N/2\) 点DFT合成 \(N\) 点DFT时,还需要 \(\frac{N}{2}\) 次复数乘法和 \(N\) 次复数加法。所以总共需要 \(\frac{N^2}{2} + \frac{N}{2} = \frac{N(N+1)}{2} \approx \frac{N^2}{2}\) 次复数乘法和 \(N(\frac{N}{2} - 1) + N = \frac{N^2}{2}\) 次复数加法。因些,通过把一个 \(N\) 点的DFT分解为两个 \(N/2\) 点DFT后,运算的计算量差不多减少到一半。

显然,对于分解的两个 \(N/2\) 点DFT,我们不用直接计算,可以再次利用FFT。从而 \(N\) 点FFT的运算量会进一步减少。通过更细致的分析可以FFT的复数乘法次数是 \(\frac{N}{2} \log_2 N\) ,从而直接计算DFT与FFT算法的计算量(只考虑了复数乘法)之比为:
\[\frac{N^2}{\frac{N}{2} \log_2 N} = \frac{2N}{\log_2 N}\]
直接计算DFT与FFT算法的计算量(只考虑了复数乘法)与序列点数 \(N\) 的关系曲线如图 13 所示。

fourier_FFT_vs_DFT.jpg

Figure 13: 直接计算DFT与FFT算法的计算量(只考虑了复数乘法)的比较

4.3.5 FFT逆变换(快速计算离散傅里叶逆变换)

如何快速地求解离散傅里叶逆变换(IDFT)呢?

我们观察一下离散傅里叶逆变换(IDFT)的公式:
\[x(n) = IDFT[X(k)] = \frac{1}{N} \sum_{k=0}^{N-1} X(k) W^{-nk} \qquad (0 \le n \le N-1)\]
与离散傅里叶变换(DFT)公式:
\[X(k) = DFT[x(n)] = \sum_{n=0}^{N-1} x(n) W^{nk} \qquad (0 \le k \le N-1)\]
比较两个公式可知,只要把DFT运算中的每一个系数 \(W^{nk}\) 换成 \(W^{-nk}\) ,最后再乘以常数 \(1/N\) ,则FFT算法就可以用来运算IDFT了。但这需要稍微调整FFT程序代码。

下面讨论一种完全不用改变FFT的程序就可以快速计算IDFT的方法。对IDFT公式取共轭:
\[x^{*}(n) = \frac{1}{N} \sum_{k=0}^{N-1} X^{*}(k) W^{nk}\]
从而:
\[x(n) = \frac{1}{N} \left[ \sum_{k=0}^{N-1} X^{*}(k) W^{nk} \right]^{*} = \frac{1}{N} \left[ DFT[X^{*}(k)] \right]^{*}\]

这说明, 只要先将 \(X(k)\) 取共轭,就可以直接利用FFT程序,最后再将运算结果取一次共轭,并乘以 \(1/N\) ,即得 \(x(n)\) 值。 因此,FFT运算和IFFT(Inverse Fast Fourier Transform)运算就可以共用一个子程序模块了,这样就很方便了。


Author: cig01

Created: <2016-07-02 Sat 00:00>

Last updated: <2018-01-02 Tue 16:56>

Creator: Emacs 25.3.1 (Org mode 9.1.4)