XMU《计算方法》实验二 FFT

实验二 FFT

一、MATLAB 代码

clear;N = 32; TIME = 5;X = linspace(-pi, pi, 33);
X = X(1 : 32);
A = X .^ 2 .* cos(X);for m = 0 : N-1w(m+1) = exp(1i * 2 * pi / 32 * m);
endi = 1;
while i < NB = A;for j = 0 : i*2 : N-1for k = 0 : i-1A(j + k + 1) = B(j / 2 + k + 1) + B(j / 2 + N / 2 + k + 1);A(j + k + i + 1) = (B(j / 2 + k + 1) - B(j / 2 + N / 2 + k + 1)) * w(1 + j / 2);endendi = i * 2;% A
endC = A;a = [];
b = [];
for k = 1 : 17a = [a, 1/16 * real((-1)^k * C(k))];b = [b, 1/16 * imag((-1)^k * C(k))];
enda, bx = [-pi : 0.01 : pi];
for k = 1 : length(x)y(k) = a(1)/2;
end
for k = 2 : 17y = y + cos((k-1) * x) * a(k) + sin((k-1) * x) *  b(k);
end
plot(x, y)hold on
yy = x .^ 2 .* cos(x);
plot(x, yy, 'r')legend('拟合图像', '原始函数')

二、运行结果

答案的三角插值多项式为 \(f(x) = \frac 12a_0 + \sum_{k=1}^{16} (a_i\cos(kx) + b_i\sin(kx))\)，其中 \(a_0\) 对应程序中的 \(a(0)\)，\(a_i = a(i+1), b_i = b(i+1)\)。

三、傅里叶变换

傅里叶变换是一种将函数从时域转换到频域的数学工具。它通过计算不同频率成分的振幅和相位来表示一个函数，从而将一个时域信号转化为频域信号。很多在时域上不可分解的东西，转成频域后很容易过滤和识别。

在离散傅里叶变换的过程中，我们将时域信号的离散序列求得的频域序列。

比如卷积运算，在数字信号处理、模式识别等领域，很多时候我们都需要对于两个时域信号求卷积，而且通常我们需要求卷积的都是离散时域序列，朴素的卷积算法的时间复杂度是 \(\Theta(n^2)\)。通过时域卷积定理，我们有

\[F[f_1(x) * f_2(x)] = F[f_1(x)]\cdot F[f_2(x)] \]

进而可以将时域的卷积转化为频域的相乘，而通过快速傅里叶变换，可以在 \(\Theta(n\log n)\) 的时间复杂度中完成这个任务。

四、离散傅里叶变换（DFT）

假设我们有 \(f(x)\) 是一个以 \(2\pi\) 为周期的复函数，\(N\) 个等分点上的函数值：

\[y_i = f(x_i) \]

那么它的傅里叶逼近可以表示为

\[S(x) =\frac 12a_0 + \sum_{k=1}^{16} (a_i\cos(kx) + b_i\sin(kx))\\= \sum_{k=0}^{N-1} c_ke^{ikx} \]

其中

\[c_j = \frac 1N \sum_{k=0}^{N-1}y_ke^{-ijk\frac{2\pi}N}\\ y_j = \sum_{k=0}^{N-1}c_k e^{ik\frac{2\pi}Nj} \]

五、快速傅里叶变换（FFT）

我们假设 \(\omega_N = e^{i\frac{2\pi}N} = \cos\frac{2\pi}N + i\sin\frac{2\pi}N\)。

DFT 本质上就是要计算 \(c_j = \sum_{k=0}^{N-1} x_kw_N^{jk}\)。

对于 \(\omega_N\) 数组，我们可以将之看作将一个单位圆分为 \(N\) 等份，因此很显然 \(\omega_N\) 满足如下的一些性质：

\[\omega_N^j\omega_N^k = \omega_N^{j+k}\\ \omega_N^{jN+k} = \omega_N^k\\ \omega_N^{jk+N/2} = -\omega_N^{jk}\\ \omega_{jN}^{jk} = \omega_N^k \]

如果 \(N\) 是 \(2\) 的倍数，那么我们可以将计算 \(c_j\) 的式子继续划分：

