快速傅立叶变换FFT学习笔记-Toy模板网

这篇具有很好参考价值的文章主要介绍了快速傅立叶变换FFT学习笔记。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

什么是FFT？

FFT（Fast Fourier Transformation）是离散傅氏变换（DFT）的快速算法，即快速傅氏变换。FFT使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数N越多，FFT算法计算量的节省就越显著。FFT可以将多项式乘法的复杂度从 $O(n^2)$ 降到 $O (n l o g n)$ 。

下图是FFT的整体计算流程，FFT变换的复杂度为 $O (n l o g n)$ ，FFT域上的pointwise乘法的复杂度为 $O (n)$ ，逆FFT变换的复杂度为 $O (n l o g n)$ ，总体复杂度为 $O (n l o g n)$ 。

快速傅立叶变换FFT学习笔记,隐私计算及密码学基础,傅立叶分析,信号处理

多项式的表示方法

方法1：系数表示

从多项式函数的定义，将所有系数视为系数向量，而由全部系数组成的向量 $a$ 叫做该多项式的系数表达:

$f(x)=\sum_{i=0}^{n-1}a_ix^i$

$a=(a_0 ,a_2, \dots, a_{n-1})$

举个简单的例子： $f(x)=5x^0+6x^1+7x^2$ 的系数表示为 ${5, 6, 7\}$ 。反之， ${5, 6, 7\}$ 的系数编码结果为 $f(x)=5x^0+6x^1+7x^2$ 。

系数表示特点是对多项式加法友好，时间复杂度是 $O (n)$ 。但是对多项式乘法不友好，采用多项式逐项相乘，时间复杂度仍为 $O(n^2)$ 。
注：系数表示的乘法表示卷积操作。

方法2：点值表示

任意选取 $n$ 个不同的自变量 $x$ 带入多项式函数 $f (x)$ 进行求值运算，将得到 $n$ 个不同的结果。于是，多项式的点值表达就是由这 $n$ 个数值点组成的集合：

$\{(x_0, f(x_0)), (x_1, f(x_1)), \dots, (x_{n-1}, f(x_{n-1}))\}$

举个简单的例子： $f(x)=5x^0+6x^1+7x^2$ 的点值表示为 ${(0, f(0)), (1, f(1)), (2, f(2))\}$ 。反之， ${5, 6, 7\}$ 的点值编码应该满足 $f (0) = 5, f (1) = 6, f (2) = 7$

点值表示的特点是对多项式乘法友好，时间复杂度是 $O (n)$ 【因为可以做element-wise乘法（EWMM）】，但是多项式加法不友好。

复数

复数的定义：设 $a, b$ 为实数，则 $z = a + bi$ 的数称为复数，其中 $a$ 称为实部， $b$ 称为虚部。
复数的模为： $|z|=\sqrt{a^2+b^2}$ 。一个复数的共轭复数为： $\overline{z}=a-bi$ ，即改变虚部的符号。

单位复数根
对于任意一个复数 $\omega$ ，其 $n$ 次幂的结果为1，就称复数 $\omega$ 是 $n$ 次单位复数根，即
$\omega^n=1$

可以看到， $n$ 次单位复数根有 $n$ 个，其几何意义为： $n$ 个单位复数根均匀地分布在以复平面原点为圆心的单位圆上。

快速傅立叶变换FFT学习笔记,隐私计算及密码学基础,傅立叶分析,信号处理

在几何意义的单位圆中，我们将圆周角 $2\pi$ 均分成 $n$ 份，则 $\frac{2\pi}{n}$ 叫做单位根的幅角。
由欧拉公式：
$e^{i\frac{2\pi}{n}} = cos \frac{2\pi}{n} + i sin \frac{2\pi}{n}$

定义 $\omega_n$ 表示一个 $n$ 次单位根：
$\omega_n = e^{i\frac{2\pi}{n}}$

