首页 > 其他分享 >快速傅里叶变换（FFT）

快速傅里叶变换（FFT）

时间：2024-10-29 23:31:28浏览次数：1  
标签：cnt ch frac 变换 傅里叶 FFT int cdots omega

快速傅里叶变换（FFT）

前言

本文为个人学习笔记，大量参考了 oi-wiki 以及其他博客的内容。

问题

记：

\[f(x) = c_0 + c_1 x + c_2 x^2 + \cdots + c_{n}x^{n} \\ g(x) = d_0 + d_1 x + d_2 x^2 + \cdots + d_{m}x^{m} \\ h(x) = f(x) \times g(x) \]

在 \(\mathcal O(n \log n)\) 内解决两个多项式乘法后的系数（即给定 \(f(x)\) 和 \(g(x)\) 的系数，要你求出 \(h(x)\) 的系数）。

分析

暴力显然是 \(\mathcal O(n^2)\) 的，优化的想法是先考虑点值表示，再考虑从点值表示转换为系数表示。

具体如下：

点值表示的意思是，你需要求出（\(\omega_n^k\) 是什么先忽略，当作是 \(n\) 个已知量即可）：

\[f(\omega_n^0), f(\omega_n^1), \cdots f(\omega_n^{n-1}) \\ g(\omega_n^0), g(\omega_n^1), \cdots g(\omega_n^{n-1}) \]

那么：

\[h(\omega_n^k) = f(\omega_n^k) \times g(\omega_n^k) \]

实际上，\(n\) 个点的点值表示法也确定了一个 \(n - 1\) 次的多项式，因此，一定存在某个算法能将点值表示法转化为系数表示（这个后面再说）。

至此，FFT 的核心思想已经说清楚了，就是考虑求出 \(f, g\) 的点值表示，那么 \(h\) 的点值表示就可以在 \(\mathcal O(n)\) 的复杂度内求出，而后再考虑从点值表示转化为系数表示。

问题一：求出某个多项式的点值表示（离散傅里叶变换 DFT）

实际上这个问题的真实含义是：怎么选取 \(\omega_n^k\) 这个已知量才能使得在一个优秀的复杂度内求出多项式的点值表示。

记 \(\omega_n^k\) 表示将复数坐标系的单位圆平均分成 \(n\) 份，从 \(x\) 轴逆时针出发的第 \(k\) 条分界箭头的复数表示。

选取这个 \(\omega_n^k\) 的原因是它有某些性质，能 " 在一个优秀的复杂度内求出多项式的点值表示 "。

性质：

\[1) \ \omega_n^k = \omega_{\frac{n}{2}}^{\frac{k}{2}} \ \ \ \ \ \ \ \\ 2) \ \omega_n^{k + \frac{n}{2}} = - \omega_n^k \]

然后开始推式子：

令

\[f_1(x) = c_0 + c_2x + \cdots + c_{n - 2}x^{\frac{n - 2}{2}} \\ f_2(x) = c_1 + c_3x + \cdots + c_{n - 1}x^{\frac{n - 2}{2}} \]

显然有

\[f(x) = f_1(x^2) + xf_2(x^2) \]

将 \(\omega _n^k(0 \leq k < \frac{n}{2})\) 代入有：

\[\begin {split} f(\omega_{n}^k) &= f_1(\omega _n^{2k}) + \omega_n^k f_2(\omega_n^{2k}) \\ &= f_1(\omega _{\frac{n}{2}}^{k}) + \omega_n^k f_2(\omega _{\frac{n}{2}}^{k}) \end {split} \]

将 \(\omega_n^{k + \frac{n}{2}}(0 \leq k < \frac{n}{2})\) 代入有：

\[\begin {split} f(\omega_n^{k + \frac{n}{2}}) &= f_1(\omega _n^{2k + n}) + \omega_n^{k + \frac{n}{2}} f_2(\omega_n^{2k + n}) \\ &= f_1(\omega _{\frac{n}{2}}^{k}) - \omega_n^k f_2(\omega _{\frac{n}{2}}^{k}) \end {split} \]

递归求解即可，有 \(\log n\) 层，时间复杂度为 \(\mathcal O(n \log n)\)，为了方便处理，一般把 \(n\) 处理为 \(\geq(n + m)\) 的二次幂，多出来的部分系数补为 \(0\) 即可。

具体实现中有以下注意事项：

1、虚数可以使用 C++ STL 库中的 complex 类型；

代码

