多项式与快速傅里叶变换

多项式

什么是多项式？

形如：

对于确定的 $P(x)$ 我们需要知道 $n+1$ 个系数，或者需要 $n+1$ 个点 $\left(x_i,P(x_i)\right)$

如果存在两个多项式 $P(x)$ 和 $Q(x)$ 都过 $n+1$ 个点，则

因为 $Q(x_i)=P(x_i)$ 并且范德蒙行列式满秩， 所以 $p_i = q_i$

什么是多项式乘法？

我所要讲的是现在已知两个多项式 $A(x)$ 与 $B(x)$ ，我们应该如何得到两个多项式的乘积$C(x)$？

上述的公式中 $a(x),b(x)$ 已知，$c_r$ 未知 通常的算法：多项式乘法（假设 $p \geqslant q$）

多项式乘法代码实现

# 定义多项式系数
a = [1, 2, 3, 4, 5]  # 表示多项式 x^4 + 2x^3 + 3x^2 + 4x + 5
b = [2, 3, 4, 5]     # 表示多项式 2x^3 + 3x^2 + 4x + 5

# 初始化结果多项式的系数数组
# 结果多项式的最大次数为 len(a) + len(b) - 2
result_length = len(a) + len(b) - 1
result = [0] * result_length

# 反转系数数组以便从低次项到高次项进行计算
a_reversed = a[::-1]
b_reversed = b[::-1]

# 多项式乘法计算
for j in range(len(b_reversed)):
    for i in range(len(a_reversed)):
        result[i + j] += a_reversed[i] * b_reversed[j]

# 反转结果数组以恢复正确的次序
result = result[::-1]

# 打印结果
print(result)

---

时间复杂度

我们讨论一下计算的复杂度 所需要的计算：

• 浮点数加法 $p \cdot q$
• 浮点数乘法 $p \cdot q$

现在我想告诉你的是有一种算法可以将算法的复杂度减到 $O(N \cdot \log_2N)$

我们即将用到的大杀器 快速傅里叶变换

傅里叶变换

泰勒级数