$\omega_n$ 也是主 $n$ 次单位根，其余 $w_{n}^{1}、w_{n}^{2}$ 等叫做 $n$ 次单位根的幂次，记为：
$\omega_{n}^{k} = e^{i\frac{2k\pi}{n}}, k=0,1,\dots,n-1$

于是，很容易知道：
$\omega_n^0=\omega_n^n=1, \omega_n^{\frac{n}{2}}=-1$

单位复数根的性质1：消去引理
$\omega_{dn}^{dk} = e^{i\frac{2dk\pi}{dn}}=e^{i\frac{2k\pi}{n}}=w_{n}^{k}$

单位复数根的性质1：折半引理
$\omega_{n}^{k+\frac{n}{2}} = \omega_{n}^{k}\omega_{n}^{\frac{n}{2}}=-\omega_{n}^{k}$
于是也可以得到：
$(\omega_{n}^{k+\frac{n}{2}})^2=(-\omega_{n}^{k})^2=(\omega_{n}^{k})^2=\omega_{n}^{2k}=\omega_{\frac{n}{2}}^{k}$
好处是将 $n$ 降到了原来的一半。

DFT：离散傅立叶变换

假设多项式：
$A(x)=\sum_{i=0}^{n-1}a_ix^i$

把 $n$ 次单位根的幂次 $x=\omega_n^k$ 分别代入多项式：
$y_k = A(\omega_n^k)=\sum_{i=0}^{n-1}a_i\omega_n^{ki}, k=0,1,\dots,n-1$

记 $y=(y_0, y_1, \dots, y_{n-1})$ 是系数向量 $a=(a_0, a_1,\dots,a_{n-1})$ 的离散傅立叶变换，即DFT。

IDFT：离散傅立叶逆变换
$a_j = \frac{1}{n}\sum_{k=0}^{n-1}y_k\omega_n^{-ki}$

DFT对应多项式求值
IDFT对应插值，求多项式系数

值得注意的是，DFT的复杂度仍然是 $O(n^2)$ 。

FFT和蝶形计算

FFT的原理是将多项式分解成奇偶两部分，并用分治的思想依次计算下去。
$A(x)=a_0+a_1x^1+\dots+a_{n-1}x^{n-1}$

$A_0(x)=a_0+a_2x^1+\dots+a_{n-2}x^{\frac{n-2}{2}}$

$A_1(x)=a_1+a_3x^1+\dots+a_{n-1}x^{\frac{n-2}{2}}$

$A(x)=A_0(x^2)+xA_1(x^2)$

证明：
$A(x)=A_0(x^2)+xA_1(x^2)=a_0+a_2x^2+a_4x^4+\dots+a_{n-2}x^{n-2}+\\a_1x^1+a_3x^3+a_5x^5+\dots+a_{n-1}x^{n-1}$
得证！

将 $x=\omega_n^k$ 代入 $A(x)=A_0(x^2)+xA_1(x^2)$ 中，得到：
$A(\omega_n^k)=A_0(\omega_n^{2k})+\omega_n^kA_1(\omega_n^{2k})=A_0(\omega_{\frac{n}{2}}^{k})+\omega_n^kA_1(\omega_{\frac{n}{2}}^{k})$

将 $x=\omega_n^{k+\frac{n}{2}}$ 代入 $A(x)=A_0(x^2)+xA_1(x^2)$ 中，得到：
$A(\omega_n^{k+\frac{n}{2}})=A_0(\omega_n^{2k+n}) + \omega_n^{k+\frac{n}{2}}A_1(\omega_n^{2k+n})=A_0(\omega_n^{2k})-\omega_n^kA_1(w_n^{2k})=A_0(\omega_{\frac{n}{2}}^{k})-\omega_n^kA_1(w_{\frac{n}{2}}^{k})$

我们发现， $A(\omega_n^k)$ 和 $A(\omega_n^{k+\frac{n}{2}})$ 的第一项完全相同，仅第二项为相反数。因此，如果知道 $A_0(\omega^k_{\frac{n}{2}})$ 和 $A_1(\omega^k_{\frac{n}{2}})$ 的值，我们就可以同时知道 $A(\omega^k_n)$ 和 $A(\omega^{k+{n\over 2}}_n)$ ，所以可以用分治思想计算FFT，原问题的规模缩减了一半。

