很详细的FFT(快速傅里叶变换)概念与实现

FFT

首先要说明一个误区,很多人认为FFT只是用来处理多项式乘的,其实FFT是用来实现多项式的系数表示法和点值表示法的快速转换的,所以FFT的用处远不止多项式乘。

FFT的前置知识:点值表示法,复数运算,三角函数。

多项式的系数表示法和点值表示法

系数表示法

[A(x)=\sum_{i=0}^n a_i*x^i ]

点值表示法

不妨将A视为关于x的函数,点值表示法就是在A的图像上取n个点,则该多项式可以被这n个点唯一确定。

[a_i=A(x_i)\ A(x)={(x_1,a_1),(x_2,a_2),\dots,(x_n,a_n)} ]

点值表示法有什么好处呢?

我们知道系数表示法下两多项式相乘是(O(n^2)),但在点值表示法下奇迹出现了:

[A(x)={(x_1,a_1),(x_2,a_2),\dots,(x_n,a_n)} \ B(x)={(x_1,b_1),(x_2,b_2),\dots,(x_n,b_n)} \ A(x)B(x)={(x_1,a_1b_1),(x_2,a_2b_2),\dots,(x_n,a_nb_n)} ]

显然这个是可以O(n)实现的。虽然但是,我们几乎不会在计算中用到点值表示法,但这也给了我们一个解决多项式乘的思路。系数转点值,相乘,点值转系数。

又很显然,我们可以随便取n个数往函数里带,可惜这样又使复杂度回到了(O(n^2))。

于是FFT出现了,FFT使我们可以用(O(n\log n))的复杂度将系数转换成一组特殊的点值,并再把点值转回系数。

复数的计算

简单的理解,复数就是实数加虚数,多少都知道点虚数吧,没错,知道点就够了。

[i=\sqrt{-1}\ z_1=a+bi\ z_2=c+di\ z_1+z_2=a+c+(b+d)i\ z_1-z_2=a-c+(b-d)i\ z_1*z_2=ac-bd+(da+bc)i ]

单位根

百度一下

记住以下性质(感兴趣可以自己推推,就是基础的三角函数)

[\omega_n^k=cosk\frac{2\pi}{n}+sink\frac{2\pi}{n} \ \omega_n^0=\omega_n^n=1\ \omega_{2n}^{2k}=\omega_n^k ]

好了,现在你已经掌握了所有FFT的前置知识了,自己来推推FFT吧。

正式开始FFT

将(\omega_n^0,\dots,\omega_n^{n-1}) 这n个数带入得到点值表示,于是:

[A(x)=a_0+a_1x+a_2{x^2}+a_3{x^3}+a_4{x^4}+a_5{x^5}+ \dots+a_{n-2}x^{n-2}+a_{n-1}x^{n-1}\ A(x)=(a_0+a_2{x^2}+a_4{x^4}+\dots+a_{n-2}x^{n-2})+(a_1x+a_3{x^3}+a_5{x^5}+ \dots+a_{n-1}x^{n-1})\ A_1(x)=a_0+a_2{x}+a_4{x^2}+\dots+a_{n-2}x^{\frac{n}{2}-1}\ A_2(x)=a_1+a_3{x}+a_5{x^2}+ \dots+a_{n-1}x^{\frac{n}{2}-1}\ A(x)=A_1(x^2)+xA_2(x^2)\ ]

我们将(ωkn(k

[A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})\ A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}(\omega_n^{2k+n})=A_1(\omega_n^{2k}\omega_n^n)-\omega_n^kA_2(\omega_n^{2k}\omega_n^n)=A_1(\omega_n^{2k})-\omega_n^kA_2(\omega_n^{2k}) ]

摘自attak的blog

观察这两个式子

[A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})\\ A(\omega_n^{k+\frac{n}{2}})=A_1(\omega_n^{2k})-\omega_n^kA_2(\omega_n^{2k})\ ]