#include <bits/stdc++.h>
template < typename T >
inline void read(T &cnt) {
    cnt = 0; char ch = getchar(); bool op = 1;
    for (; ! isdigit(ch); ch = getchar())
        if (ch == '-') op = 0;
    for (; isdigit(ch); ch = getchar())
        cnt = cnt * 10 + ch - 48;
    cnt = op ? cnt : - cnt;
}

const int N = (1 << 22) + 5;
const double PI = acos(-1);

inline void FFT(std::complex < double > *A, int n) {
    if (n == 1) return;
    int m = (n >> 1);
    std::complex < double > A0[m], A1[m];
    for (int i = 0; i < m; ++ i) {
        A0[i] = A[i * 2];
        A1[i] = A[i * 2 + 1];
    }
    FFT(A0, m); FFT(A1, m); // 递归处理
    auto W = std::complex < double > (cos(2.0 * PI / n), sin(2.0 * PI / n)),
         w = std::complex < double > (1.0, 0.0); // 从 w_n^0 出发
    for (int i = 0; i < m; ++ i) { // 根据式子计算 A 即可
        A[i] = A0[i] + w * A1[i];
        A[i + m] = A0[i] - w * A1[i];
        w *= W; // 等价于 w_n^k -> w_n^{k + 1}
    }
}

int n, m;
std::complex < double > F[N], G[N];

int main() {
    read(n), read(m);
    for (int i = 0; i <= n; ++ i) {
        int x; read(x);
        F[i] = x;
    }
    for (int i = 0; i <= m; ++ i) {
        int x; read(x);
        G[i] = x;
    }
    int sum = 1;
    while (sum <= n + m) sum *= 2; // 补齐成二次幂
    FFT(F, sum);
    FFT(G, sum);
    for (int i = 0; i < sum; ++ i)
        F[i] *= G[i];
    return 0;   
}

优化

递归实在太慢了！

以 \(8\) 项多项式为例，模拟拆分的过程：

  • 初始序列为 \(\{x_0, x_1, x_2, x_3, x_4, x_5, x_6, x_7\}\)
  • 一次二分之后 \(\{x_0, x_2, x_4, x_6\},\{x_1, x_3, x_5, x_7 \}\)
  • 两次二分之后 \(\{x_0,x_4\} \{x_2, x_6\},\{x_1, x_5\},\{x_3, x_7 \}\)
  • 三次二分之后 \(\{x_0\}\{x_4\}\{x_2\}\{x_6\}\{x_1\}\{x_5\}\{x_3\}\{x_7 \}\)

规律：其实就是原来的那个序列，每个数用二进制表示，然后把二进制翻转对称一下，就是最终那个位置的下标。比如 \(x_1\) 是 001，翻转是 100，也就是 4，而且最后那个位置确实是 4。我们称这个变换为位逆序置换（bit-reversal permutation），证明留给读者自证。

实际上，位逆序置换可以 \(\mathcal O(n)\) 从小到大递推实现，设 \(len=2^k\)，其中 k 表示二进制数的长度，设 \(R(x)\) 表示长度为 \(k\) 的二进制数 \(x\) 翻转后的数（高位补 0）。我们要求的是 \(R(0),R(1),\cdots,R(n-1)\)。

首先 \(R(0)=0\)。

我们从小到大求 \(R(x)\)。因此在求 \(R(x)\) 时，\(R\left(\left\lfloor \dfrac{x}{2} \right\rfloor\right)\) 的值是已知的。因此我们把 \(x\) 右移一位（除以 \(2\)），然后翻转，再右移一位，就得到了 \(x\) 除了（二进制）个位之外其它位的翻转结果。

考虑个位的翻转结果：如果个位是 0，翻转之后最高位就是 0。如果个位是 1，则翻转后最高位是 1，因此还要加上 \(\dfrac{len}{2}=2^{k-1}\)。综上

\[R(x)=\left\lfloor \frac{R\left(\left\lfloor \frac{x}{2} \right\rfloor\right)}{2} \right\rfloor + (x\bmod 2)\times \frac{len}{2} \]

举个例子：设 \(k=5\)，\(len=(100000)_2\)。为了翻转 \((11001)_2\)：

  1. 考虑 \((1100)_2\)，我们知道 \(R((1100)_2)=R((01100)_2)=(00110)_2\)，再右移一位就得到了 \((00011)_2\)。
  2. 考虑个位，如果是 \(1\)，它就要翻转到数的最高位，即翻转数加上 \((10000)_2=2^{k-1}\)，如果是 \(0\) 则不用更改。