总结一下，FFT的计算如下：
$A(\omega_n^k)=A_0(\omega_{\frac{n}{2}}^{k})+\omega_n^kA_1(\omega_{\frac{n}{2}}^{k})$

$A(\omega_n^{k+\frac{n}{2}})=A_0(\omega_{\frac{n}{2}}^{k})-\omega_n^kA_1(w_{\frac{n}{2}}^{k})$

可以通过这样的方式将一个多项式一直分解下去，如下图是对16点输入的分解：

快速傅立叶变换FFT学习笔记,隐私计算及密码学基础,傅立叶分析,信号处理

在计算FFT时，需要成对的点做蝶形运算，这里成对的点就是0和8、4和12等，这个分组的过程可以用bit reverse实现。

8点FFT计算图示：

快速傅立叶变换FFT学习笔记,隐私计算及密码学基础,傅立叶分析,信号处理

每一对数的蝶形运算：

快速傅立叶变换FFT学习笔记,隐私计算及密码学基础,傅立叶分析,信号处理

Bit Reverse确定蝶形运算对

从下面这个8点FFT可以很清楚地看到，FFT蝶形运算时打乱了输入的顺序（倒位序），倒位序是由bit reverse操作得到的。

快速傅立叶变换FFT学习笔记,隐私计算及密码学基础,傅立叶分析,信号处理

FFT的输入为倒位序，输出为自然顺序。

快速傅立叶变换FFT学习笔记,隐私计算及密码学基础,傅立叶分析,信号处理

Bit reverse的原理其实并不复杂，从上文中16点输入的奇偶分解那个图就很容易看出来。

用PyTorch实现bit reverse

import numpy as np
import torch


def bit_reverse(N, data):
    num_bits = int(np.log2(N))

    reverse_indices = torch.zeros(N, dtype=int)
    for i in range(N):
        reverse_indices[i] = int(f"{i:0{num_bits}b}"[::-1], 2)

    reverse_indices = reverse_indices.unsqueeze(0)
    reversed_data = torch.take_along_dim(data, reverse_indices, dim=-1)  
    # 对于一个多维的输入数据，只在dim=-1，也就是N维度上进行排列

    return reverse_indices, reversed_data


N = 8
# data = torch.tensor([10, 11, 12, 13, 14, 15, 16, 17])
data = torch.tensor([[10, 11, 12, 13, 14, 15, 16, 17],
                     [20, 21, 22, 23, 24, 25, 26, 27]])
reverse_indices, reversed_data = bit_reverse(N, data)
print('original data:\n', data)
print('reverse indices:', reverse_indices)
print('reversed data:\n', reversed_data)

_, reversed_data_back = bit_reverse(N, reversed_data)
print('reverse back:\n', reversed_data_back)

original data:
 tensor([[10, 11, 12, 13, 14, 15, 16, 17],
        [20, 21, 22, 23, 24, 25, 26, 27]])
reverse indices: tensor([[0, 4, 2, 6, 1, 5, 3, 7]])
reversed data:
 tensor([[10, 14, 12, 16, 11, 15, 13, 17],
        [20, 24, 22, 26, 21, 25, 23, 27]])
reverse back:
 tensor([[10, 11, 12, 13, 14, 15, 16, 17],
        [20, 21, 22, 23, 24, 25, 26, 27]])

RFFT和FFT

RFFT中的R是实数的意思，RFFT是FFT的特殊版本，为实数输入设计。RFFT利用了实数的傅立叶变换为共轭对称这个事实，因此RFFT只需要计算一半的傅立叶变换系数。所以RFFT效率明显高于FFT，并且也只有一半的存储开销。

快速傅立叶变换FFT学习笔记,隐私计算及密码学基础,傅立叶分析,信号处理