显然(怎么又是显然)只要求出(A_1(\omega_n^{2k}))和(\omega_n^kA_2(\omega_n^{2k}))就可以得出(A(\omega_n^k))和(A(\omega_n^{k+\frac{n}{2}})),而(A_1(\omega_n^{2k}))和(\omega_n^kA_2(\omega_n^{2k}))又可以进一步递归,正是这一点性质使FFT的复杂度达到了优秀的(O(n\log n))

IFFT

前面提到了,系数转点值,相乘,点值转系数,所以我们还需要能在(O(n\log n))复杂度下完成点值转系数,这就是IFFT,快速傅利叶逆变换。

可以将傅里叶变换和傅立叶逆变换的公式表示写出来(我也不会推逆变换)。

A(x) 表示多项式的系数表示法 B(x)表示多项式的点值表示法

[FFT\ B(x)=\sum_{i=0}^{N-1} A(k)\omega_N^{ik}\ IFFT\ A(x)=\frac 1 N \sum _{i=0} ^{N-1} B(k)\omega_n^{-ik}\ ]

IFFT就是在把FFT的(\omega_n^{ik}改成\omega_n^{-ik}),然后再乘个(\frac 1 N)

考虑(\omega_n^k)的几何意义(高一三角函数)可以得到

[\omega_n^{ik}=a+bi\ \omega_n^{-ik}=a-bi ]

所以IFFT只需要在FFT上做一点改动。

千万别忘了最后乘个(\frac 1 N)

实现

个人认为这是所有算法最难也最重要的部分,然而很多blog都是将这部分一笔带过,所以我决定来详细的讲讲。

递归写法

竟然是递归,那必然有递归极限。每次递归多项式的项数剩下一半,只剩一项时(A(x)=a_1)与(x)无关,所以直接返回就好了。

在求出了(A_1(\omega_n^{2k}))和(\omega_n^kA_2(\omega_n^{2k}))后做两次多项式加减可以得出(A(\omega_n^k))和(A(\omega_n^{k+\frac{n}{2}}))

由于多次的递归和数组新建和赋值使递归写法的常数大的出奇,所以我们需要更好的写法。

重难点来了!

递推做法

首先观察这张图。

很详细的FFT(快速傅里叶变换)概念与实现

我们的递归做法就是从上向下将原数列对半拆,再合并。

我们的合并顺序就是从下到上合并,详细的说先合并(a_0和a_4,a_2和a_6,a1和a_5,a_3和a_7,然后a_0,a_4合并好的整体与a_2和a_6合并好的整体再合并\dots)

观察最上和最下层的数列的二进制表示,发现就是将二进制翻转了。我们有这个神仙操作可以快速的翻转二进制。