蝶形运算优化

已知 \(f_1(\omega_{n/2}^k)\) 和 \(f_2(\omega_{n/2}^k)\) 后，需要使用下面两个式子求出 \(f(\omega_n^k)\) 和 \(f(\omega_n^{k+n/2})\)：

\[\begin{aligned} f(\omega_n^k) & = f_1(\omega_{n/2}^k) + \omega_n^k \times f_2(\omega_{n/2}^k) \\ f(\omega_n^{k+n/2}) & = f_1(\omega_{n/2}^k) - \omega_n^k \times f_2(\omega_{n/2}^k) \end{aligned} \]

使用位逆序置换后，对于给定的 \(n, k\)：

  • \(f_1(\omega_{n/2}^k)\) 的值存储在数组下标为 \(k\) 的位置，\(f_2(\omega_{n/2}^k)\) 的值存储在数组下标为 \(k + \dfrac{n}{2}\) 的位置。
  • \(f(\omega_n^k)\) 的值将存储在数组下标为 \(k\) 的位置，\(f(\omega_n^{k+n/2})\) 的值将存储在数组下标为 \(k + \dfrac{n}{2}\) 的位置。

因此可以直接在数组下标为 \(k\) 和 \(k + \frac{n}{2}\) 的位置进行覆写，而不用开额外的数组保存值。此方法即称为 蝶形运算，或更准确的，基 - 2 蝶形运算。

再详细说明一下如何借助蝶形运算完成所有段长度为 \(\frac{n}{2}\) 的合并操作：

1、令段长度为 \(s = \frac{n}{2}\)；
2、同时枚举序列 \(\{f_1(\omega_{n/2}^k)\}\) 的左端点 \(l_g = 0, 2s, 4s, \cdots, N-2s\) 和序列 \(\{f_2(\omega_{n/2}^k)\}\) 的左端点 \(l_h = s, 3s, 5s, \cdots, N-s\)；
3、合并两个段时，枚举 \(k = 0, 1, 2, \cdots, s-1\)，此时 \(f_1(\omega_{n/2}^k)\) 存储在数组下标为 \(l_g + k\) 的位置，\(f_2(\omega_{n/2}^k)\) 存储在数组下标为 \(l_h + k\) 的位置；
4、使用蝶形运算求出 \(f(\omega_n^k)\) 和 \(f(\omega_n^{k+n/2})\)，然后直接在原位置覆写。

代码

#include <bits/stdc++.h>
template < typename T >
inline void read(T &cnt) {
    cnt = 0; char ch = getchar(); bool op = 1;
    for (; ! isdigit(ch); ch = getchar())
        if (ch == '-') op = 0;
    for (; isdigit(ch); ch = getchar())
        cnt = cnt * 10 + ch - 48;
    cnt = op ? cnt : - cnt;
}

const int N = (1 << 22) + 5;
const double PI = acos(-1);

int rev[N];

inline void change(std::complex < double > *A, int n) {
    for (int i = 0; i < n; ++ i) { // 求 R 数组
        rev[i] = rev[i >> 1] >> 1;
        if (i & 1) {
            rev[i] |= (n >> 1);
        }
    }

    for (int i = 0; i < n; ++ i) // 将原序列 变为 底层对应的序列
        if (i < rev[i]) std::swap(A[i], A[rev[i]]);
} 
inline void FFT(std::complex < double > *A, int n) {
    change(A, n);
    for (int m = 2; m <= n; m *= 2) { // m 是当前处理的每段长度
        auto W = std::complex < double > 
          (cos(2.0 * PI / m), sin(2.0 * PI / m));  
        for (int x = 0; x < n; x += m) { // x 是每段的开头
            auto w = std::complex < double > (1.0, 0.0);
            for (int i = x; i < x + m / 2; ++ i) { // 求出每段的点值表示 根据公式求即可
                auto A0 = A[i], A1 = A[i + m / 2];
                A[i] = A0 + w * A1;
                A[i + m / 2] = A0 - w * A1;
                w *= W;
            }
        }
    }
}


int n, m;
std::complex < double > F[N], G[N];

int main() {
    change(F, 8);
    read(n), read(m);
    for (int i = 0; i <= n; ++ i) {
        int x; read(x);
        F[i] = x;
    }
    for (int i = 0; i <= m; ++ i) {
        int x; read(x);
        G[i] = x;
    }
    int sum = 1;
    while (sum <= n + m) sum *= 2; // 补齐成二次幂
    FFT(F, sum);
    FFT(G, sum);
    for (int i = 0; i < sum; ++ i)
        F[i] *= G[i];
    FFT(F, sum);
    return 0;   
}

