首页 > 其他分享 >快速傅立叶变换

快速傅立叶变换

时间：2024-02-27 18:03:30浏览次数：20  
标签：frac 变换 omega rev len int 傅立叶 n2 快速

快速傅立叶变换

1 引入

现在有两个多项式 $f(x)$，$g(x)$：
$$
f(x)=5x^2+3x+7\
g(x)=7x^2+2x+1
$$
我们要求出两者相乘的结果，按照多项式相乘的运算法则，把每一项相乘，总复杂度为 $O(n^2)$。

这是传统算法能到达的最好复杂度。能不能进行优化呢？

使用快速傅立叶变换，我们可以实现 $O(n\log n)$ 的复杂度。

1.1 概念

快速傅里叶变换（FFT），即运用计算机计算离散傅立叶变换 （DFT）的快速算法的统称。于 1965 年提出。

2 多项式表示

多项式的常见表示方法有两种：系数表示和点值表示。

系数表示就是最常见的表示方式，例如上面的 $f(x)$，可以表示为 $(5,3,7)$。

而点值表示法，则是在平面上取函数上的点。有一个定理：

$n$ 次多项式可以有 $n+1$ 个点唯一确定。

于是上面的 $f(x)$ 也可表示为 ${(0,7),(1,15),(-1,9)}$。

而利用点值表示求多项式乘积很方便，取相同的 $x$，然后将多项式的值相乘，得到新的点值，重复 $n+m$ 次可以得到最终多项式的点值表示。

我们发现，点值表示是 $O(n)$ 的，非常优秀。如果我们能在较短时间内将系数表示转化为点值表示，求出结果在转化回系数表示，就能快速求出多项式乘法。

3 单位根

关于复数：虚数与复数与欧拉公式 - 知乎 (zhihu.com)

以复平面上单位圆为起点，单位圆的 $n$ 等分点为终点，可以唯一得到 $n$ 个向量，也就是 $n$ 个复数。设幅角为正数且最小的复数为 $\omega_n$，称其为 $n$ 次单位根。则有：
$$
\omega_n=\cos\frac{2\pi}{n}+i\sin\frac{2\pi}{n}
$$
又由欧拉公式得：
$$
\omega_n^k=\cos\frac{2k\pi}{n}+i\sin\frac{2k\pi}{n}
$$
特别的，$\omega_n0=\omega_nn=1$。

4 快速傅立叶变换

4.1 分治法实现

FFT 的基本想法是分治。对于 DFT 来说，它分治的求解 $x=\omega_n^k$ 时 $f(x)$ 的值。

对于多项式 $f(x)=\sum\limits_{i=0}{n-1}a_ixi$，同时不妨设 $n=2^k$（缺失的部分用 $a_i=0$ 补齐）。我们将其按 $a_i$ 下标的奇偶性分开，也就是：
$$
f(x)&=(a_0+a_2x2+a_4x4+\cdots+a_{n-2}x{n-2})+\&(a_1x+a_3x3+a_5x5+\cdots+a_{n-1}x)
$$
令
$$
f_1(x)=a_0+a_2x1+a_4x2+\cdots+a_{n-2}x^{\frac n2-1}\
f_2(x)=a_1+a_3x1+a_5x2+\cdots+a_{n-1}x^{\frac n2-1}
$$
则有：
$$
f(x)=f_1(x2)+xf_2(x2)
$$
令 $x=\omega_n^k$，利用偶次单位根性质 $\omega_ni=-\omega_n$，以及 $f_1(x)$ 与 $f_2(x)$ 都为偶函数，可以知道当 $x=\omega_n^k$ 与 $x=\omega_n^{i+\frac n2}$ 对应值相同，于是有：
$$
f(\omega_nk)&=&f_1((\omega_nk)2)+\omega_nk\times f_2((\omega_nk)2)\
&=&f_1(\omega_n{2k})+\omega_nk\times f_2(\omega_n^{2k})\
&=&f_1(\omega_{\frac n2}k)+\omega_nk\times f_2(\omega_{\frac n2}^k)
$$
同时也可以得到：
$$
f(\omega_n^{k+\frac n2})&=&f_1(\omega_n{2k+n})+\omega_n\times f_2(\omega_n^{2k+n})\
&=&f_1(\omega_n{2k})-\omega_nk\times f_2(\omega_n^{2k})\
&=&f_1(\omega_{\frac n2}k)-\omega_nk\times f_2(\omega_{\frac n2}^k)
$$
因此，在求出 $f_1(\omega_{\frac n2}^k)$ 和 $f_2(\omega_{\frac n2}^k)$，可以同时求出 $f(\omega_n^k)$ 和 $f(\omega_n^{k+\frac n2})$。