其中$\{x^n| n\in\{0,1,2,3\cdots\}$是一组基（虽然不是正交基）$a_i$是这一组基坐标系下的坐标 。

可以使用一些判别法求出泰勒级数的收敛半径

傅里叶级数

在讨论傅里叶级数之前我们先证明一些之后可能用到的小推论

核函数 $D_N(x)$

核函数 $D_N(x)$ 的积分

三角函数的正交

$\cos(mx)$ 和 $\cos(nx)$ 当 $m \ne n$ 时正交
$\sin(mx)$ 和 $\sin(nx)$ 当 $m \ne n$ 时正交
$\cos(mx)$ 和 $\sin(nx)$ 当 $m \ne n$ 时正交

验证三角函数族$\{ \left. 1,sin(nx),cos(nx) \right | n\in\{1,2,3\cdots\}$在区间$[-\pi,\pi]$上的正交性

我们使用以下内积定义：

证明：

---

计算傅里叶级数的系数

若函数$f(x)$可以展开成为$a_0, a_{n}, b_{n}$成立

计算$a_0,a_n,b_n$分别是多少

贝塞尔不等式的证明

贝塞尔不等式: 其中$f(x)$为任意函数

要证明贝塞尔不等式，我们可以考虑误差项。 定义 $S_N(x)$ 为傅里叶级数的第 $N$ 部分和：

误差项 $R_N(x)$ 为：

误差项 $R_N^2(x)$ 为：

因为贝塞尔不等式成立 如果 $f(x)$ 可积 $\to \lim_{n \to \infty}a_n,b_n = 0,0$

黎曼引理的另一种证明方式

傅里叶级数收敛性的证明

将 $a_0, a_n, b_n$ 带入：

所以：

将积分中的变量用别的字母替代 $x \to t$：

将 $\cos nx, \sin nx$ 放入积分内：

交换积分与求和顺序：

合并 $f(t)$ 和 $\cos nt \cos nx$, $\sin nt \sin nx$：

将 $\frac{1}{2\pi}\int_{-\pi}^{\pi}f(t)\,\text dt$ 拿进去，$f(x)$ 拿出来：

换元，令 $t-x=u \to t = u+x$：

因为 $\{\cos nu | n \ge 1\}$ 都是周期函数，所以 $\int_{-\pi-x}^{\pi-x} \to \int_{-\pi}^{\pi}$：

由$D_N(x)$：

现在让 $N \to +\infty$, $x = x_0$：

将 $\int_{-\pi}^{\pi} \to \int_{-\pi}^{0} + \int_{0}^{\pi}$:

对 $[-\pi,0]$ 令 $-u = v$：

合并 $u, v \to t$：

同时减去 $s = \frac{f_+(x_0) + f_-(x_0)}{2}$，其中 $s$ 是 $f(x_0)$ 的左右极限除以 2：

将 $\int_{0}^{\pi} = \int_{0}^{\delta} + \int_{\delta}^{\pi}$:

由于 $\int_{\delta}^{\pi} = 0$:

由于 $f(x)$ 可积
或者说 $\int_{0}^{\pi} \frac{f(x_0+t) - f_+(x_0) + f(x_0-t) - f_-(x_0)}{t} \,\text dt$ 收敛
或者 $\lim_{\delta \to 0^+} \int_{0}^{\delta} \frac{f(x_0+t) - f_+(x_0) + f(x_0-t) - f_-(x_0)}{t} \,\text dt = 0$:

最终得出结论：

欧拉公式

根据欧拉公式 $e^{i\theta}=\cos\theta + i\sin\theta$，我们可以对其进行变形来得到 $\cos\theta$ 和 $\sin\theta$ 的表达式。
已知 $e^{i\theta}=\cos\theta + i\sin\theta$，同时 $e^{-i\theta}=\cos(-\theta)+i\sin(-\theta)=\cos\theta - i\sin\theta$

傅里叶级数的复数表示

将 $\cos n\theta, \sin n\theta$ 带入傅里叶级数：

合并 $e^{in\theta}, e^{-in\theta}$ 的系数 ：

合并 $\sum$ ：

$c_n$ 的复数表示:

傅里叶级数的复数表示:

周期为T的傅里叶级数的复数表示

求 $a_n, b_n$ :

$c_n$ 的复数表示:

e 形式的正交基

要判断一组复指数函数是否正交，我们可以考虑以下形式的函数：

其中 $n$ 和 $m$ 是不同的整数。为了判断 $e^{i\frac{2\pi}{T} nt}$ 和 $e^{i\frac{2\pi}{T} mt}$ 的正交性，我们需要计算它们在一个周期 $T$ 内的内积：

计算这个积分：

当 $n = m$ 时：

当 $n \neq m$ 时：

连续傅里叶变换

从离散到连续，我们考虑 $f(t)$ 是一个非周期信号，即周期 $T \to \infty$。在这种情况下，傅里叶级数变成傅里叶变换。具体做法是将离散的频率分量变为连续的频率分量。

频率分量的连续化

为了将傅里叶级数的频率分量从离散变为连续，我们用 $\omega = n \omega_0$ 替换 $n$。注意到 $\omega_0 = \frac{2\pi}{T}$，当 $T \to \infty$ 时，$\omega_0 \to 0$，频率分量变得连续。

所以，我们有： $\omega = n \omega_0 = n \frac{2\pi}{T}$

这时，傅里叶系数 $c_n$ 可以表示为： $c_n = \frac{1}{T} \int_{0}^{T} f(t) e^{-i n \omega_0 t} \, dt$ 由于 $\omega = n \omega_0$，我们可以改写为：

极限过程

当 $T \to \infty$ 时，$\omega_0 \to 0$，频率分量变得连续。我们可以将傅里叶系数 $c_n$ 变为傅里叶变换中的连续频谱 $F(\omega)$。

考虑到 $\omega$ 是连续的，我们需要引入一个密度函数来替代 $c_n$。傅里叶系数 $c_n$ 的表达式变为：

当 $T \to \infty$ 时，这个积分的上下限变为 $-\infty$ 到 $\infty$，所以我们得到傅里叶变换的公式：

连续傅里叶逆变换

傅里叶级数的逆变换公式为：

当 $T \to \infty$ 时，离散的和 (sum) 变为积分 (integral)，并且 $c_n$ 变为 $F(\omega)$：

通过将傅里叶级数的周期 $T$ 取到无限大，我们可以将离散的频率分量变为连续的频率分量，从而推导出傅里叶变换。这展示了傅里叶级数和傅里叶变换在处理周期性和非周期性信号时的内在联系。

离散傅里叶变换

对于周期信号 $x[n]$，我们可以使用傅里叶级数来表示它：

其中 $X[k]$ 是傅里叶系数。

复指数函数的正交性

离散傅里叶变换系数

为了找到傅里叶变换系数 $X[k]$，我们将上述傅里叶级数公式两边乘以 $e^{-i \frac{2\pi}{N} kn}$ 并对 $n$ 求和：

交换求和次序：

这样 $k=0,1,2,3 \cdots$ 代入上式

离散傅里叶逆变换

DFT:

$X[k] = \sum_{n=0}^{N-1} x[n] \cdot e^{-i \frac{2\pi}{N} kn}$

IDFT:

$x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] \cdot e^{i \frac{2\pi}{N} kn}$

离散傅里叶变换可以用矩阵形式来表示，这种表示方式在计算机实现中非常有用，因为它将DFT视为矩阵与向量的乘积。

离散傅里叶变换矩阵

假设我们有一个长度为 $N$ 的离散序列 $\mathbf{x} = [x[0], x[1], \ldots, x[N-1]]^T$。DFT 可以表示为矩阵乘法：

其中：

• $\mathbf{X} = [X[0], X[1], \ldots, X[N-1]]^T$ 是频域表示的结果向量。
• $\mathbf{W}_N$ 是 $N \times N$ 的DFT矩阵。

DFT 矩阵 $\mathbf{W}_N$ 的元素由以下公式给出：

其中 $k, n = 0, 1, 2, \ldots, N-1$。具体来说，矩阵 $\mathbf{W}_N$ 的第 $k$ 行，第 $n$ 列的元素是：

由此可知当n为偶数时， $W_{N}[k, n]=W_N[k+\frac{N}{2},n]$

矩阵形式的DFT计算

将信号 $\mathbf{x}$ 与 DFT 矩阵 $\mathbf{W}_N$ 相乘，即可得到频域表示 $\mathbf{X}$：

其中 $W = e^{-i \frac{2\pi}{N}}$ 是单位根。

逆离散傅里叶变换矩阵

逆离散傅里叶变换也可以用矩阵形式表示，其矩阵 $\mathbf{W}_N^*$ 是DFT矩阵的共轭转置，并且需要乘以 $\frac{1}{N}$：

其中 $\mathbf{W}_N^*$ 的元素为：

快速傅里叶变换

快速傅里叶变换是一种 求解离散傅里叶变换的 算法

离散傅里叶变换

利用 $e^{−i\frac{2π}{N}kn}$ 的周期性,通过观察可以发现，当$n$为偶数时（我们用$2n$表示）

由此可知当n为偶数时

当$n$为奇数时，我们可以提取一个 $e^{−i\frac{2π}{N}k}$ 公因数，让剩下的 $W_{N}[k, n]$ 满足$n$为偶数的条件。

将 $X[k]$ 分为 $n$ 为偶数以及 $n$ 为奇数两个部分。

我们令偶数部分为 $E[k]$ ，奇数部分提取公因数后的部分为 $O[k]$ ，则 

对于上式 $k \lt \frac{N}{2}$ 成立，当 $\frac{N}{2} \leqslant k \lt N$ 反转一圈

$W_N[k,1]$反转半圈，相当于取反

快速傅里叶变换公式 $k \lt \frac{N}{2}$

代码实现

import cmath

def fft(x):
    N = len(x)
    if N <= 1:
        return x
    elif N % 2 != 0:
        raise ValueError("Size of x must be a power of 2")

    # FFT of even and odd indexed elements
    even = fft(x[0::2])
    odd = fft(x[1::2])

    # Combine
    T = [cmath.exp(-2j * cmath.pi * k / N) * odd[k] for k in range(N // 2)]
    return [even[k] + T[k] for k in range(N // 2)] + \
           [even[k] - T[k] for k in range(N // 2)]

# 示例数据
x = [0, 1, 0, -1]
result = fft(x)
for i, val in enumerate(result):
    print(f"X[{i}] = {val}")

时间复杂度

分解到最后只有一项，只需要进行一次乘法运算即可。所以$T(1)=1$。

多项式与快速傅里叶变换

用时间复杂度 $O(N \cdot \log_2N)$ 求解 $c_r$ 需要三步骤

应用快速傅里叶变换转换到点值表示

对于多项式 $A(x) = a_nx^n+a_{n-1}x^{n-1}+\cdots+a_0$ 选出大于 $p+q$ 的最小的 2 的幂次 $N$ 。 对 $A(x),B(x)$ 分别做快速傅里叶变换

逐点相乘

$C[n] = A[n] \cdot B[n]$

逆快速傅里叶转换回系数表示

代码实现

import cmath


def fft(x):
    N = len(x)
    if N <= 1:
        return x
    elif N % 2 != 0:
        raise ValueError("Size of x must be a power of 2")

    # FFT of even and odd indexed elements
    even = fft(x[0::2])
    odd = fft(x[1::2])

    # Combine
    T = [cmath.exp(-2j * cmath.pi * k / N) * odd[k] for k in range(N // 2)]
    return [even[k] + T[k] for k in range(N // 2)] + [
        even[k] - T[k] for k in range(N // 2)
    ]


def ifft(x):
    N = len(x)
    if N <= 1:
        return x
    elif N % 2 != 0:
        raise ValueError("Size of x must be a power of 2")

    # FFT of even and odd indexed elements
    even = ifft(x[0::2])
    odd = ifft(x[1::2])

    # Combine
    T = [cmath.exp(2j * cmath.pi * k / N) * odd[k] for k in range(N // 2)]
    return [even[k] + T[k] for k in range(N // 2)] + [
        even[k] - T[k] for k in range(N // 2)
    ]


def next_power_of_two(x):
    return 1 if x == 0 else 2 ** (x - 1).bit_length()


def round_and_convert(C, n, m):
    # 初始化结果列表
    result = []

    # 遍历复数数组 C
    for c in C:
        # 取复数的实部，四舍五入并转换为整数
        rounded_int = float(f"{c.real:.4f}")
        # 将结果添加到列表中
        result.append(rounded_int)

    # 返回结果列表的前n+m-1个元素
    return result[: n + m - 1]


def fft_multiply(A, B):
    # 确定输出多项式的长度
    n = len(A)
    m = len(B)
    N = next_power_of_two(n + m - 1)

    # 将 A 和 B 扩展到长度 N
    A.extend([0] * (N - n))
    B.extend([0] * (N - m))

    # 应用 FFT
    A_fft = fft(A)
    B_fft = fft(B)

    # 逐点相乘
    C_fft = [A_fft[i] * B_fft[i] for i in range(N)]

    # 应用逆 FFT
    C = [i / N for i in ifft(C_fft)]

    # 由于数值误差，结果可能是复数，取实部并四舍五入
    return round_and_convert(C, n, m)


# # 示例多项式
A = [1.5, 2, 3]  # 表示多项式 1.5 + 2x + 3x^2
B = [4, 5.3, 6]  # 表示多项式 4 + 5.3x + 6x^2
C = fft_multiply(A, B)
print(C) # [6.0, 15.95, 31.6, 27.9, 18.0]