问题二：将点值表示转化为系数表示（傅里叶反变换 IDFT）

点值表示的矩阵形式为：

\[\begin{bmatrix}f(\omega_n^0) \\ f(\omega_n^1) \\ f(\omega_n^2) \\ f(\omega_n^3) \\ \vdots \\ f(\omega_n^{n-1}) \end{bmatrix} = \begin{bmatrix}1 & 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_n^1 & \omega_n^2 & \omega_n^3 & \cdots & \omega_n^{n-1} \\ 1 & \omega_n^2 & \omega_n^4 & \omega_n^6 & \cdots & \omega_n^{2(n-1)} \\ 1 & \omega_n^3 & \omega_n^6 & \omega_n^9 & \cdots & \omega_n^{3(n-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \omega_n^{3(n-1)} & \cdots & \omega_n^{(n-1)^2} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_{n-1} \end{bmatrix} \]

怎么求系数 \(a\) 呢？根据线性代数的知识：

\[Ax = b \\ x = A^{-1}b \]

如果能求出 \(A^{-1}\)，那么 \(A^{-1} b\) 也是两个多项式相乘的结果，FFT 即可。

唯一的问题变为怎么求解 \(A^{-1}\)。

根据矩阵的逆的定义，有

\[A^{-1} \cdot A = E \]

设 \(V\) 为原矩阵，\(G\) 为逆矩阵，考虑最终落在 \(E(i, j)\) 的值：

\[E(i, j) = \sum_{k=0}^{n-1} G(i, k) \cdot V(k, j) = \sum_{k=0}^{n-1} G(i, k) \cdot \omega_n^{kj} = [i == j] \]

引理

当 \(k\) 不是 \(n\) 的倍数时，

\[\sum_{i=0}^{n-1}\omega_n^{ki} = 0 \]

证明如下：

\[\sum_{i=0}^{n-1}\omega_n^{ki} = \frac{\omega_n^{kn} - 1}{1 - \omega_n^{k}} = \frac{1 - 1}{1 - \omega_n^{k}} = 0 \]

令 \(G(i, k) = \omega_n^{-ik}\)，则：

\[\sum_{k=0}^{n-1} G(i, k) \cdot \omega_n^{kj} = \sum_{k=0}^{n-1} \omega_n^{-ik} \cdot \omega_n^{kj} = \sum_{k=0}^{n-1} \omega_n^{k(j-i)} \]

当 \(j-i\) 不为 \(n\) 的倍数（0）时，上式为 0；

反之，有：

\[\sum_{k=0}^{n-1} \omega_n^{k(j-i)} = \sum_{k=0}^{n-1} \omega_n^{0} = n \]

再前面补个系数 \(\frac{1}{n}\) 即可，故：

\[G(i, k) = \frac{1}{n}\omega_n^{-ik} \]

int main() {
    change(F, 8);
    read(n), read(m);
    for (int i = 0; i <= n; ++ i) {
        int x; read(x);
        F[i] = x;
    }
    for (int i = 0; i <= m; ++ i) {
        int x; read(x);
        G[i] = x;
    }
    int sum = 1;
    while (sum <= n + m) sum *= 2; // 补齐成二次幂
    FFT(F, sum);
    FFT(G, sum);
    for (int i = 0; i < sum; ++ i)
        F[i] *= G[i];
    FFT(F, sum);
    std::reverse(F + 1, F + sum); // 从第一位开始翻转
								  // 翻转后变为 0 1-n, 2-n, ..., -1 
    							  // 实际上等价于 0, 1, 2, ..., n-1
    for (int i = 0; i <= n + m; ++ i) { // 四舍五入
        std::cout << (int)(F[i].real() / sum + 0.5) << ' ';
    }
    return 0;   
}

标签：cnt,ch,frac,变换,傅里叶,FFT,int,cdots,omega
From： https://www.cnblogs.com/chzhc-/p/18514743

相关文章

  • 手写fft算法，和内置fft算法对比
    好的，下面我将提供一个完整的Python示例，包括手写FFT算法（快速傅里叶变换）和使用NumPy的内置FFT算法，然后对比两者的结果并绘制图形。1.手写FFT算法我们将实现一个简单的FFT算法，即Cooley-Tukey算法，这是一种最常用的FFT算法。2.NumPy内置FFT算法NumPy库提供了一个非常高效的FFT实......
  • 三维变换矩阵知识点
    一、行矩阵与列矩阵、行主序与列主序行主序和列主序影响如何在内存中访问和存储矩阵数据，与行矩阵和列矩阵的概念没有直接的关系行矩阵与列矩阵：行矩阵：        $$\begin{pmatrix}a_{11}&a_{12}\\a_{21}&a_{22}\end{pmatrix}$$列矩阵：指的是矩阵的元......
  • 不做行变换证明矩阵的行秩等于列秩
    设\(A\inM_{m\timesn}(\mathbbF)\)，令\(r=\dimR(A),s=\dimC(A)\).不妨设\(A\)的基行为前\(r\)行，令\(\tilde{A}\)为截取\(A\)的前\(r\)行所得矩阵，令\(t=\dimC(\tilde{A})\)，不妨设\(\tilde{A}\)的基列为前\(t\)列.任取\(1\lek\len\)，则存在\(\l......
  • 点跟踪论文—RAFT: Recurrent All-Pairs Field Transforms for Optical Flow-递归的全
    点目标跟踪论文—RAFT:RecurrentAll-PairsFieldTransformsforOpticalFlow-递归的全对场光流变换读论文RAFT密集光流跟踪的笔记RAFT是一种新的光流深度网络结构，由于需要基于点去做目标的跟踪，因此也是阅读了像素级别跟踪的一篇ECCV2020的经典论文——RAFT,递归的......
  • FFT在实际工程应用中有什么作用
    快速傅里叶变换（FFT）在实际工程应用中起到三大作用：1、信号处理与分析；2、图像处理领域；3、声学工程。其中，信号处理是FFT应用最广泛的领域。FFT算法能够将时域信号转换为频域信号，分析信号的频率组成，这对于去噪、信号压缩、特征提取等具有重要意义。在现代通讯系统中，FFT通过快速转换帮......
  • 拉普拉斯变换10.24
    目录1.拉普拉斯变换2.拉普拉斯收敛域3.导数的拉普拉斯变换推导过程5.传递函数6.电感电阻电路动态方程拉氏变换常数输入L逆变换7.控制系统传递函数8.非零初始状态的传递函数1.拉普拉斯变换\[\mathscr{L}[f(t)]=F(s)=\int^\infty_0f(t)e^{-st}dt\]$s=\sigma+j\ome......
  • matlab求冲击输入拉普拉斯变换，和响应
    在MATLAB中，求解冲击输入的拉普拉斯变换以及系统的冲击响应通常涉及以下几个步骤：定义冲击输入：在MATLAB中，冲击输入通常用狄拉克δ函数（Diracdeltafunction）表示，可以使用dirac函数来创建。计算拉普拉斯变换：使用laplace函数计算冲击输入的拉普拉斯变换。计算系统冲击响应：如果......
  • LTE 利用FFT 实现PSS的快速相关
    本期介绍一下怎么利用快速傅里叶变换来实现LTEPSS的快速相关。看下数字信号处理书本上线性卷积的数学表达式：假设h(n)和x(n)的长度分别为N和M，线性卷积的结果用yline(n)表示，则对比一下可以发现LTEPSS相关可以用线性卷积来实现，只需要把本地序列共轭翻转。我们知道用FFT方......
  • LTE 基于快速哈达玛hadamard变换SSS辅同步信号SSS检测之hadamard变换公式推导
    LTESSSs序列生成的阶数为31阶，所以hadamard矩阵的阶数为32阶，定义一个32阶的hadamard矩阵H32，下面进行hadamard快速变换公式推导继续分解后面公式的推导小编还在继续......
  • 基于FFT + CNN - BiGRU-Attention 时域、频域特征注意力融合的电能质量扰动识别模型
    往期精彩内容：Python-电能质量扰动信号数据介绍与分类-CSDN博客Python电能质量扰动信号分类(一)基于LSTM模型的一维信号分类-CSDN博客Python电能质量扰动信号分类(二)基于CNN模型的一维信号分类-CSDN博客Python电能质量扰动信号分类(三)基于Transformer的一维信号分类模型-......

赞助商

阅读排行

网站主页关于我们联系我们网站地图

本网站内容转载自其他媒体，侵权联系[admin##ips99.com]。

Copyright © 2020-2023 IPS99 版权所有 IPS99