然后，我们继续对 $f_1$ 和 $f_2$ 递归求解即可。这是 FFT 的一种实现：分治 DFT。

我们发现，递归求解时，必须满足 $n=2^k$ 才能实现，这正是我们早在开头就提到的。

在带入的时候，我们要带入 $n$ 个不同的值，所以直接带入 $\omega_n0,\omega_n1,\omega_n2,\cdots,\omega_n(n=2^k)$ 总共 $2^k$ 个不同值。

代码实现上， STL 给出了复数模板 <complex>。直接食用即可。

以上就是 FFT 中对于 DFT 的介绍，完成了我们在 2 结尾处所提到的第一步：系数表示转化为点值表示。

代码：

typedef complex<double> comp;

const comp i(0, 1);
const int Maxn = 1 << 20;
const double pi = acos(-1);

comp tmp[Maxn];

void DFT(comp *f, int n) {
	if(n == 1) return ;
	for(int i = 0; i < n; i++) {
		tmp[i] = f[i];
	}
	for(int i = 0; i < n; i++) {
		if(i & 1) {//偶数放左边，奇数放右边 
			f[n / 2 + i / 2] = tmp[i];
		}
		else {
			f[i / 2] = tmp[i];
		}
	}
	comp *g = f, *h = f + n / 2;//递归求解
	DFT(g, n / 2);
	DFT(h, n / 2);//分治
    //↑：分治
    //↓：合并分治
	comp cur(1, 0), step(sin(2 * pi / n), sin(2 * pi / n));
	//当前单位根为 cur，step 为两个单位根的差。
	for(int k = 0; k < n / 2; k++) {
		tmp[k] = g[k] + cur * h[k];
		tmp[k + n / 2] = g[k] - cur * h[k];//推出的两个公式 
		cur *= step;
	} 
	for(int i = 0; i < n; i++) {
		f[i] = tmp[i];
	}
}

时间复杂度 $O(n\log n)$。

4.2 倍增法实现

我们从上面的算法继续优化。我们使用递归会耗费更多额外内存。我们在数组中模拟递归中的拆分，然后倍增进行合并。

对于拆分，使用位逆序置换；

对于合并，使用蝶形运算优化。

4.2.1 位逆序置换

以八项多项式为例，拆分过程为（括号中数字为下标）：

  • 初始序列 ${0,1,2,3,4,5,6,7}$；
  • 第一次后 ${0,2,4,6}{1,3,5,7}$；
  • 第二次后 ${0,4}{2,6}{1,5}{3,7}$；
  • 第三次后 ${0}{4}{2}{6}{1}{5}{3}{7}$。

同时我们有规律：对于原先序列，将下标用二进制表示，反转后就是最终位置的下标。例如 $1$ 是 $(001)_2$，反转后是 $(100)_2$，也就是 $4$。而最终 $1$ 确实在 $4$ 上。

这个非常巧妙的东西，就叫做位逆序置换，其实名字就是算法本身了。

（证明先咕）

位逆序置换有很简单的 $O(n \log n)$ 方式，但是有一种更好的 $O(n)$ 做法。

设 $R(x)$ 为我们二进制反转后的数。首先有 $R(0)=0$。再设 $len=2^k$，其中 $k$ 为二进制数的长度。

从小到大求解 $R(x)$。这样保证了求 $R(x)$ 时，$R(\left\lfloor\dfrac{x}{2}\right\rfloor)$ 已经求出。因此，我们把 $x$ 右移一位（除以 $2$），然后反转，再右移一位，就得到了 $x$ 除二进制个位之外，其它位的反转结果。

考虑个位的反转结果：如果个位是 $0$，那么最高位就是 $0$；如果个位是 $1$，那么最高位是 $1$，此时还要加上 $\dfrac{len}2=2^{k-1}$。综上有：
$$
R(x)=\left\lfloor \dfrac{R(\left\lfloor\dfrac{x}{2}\right\rfloor)}{2}\right\rfloor+(x\bmod{2})\times\dfrac{len}2
$$

举个例子：设 $k=5$，$len=32=(100000)_2$。此时翻转 $(11001)_2$。

  1. 右移一位，即 $(1100)_2$，补齐后是 $(01100)_2$，反转后$(00110)_2$，再右移一位得到 $(0011)_2$。
  2. 由于个位为 $1$，所以加上 $2^{k-1}=(10000)_2$，即 $(10011)_2$.