因此，当我们的输入为实数时（比如图像卷积任务），我们就可以利用实数的傅立叶变换为共轭对称这个特性，用RFFT替换FFT来提高计算效率。

为什么以实数为输入的FFT结果是共轭对称？

共轭对称示意图如下：

快速傅立叶变换FFT学习笔记,隐私计算及密码学基础,傅立叶分析,信号处理

可以看到以实数为输入的FFT的显著特征：在FFT域上，X[0]和X[N/2]都是实数（虚部为0），其他部分均为复数。复数部分成共轭对称，即X[i]和X[N-i]是共轭对称的。

我们直接看DFT的计算公式（把上文中的索引i改成了t，方便和复数i区分开）：
$A(\omega_n^k)=\sum_{t=0}^{n-1}a_t\omega_n^{kt}$

其中,
$\omega_{n}^{k} = e^{i\frac{2k\pi}{n}}$

于是，代入得到：
$A(\omega_n^k)=\sum_{t=0}^{n-1}a_t e^{i\frac{2kt\pi}{n}}~~~~(1)$

同时，我们可以计算出与上面点对称的点：
$\omega_n^{n-k}=e^{i\frac{2(n-k)\pi}{n}}$

同样，代入得到：
$A(\omega_n^{n-k})=\sum_{t=0}^{n-1}a_t \omega_n^{(n-k)t}$

其中，
$\omega_n^{(n-k)t}=\omega_n^{nt-kt}=\omega_n^{nt}/\omega_n^{kt}=\omega_n^{-kt}$

于是，
$A(\omega_n^{n-k})=\sum_{t=0}^{n-1}a_t \omega_n^{-kt}=\sum_{t=0}^{n-1}a_t e^{-i\frac{2kt\pi}{n}}~~~~(2)$

容易发现，式（1）和（2）只是在 $e$ 的指数上为相反数关系！
根据欧拉公式，对于（1）：
$e^{i\frac{2kt\pi}{n}}=cos\frac{2kt\pi}{n}+isin\frac{2kt\pi}{n}$

对于（2）：
$e^{-i\frac{2kt\pi}{n}}=cos\frac{2kt\pi}{n}-isin\frac{2kt\pi}{n}$

证毕！

用PyTorch实现FFT和RFFT

import torch
torch.manual_seed(55)

N = 8
x = torch.randint(0, 5, (1, N))
print('x:', x)

x_fft = torch.fft.fft(x)  # 中心共轭对称的fft值
print('x_fft:', x_fft)

x_ifft = torch.fft.ifft(x_fft)
print('x_ifft:', x_ifft)

x_rfft = torch.fft.rfft(x)  # 只存前半部分的fft值
print('x_rfft:', x_rfft)

x_irfft = torch.fft.irfft(x_rfft)
print('x_irfft:', x_irfft)

x: tensor([[2, 4, 0, 2, 3, 3, 4, 2]])
x_fft: tensor([[20.0000+0.0000j, -0.2929+3.2929j,  1.0000-3.0000j, -1.7071-4.7071j,
         -2.0000+0.0000j, -1.7071+4.7071j,  1.0000+3.0000j, -0.2929-3.2929j]])
x_ifft: tensor([[2.+0.j, 4.+0.j, 0.+0.j, 2.+0.j, 3.+0.j, 3.+0.j, 4.+0.j, 2.+0.j]])
x_rfft: tensor([[20.0000+0.0000j, -0.2929+3.2929j,  1.0000-3.0000j, -1.7071-4.7071j,
         -2.0000+0.0000j]])
x_irfft: tensor([[2., 4., 0., 2., 3., 3., 4., 2.]])

在这个结果中可以观察到两点为实数、其余点共轭对称、RFFT仅存储FFT前一半结果的性质。
注意：实际上，在PyTorch的RFFT实现中，就是先计算的FFT，然后截掉后半部分对称的结果。而在实际的FFT硬件计算架构中，做了更复杂的优化，并不需要计算完整的N点FFT。文章来源地址https://www.toymoban.com/news/detail-816835.html