rev[i]=(rev[i>>1]>>1)|((i&1)<

于是我们可以(O(n))得到最下层的数列,再依次往上推,即不用递归也不用开大量数组,让代码变得飞快!

void fft(cp *a,int n,int inv)
{
    int bit=0;
    while((1<>1]>>1)|((i&1)<

luogu FFT板子

#include
using namespace std;
typedef complex cp;
const int N=(1<>n>>m;
    for(int i=0;i>a[i];
    for(int i=0;i>b[i];
    for(lim=1;lim>1]>>1)|((i&1)<

继续学习NTT

Original: https://www.cnblogs.com/29taorz/p/15712742.html
Author: T_X蒻
Title: 很详细的FFT(快速傅里叶变换)概念与实现

原创文章受到原创版权保护。转载请注明出处:https://www.johngo689.com/644121/

转载文章受原作者版权保护。转载请注明原作者出处!

(0)

大家都在看

  • Locust的使用二

    Locust的使用一 通过命令参数可以配置Locust运行方式 文档 https://docs.locust.io/en/stable/configuration.html#com…

    技术杂谈 2023年5月31日
    072
  • volatile关键字的用法

    volatile关键字 什么是可见性? 可见性是指线程A改变变量的值后,线程B可以马上看到更改后变量的值 volatile的作用 关键字volatile提示线程每次从共享内存中读取…

    技术杂谈 2023年7月24日
    060
  • python 函数闭包

    def power(exp): def exp_of(num): return num ** exp return exp_of # 返回函数对象 square = power(2…

    技术杂谈 2023年7月25日
    065
  • go 单元测试testify

    testify用go实现的一个assert风格的测试框架,这个包提供了我们需要的断言的功能,提供了非常丰富的断言方法。 提供了测试suite、断言、mock三种功能。 官方文档:h…

    技术杂谈 2023年5月31日
    073
  • [深度学习]-Dataset数据集加载

    加载数据集dataloader from torch.utils.data import DataLoader form &#x81EA;&#x5DF1;&…

    技术杂谈 2023年7月10日
    085
  • Mysql之Binlog

    1、简述 binlog 二进制日志文件,这个文件记录了MySQL所有的DML操作。通过binlog日志我们可以做数据恢复,增量备份,主主复制和主从复制等等。 2、Docker中无法…

    技术杂谈 2023年6月21日
    0163
  • python动态网站爬虫实战(requests+xpath+demjson+redis)

    404. 抱歉,您访问的资源不存在。 可能是网址有误,或者对应的内容被删除,或者处于私有状态。 代码改变世界,联系邮箱 contact@cnblogs.com 园子的商业化努力-困…

    技术杂谈 2023年6月21日
    073
  • JPA作持久层操作

    JPA(Hibernate是jpa的实现) jpa是对实体类操作,从而通过封装好的接口直接设置数据库的表结构。虽然jpa可以直接通过编写java代码来操作数据库表结构,避免了sql…

    技术杂谈 2023年7月11日
    065
  • pip下载慢问题解决方案

    在使用Python开发过程中,经常要用pip安装软件包,但是直接使用pip install packagename经常慢得要死,而且慢就算了很多时候还下载完成安装失败。 pip默认…

    技术杂谈 2023年7月11日
    064
  • Spring mvc源码分析系列–Servlet的前世今生

    Spring mvc源码分析系列–Servlet的前世今生 概述 上一篇文章Spring mvc源码分析系列–前言挖了坑,但是由于最近需求繁忙,一直没有时间…

    技术杂谈 2023年7月25日
    093
  • 【原创】k8s微服务滚动发布(服务持续可用)实践笔记

    背景 对于业务和产品来讲,随时都有紧急小版本功能上线,对于研发人员来讲,线上如果有一些紧急的bug,都需要随时发版修正;而对于产品使用用户来讲,任何的功能和版本发布,要尽可能对用户…

    技术杂谈 2023年7月24日
    082
  • edraw max for mac(亿图图示 mac)中文版

    Original: https://www.cnblogs.com/123ccy/p/16548105.htmlAuthor: -Mac123-Title: edraw max f…

    技术杂谈 2023年5月31日
    0103
  • [转]私有笔记部署

    故事的起源是一个由于线性代数期末考几道计算题卡住算不出来折腾半天而考后看某课代表提前交卷又感觉人均 AK 了以致十分 emo 想要暂时逃避学习的下午。 TL;DR 思源笔记最好。快…

    技术杂谈 2023年5月30日
    0127
  • 【LEETCODE】70、字符匹配1023 Camelcase Matching

    最近做leetcode总感觉自己是个智障,基本很少有题能自己独立做出来,都是百度。。。 不过终于还是做出了一题。。。而且速度效率还可以 哎,加油吧,尽量锤炼自己 package y…

    技术杂谈 2023年7月24日
    067
  • 脚本之美│VBS 入门交互实战

    404. 抱歉,您访问的资源不存在。 可能是网址有误,或者对应的内容被删除,或者处于私有状态。 代码改变世界,联系邮箱 contact@cnblogs.com 园子的商业化努力-困…

    技术杂谈 2023年7月11日
    082
  • 下载Wistia视频

    视频上右键后会出现类似下图的对话框的视频就属于本日记讨论的对象。 右键会有About Wistia2. Copy link and thumbnail 点击,粘贴到记事本。会得到类…

    技术杂谈 2023年5月31日
    087
亲爱的 Coder【最近整理,可免费获取】👉 最新必读书单  | 👏 面试题下载  | 🌎 免费的AI知识星球