代码：

int rev[Maxn];//R(x)

void change(comp y[], int len) {
	for(int i = 0; i < len; i++) {//求翻转后数组 
		rev[i] = rev[i >> 1] >> 1;
		if(i & 1) {
			rev[i] += (len >> 1);
		}
	}
	for(int i = 0; i < len; i++) {//下标交换 
		if(i < rev[i]) {//保证只交换一次 
			swap(y[i], y[rev[i]]);
		}
	}
}

4.2.2 蝶形运算优化

已知 $f_1(\omega_{\frac n2}^k)$ 和 $f_2(\omega_{\frac n2}^k)$ 后，用下面两式求出 $f(\omega_n^k)$ 和 $f(\omega_n^{k+\frac n2})$ ：
$$
f(\omega_n^k)=f_1(\omega_{\frac n2}k)+\omega_nk\times f_2(\omega_{\frac n2}^k)\
f(\omega_n^{k+\frac n2})=f_1(\omega_{\frac n2}k)-\omega_nk\times f_2(\omega_{\frac n2}^k)
$$
使用位逆序置换后，对于给定 $n,k$：

  • $f_1(\omega_{\frac n2}^k)$ 存储在数组下标为 $k$ 的位置，$f_2(\omega_{\frac n2}^k)$ 存储在数组下标为 $k+\dfrac n2$ 的位置。
  • $f(\omega_n^k)$ 存储在数组下标为 $k$ 的位置，$f(\omega_n^{k+\frac n2})$ 存储在数组下标为 $k+\dfrac n2$ 的位置。

因此可以直接再两个数组下标处进行覆盖，不避开额外数组。这就是蝶形运算。

最后，详细讲解一下倍增法求 FFT 的具体实现方式：

  1. 令当前段长为 $s=\dfrac n2$。
  2. 同时枚举序列 ${f_1(\omega_{\frac n2}^k)}$ 的左端点 $l_1=0,2s,4s,\cdots,n-2s$ 和序列 ${f_2(\omega_{\frac n2}^k)}$ 的左端点 $l_2=s,3s,5s,\cdots,n-s$。
  3. 合并时，枚举 $k=0,1,2,\cdots,s-1$，此时 $f_1(\omega_{\frac n2}^k)$ 存储在数组下标为 $l_1+k$ 的位置，$f_2(\omega_{\frac n2}^k)$ 存储在数组下标为 $l_2+k$ 的位置。
  4. 使用蝶形运算算出 $f(\omega_n^k)$ 和 $f(\omega_n^{k+\frac n2})$，然后直接覆盖。

代码在文末有。

5 快速傅里叶逆变换

我们在上文已经知道了求解 DFT 的 FFT，成功将系数表示转化为了点值表示。

接下来考虑求解 IDFT，即把点值表示转化为系数表示。

我们从线性代数的角度理解 IDFT 的过程。

首先，DFT 是一个线性变换，理解为目标多项式为向量，左乘一个矩阵得到变换后的向量，模拟带入单位根的过程：
$$
\begin{bmatrix}
y_0\
y_1\
y_2\
y_3\
\vdots\
y_{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^1 & \omega_n^1 & \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}\times
\begin{bmatrix}
a_0\
a_1\
a_2\
a_3\
\vdots\
a_{n-1}
\end{bmatrix}
$$
现在我们已经得到左边结果了，中间的 $x$ 值在目标多项式的点值中也一一对应。所以，根据矩阵的基础知识，我们只要在狮子两边同时左乘中间的逆矩阵即可。

而这个矩阵的逆矩阵也很特殊，就是每一项取倒数，在除以变换长度 $n$，就能得到逆矩阵。

为了使计算结果为原来的倒数，根据欧拉公式，有：
$$
\dfrac{1}{\omega_k}=\omega_k{-1}=e{k}}=\cos(\dfrac{2\pi}k)+i\sin(-\dfrac{2\pi}k)
$$
因此我们可以把单位根 $\omega_k$ 取成 $e^{-\frac{2\pi i}{k}}$，这样计算结果就为原先倒数，之后唯一多的操作就只有再除以长度 $n$，其他操作过程与 DFT 完全一致。我们可以定义一个函数同时完成 DFT 和 IDFT，用一个参数判断即可。

例如前面给到的递归版 FFT，就可以写成：

typedef complex<double> comp;

const comp i(0, 1);
const int Maxn = 1 << 20;
const double pi = acos(-1);