\[\begin{align*} c_j &= \sum_{k=0}^{N-1}x_k\omega_N^{jk}\\ &= \sum_{k=0}^{\frac N2 - 1} x_k\omega_N^{jk} + \sum_{k=0}^{\frac N2-1} x_{\frac N2 + k}\omega_N^{j(\frac N2 + k)}\\ &= \sum_{k=0}^{\frac N2 - 1} x_k\omega_N^{jk} + (-1)^j\sum_{k=0}^{\frac N2-1} x_{\frac N2 + k}\omega_N^{jk}\\ &= \sum_{k=0}^{\frac N2 - 1} (x_k +(-1)^j x_{\frac N2 + k})\omega_N^{jk} \end{align*} \]

对于 \(j\) 讨论其奇偶性可得，

\[c_{2j} = \sum_{k=0}^{\frac N2 - 1} (x_k + x_{\frac N2 + k})\omega_{\frac N2}^{jk}\\ c_{2j+1} = \sum_{k=0}^{\frac N2 - 1} (x_k - x_{\frac N2 + k})\omega_N^k\omega_{\frac N2}^{jk} \]

令 \(y_k = x_k + x_{\frac N2 + k}, y_{\frac N2 + k} = x_k - x_{\frac N2 + k}\)。

则上式可化为

\[c_{2j} = \sum_{k=0}^{\frac N2 - 1} y_k\omega_{\frac N2}^{jk}\\ c_{2j+1} = \sum_{k=0}^{\frac N2 - 1} y_{\frac N2 +k}\omega_N^k\omega_{\frac N2}^{jk} \]

我们只需要反复使用这个公式，将计算的过程分治下去，当达到 \(N=1\) 时，我们就不需要继续递归了。

显然，如果我们的递归过程模拟正确，那么我们很容易就可以求出正确的离散傅里叶变换的结果。

六、\(f(x) = x^2\cos x\) 的三角插值多项式

我们如果需要计算 \(f(x) = x^2\cos x\) 的三角插值多项式，就是要计算出这个函数相关的 \(c_j\) 的值，然后通过如下方法，计算出其三角多项式。

\[\begin{align*} &\frac 1{16} c_j(-1)^j \\ =& \frac 1{16} c_j e^{-i\pi j}\\ =& \frac 1{16}\sum_{k=0}^{31} f_j[\cos j(-\pi + \frac 1{16}k) + i\sin j(-\pi + \frac 1{16}k)]\\ =& \frac 1{16}\sum_{k=0}^{31} f_j(\cos jy_k + i\sin jy_k) \end{align*} \]

由 \(a_j + ib_j = \frac 1{16} c_j(-1)^j\)，得

\[a_j = Re[\frac 1{16} c_j(-1)^j]\\ b_j = Im[\frac 1{16} c_j(-1)^j] \]

即可得到 \(a_j\) 和 \(b_j\) 的两个系数序列。

因为 \(16\) 次的三角插值多项式实际上需要 \(17\) 阶的序列，因此我们取了 \(32\) 个插值节点，只取结果的最高的 \(17\) 项。

七、理解和感悟

傅立叶变换是一种数学工具，用于将信号从时域转换到频域。它被广泛应用于信号处理、图像处理、音频处理等领域。通过傅立叶变换，我们可以将信号分解为不同频率的正弦波，从而更好地分析和理解信号的性质。

快速傅立叶变换是一种优化过的傅立叶变换算法，能够在较短的时间内计算出大量的离散傅立叶变换结果。相对于朴素的傅立叶变换算法，快速傅立叶变换具有更高的计算效率和更小的计算复杂度。它被广泛应用于信号处理、图像处理、音频处理、数据压缩等领域。

这是非常重要和强大的数学工具，傅立叶变换和快速傅立叶变换可以帮助我们更好地理解和分析信号、图像、音频等数据。同时，它们也为实际问题提供了有效的解决方案。这使得我们可以更好地探索和利用数据，推动科学和工程领域的发展。

通过学习傅立叶变换和快速傅立叶变换，我们还可以更好地了解数学在实际问题中的应用和价值。这种数学思维和方法的应用可以帮助我们更好地理解和处理现实世界中的问题，从而使我们的生活更加便捷和舒适。

傅立叶变换和快速傅立叶变换是非常重要和强大的数学工具，它们不仅可以帮助我们更好地理解和分析数据，还可以推动科学和工程领域的发展。