comp tmp[Maxn];

void DFT(comp *f, int n, int rev) { //rev=1:DFT  rev=-1:IDFT
	if(n == 1) return ;
	for(int i = 0; i < n; i++) {
		tmp[i] = f[i];
	}
	for(int i = 0; i < n; i++) {
		if(i & 1) {//偶数放左边，奇数放右边 
			f[n / 2 + i / 2] = tmp[i];
		}
		else {
			f[i / 2] = tmp[i];
		}
	}
	comp *g = f, *h = f + n / 2;//递归求解
	DFT(g, n / 2, rev);
	DFT(h, n / 2, rev);//分治
    //↑：分治
    //↓：合并分治
	comp cur(1, 0), step(sin(2 * pi / n), rev * sin(2 * pi / n));
	//当前单位根为 cur，step 为两个单位根的差。
	for(int k = 0; k < n / 2; k++) {
		tmp[k] = g[k] + cur * h[k];
		tmp[k + n / 2] = g[k] - cur * h[k];//推出的两个公式 
		cur *= step;
	} 
	for(int i = 0; i < n; i++) {
		f[i] = tmp[i];
	}
}

同时，我们也可以得出非递归版的 FFT 代码：

void fft(comp y[], int len, int rev) {//rev=1:DFT rev=-1:IDFT
	change(y, len);
	for(int h = 2; h <= len; h <<= 1) {//枚举长度 len 
		comp step(cos(2 * pi / h), rev * sin(2 * pi / h));
		for(int j = 0; j < len; j += h) {//枚举 0,s,2s,3s,...,n-2s(n-n)
			comp w(1, 0);
			for(int k = j, k < j + h / 2; k++) {//直接枚举 l_1+k 
				comp u = y[k];
				comp t = w * y[k + h / 2];//l_1+k+s=l_2+k
				y[k] = u + t;
				y[k + h / 2] = u - t;//覆写
				w *= step; 
			}
		}
	}
	if(rev == -1) {//如果是 IDFT，除以长度 len 
		for(int i = 0; i < len; i++) {
			y[i].x /= len;
		}
	}
}

6 模板代码

以 P3803 【模板】多项式乘法（FFT） 为例，给出代码：

#include <bits/stdc++.h>

using namespace std;

typedef complex<double> comp;

const comp i(0, 1);
const int Maxn = 6e6 + 5;
const double pi = acos(-1);
int rev[Maxn];//R(x)

int n, m, q = 1;

comp a[Maxn], b[Maxn];

void change(comp y[], int len) {
	for(int i = 0; i < len; i++) {//求翻转后数组 
		rev[i] = rev[i >> 1] >> 1;
		if(i & 1) {
			rev[i] += (len >> 1);
		}
	}
	for(int i = 0; i < len; i++) {//下标交换 
		if(i < rev[i]) {//保证只交换一次 
			swap(y[i], y[rev[i]]);
		}
	}
}

void fft(comp y[], int len, int rev) {//rev=1:DFT rev=-1:IDFT
	change(y, len);
	for(int h = 2; h <= len; h <<= 1) {//枚举长度 len 
		comp step(cos(2 * pi / h), rev * sin(2 * pi / h));
		for(int j = 0; j < len; j += h) {//枚举 0,s,2s,3s,...,n-2s(n-n)
			comp w(1, 0);
			for(int k = j; k < j + h / 2; k++) {//直接枚举 l_1+k 
				comp u = y[k];
				comp t = w * y[k + h / 2];//l_1+k+s=l_2+k
				y[k] = u + t;
				y[k + h / 2] = u - t;//覆写
				w *= step; 
			}
		}
	}
}

int main() {
	ios::sync_with_stdio(0);
	cin >> n >> m;
	for(int i = 0; i <= n; i++) {
		double x;
		cin >> x;
		a[i].real(x);
	}
	for(int i = 0; i <= m; i++) {
		double x;
		cin >> x;
		b[i].real(x);
	}
	q = 1;
	while(q <= n + m) q <<= 1;
	fft(a, q, 1);
	fft(b, q, 1);
	for(int i = 0; i <= q; i++) {
		a[i] *= b[i];
	}
	fft(a, q, -1);        
	for(int i = 0; i <= n + m; i++) {
		cout << (int)(a[i].real() / q + 0.5) << " ";
	}                                                                                                                       
	return 0;
}

标签：frac,变换,omega,rev,len,int,傅立叶,n2,快速
From： https://www.cnblogs.com/dzbblog/p/18037444

相关文章

  • 矩阵乘法与矩阵快速幂
    矩阵乘法与矩阵快速幂1矩阵乘法1.1定义对于两个矩阵$A,B$，其中$A$大小为$n\timesm$，$B$大小为$m\timesp$，则这两个矩阵可以做乘法，得到的矩阵$C$的大小为$n\timesp$。例如：$$A=\begin{bmatrix}a_{11}&a_{12}&a_{13}\a_{21}&a_{22}&a_{23}\end{bmatrix}......
  • Hybird App开发，一种快速实现纯血鸿蒙App开发的理念
    2024年1月18日的开发者（HDC）大会上，就官宣了“纯血鸿蒙”操作系统即将于2024年3季度正式投产。与此同时，支付宝、京东、小红书、微博、高德地图、中国移动等在内的超百个头部应用都启动了鸿蒙原生应用开发，鸿蒙开发者日新增注册量已过万，同时众多985、211高校接连开设HarmonyOS相关课程......
  • Python函数每日一讲 - 简洁快速学会globals()函数
    引言在Python中，globals()函数是一个强大的工具，它允许您访问全局命名空间中的所有变量和函数。本文将深入探讨globals()函数的语法、用法以及实际应用场景，帮助大家更好地理解和使用这个函数。语句概览globals()函数的语法如下：globals()函数实例下面是globals()函数......
  • GitHub项目如何快速稳定涨星，让Star飞一会
    GitHub现在已经成了日常开发中必不可少的网站，日常工作和学习中要用到好多上面的开源项目，评价项目质量好坏的一个重要标准就是看Star和Fork的数量，如果看到个Star超过100以上的，基本上这个项目是靠谱的，如果超过1000过，那已经算是很流行了，至于一万以上的，基本上都是如雷贯耳的存在了。......
  • 像闪电般击碎天空吧——快速轻量化模型之 SDXL-Lightning
    SDXL-Lightning是一个由ByteDance（字节跳动）开发的文本到图像的生成模型，其主要贡献在于其高速的生成能力和轻量化的设计。模型的特点快速生成：SDXL-Lightning能够在几步之内生成高质量的1024px图像，这使得它在生成速度上具有优势。这种快速生成的能力使得SDXL-Lightning可以......
  • 如何在三维地球上快速拉白模以辅助建筑规划设计？
       通过以下方法可以在三维地球上快速拉白模以辅助建筑规划设计。 方法/步骤下载三维地图浏览器http://www.geosaas.com/download/map3dbrowser.exe，安装完成后桌面上出现”三维地图浏览器“图标。 2、双击桌面图标打开”三维地图浏览器“ 3、点击“要素标......
  • ES6技巧（快速解构赋值、数组去重、数组转对象）
    1.如何将a,b的值快速互换leta=1;letb=2;[a,b]=[b,a]解析：首先，我们有变量 a 被赋值为 1，变量 b 被赋值为 2。然后，[a,b]=[b,a] 这行代码实际上是将数组 [b,a] 解构赋值给了数组 [a,b]。在解构赋值过程中，右侧的数组 [b,a] 中的值会被依次赋给左侧......
  • 系统设计中的快速估算技巧
    拿到一堆数据，去做架构也好，设计也好，可行性分析也好，工程上需要的是严谨。但是也有很多场景，比如即时的问题争辩和讨论，我们往往需要的是快速、直接的估算，这样的数据显然不需要非常精确，甚至可以说它一定会非常粗略，我们的目标往往只停留在“量级”的级别，但是我们依然可以对方案有......
  • Accurately computing running variance —— 已知两个数列各自的均值和方差，如何快速
    原内容来自：https://www.johndcook.com/blog/standard_deviation/计算公式：该种计算方式可以只保存历史数据的平方和，与历史数据的和。相关前文：已知两个数列各自的均值和方差，如何快速求出两个数列拼合后的均值和方差......
  • 已知两个数列各自的均值和方差，如何快速求出两个数列拼合后的均值和方差
    问题：数列A为[1,2,3,4,5,6,7,8,9]，已知数列A的均值和方差和个数为mean_x，var_x，size_x数列B为[20,21,22,23,24,25,26,27,28,29]，已知数列B的均值和方差和个数为mean_y，var_y，size_y现在将数列A与数列B拼合为数列Z，则数列Z为[1,2,3,4,5,6,7,8,9,20,......

赞助商

阅读排行

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

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

Copyright © 2020-2023 IPS99 版权所有 IPS99