【清橙A1084】【FFT】快速傅里叶变换

问题描述

　　离散傅立叶变换在信号处理中扮演者重要的角色。利用傅立叶变换，可以实现信号在时域和频域之间的转换。
　　对于一个给定的长度为n=2m (m为整数) 的复数序列X0, X1, …, Xn-1，离散傅立叶变换将得到另一个长度为n的复数序列Y0, Y1, …, Yn-1。其中
　　Yi=X0+X1wi+ X2w2i+ X3w3i+…+ Xn-1w(n-1)i
　　其中w=e2πI/n=cos(2π/n)+I sin(2π/n)，称为旋转因子，其中I为虚数单位，I2= –1。
　　给定输入序列X，请输出傅立叶变换后的输出序列。

　　由于所有的复数C都可以表示成C=a+Ib的形式，即由实部和虚部的和表示，所以C可以用一个二元组(a, b)来表示，用这种方法w可表示为(cos(2π/n), sin(2π/n))。复数的计算规则如下：
　　(a1, b1)+(a2, b2)=(a1+a2, b1+b2)
　　(a1, b1)(a2, b2)=(a1, b1)*(a2, b2)=(a1*a2-b1*b2, a1*b2+a2*b1)

　　对于本题，你可以按照上面的公式直接计算，也可以按下面的方法计算。使用上面的公式的复杂度为O(n2)，可以得到一半的分数，而下面的方法的复杂度为O(n log n)，可以得到全部的分数。如果要使用上面的公式直接计算，请跳过下面的算法描述部分。

算法描述

　　注意到上式中w=e2πI/n=cos(2π/n)+I sin(2π/n)，所以wn+k=wk，这个公式将在下面的计算用用到。
　　对上面的公式进行变形，得到：
　　Yi
　　=X0 + X1wi + X2w2i + X3w3i +…+ Xn-1w(n-1)i
　　=X0 + X2w2i + X4w4i +…+ Xn-2w(n-2)i + wi(X1 + X3w2i + X5w4i +…+ Xn-1w(n-2)i)
　　=(X0 + X2φi + X4φ2i +…+ Xn-2φ(n/2-1)i) + wi(X1 + X3φi + X5φ2i +…+ Xn-1φ(n/2-1)i)
　　其中φ=w2。利用这个公式可得
　　Yi+n/2=(X0 + X2φi+n/2 + X4φ2(i+n/2) +…+ Xn-2φ(n/2-1) (i+n/2)) + wi(X1 + X3φ(i+n/2) + X5φ2(i+n/2) +…+ Xn-1φ(n/2-1) (i+n/2))
　　其中φi+n/2=w2i+n=w2i=φi。
　　所以当i<n/2时，令pi=X0 + X2φi + X4φ2i +…+ Xn-2φ(n/2-1)i，qi= X1 + X3φi + X5φ2i +…+ Xn-1φ(n/2-1)i，则Yi=pi+wiqi，Yi+n/2= pi+wi+n/2qi。
　　利用上面的公式，我们可以得到一种快速计算旋转因子为w的傅立叶变换的方法如下（快速傅立叶变换）：
　　算法Y=Fourier(X, n, w)
　　参数：X={Xi}为待变换的序列，n为序列的长度（2的整数次幂），w为旋转因子
　　结果：X的傅立叶变换Y={Yi}
　　1. 计算{X0, X2, X4, …, Xn-2}在旋转因子为φ=w2时的傅立叶变换序列{pi}
　　2. 计算{ X1, X3, X5, …, Xn-1}在旋转因子为φ=w2时的傅立叶变换序列{qi}
　　3. 对于0<=i<n，计算Yi=pi+wiqi。其中w0=(1, 0)，wi=wi-1*w，你要设置一个临时变量进行累乘而不能每次重新计算wi。
　　使用这种方法，即可在O(n log n)的时间内计算傅立叶变换。

输入格式

　　输入的第一行包含一个整数n，表示待变换的序列的长度，n是2的整数次幂。n<=16384。
　　接下来n行，每行2个实数a, b表示序列中的一个元素为(a, b)。

输出格式

　　输出n行，每行两个实数，表示经过变换后的序列。为了防止运算时的误差影响结果的输出，请将结果中的每一个实数除以n后保留两位小数。

样例输入

4
1.0 0.0
1.0 0.0
0.0 0.0
0.0 0.0

样例输出

0.50 0.00
0.25 0.25
0.00 0.00
0.25 -0.25

【分析】

也是裸题。

  1 /*
  2 宋代苏轼
  3 《临江仙·夜饮东坡醒复醉》
  4 夜饮东坡醒复醉，归来仿佛三更。家童鼻息已雷鸣。敲门都不应，倚杖听江声。
  5 长恨此身非我有，何时忘却营营。夜阑风静縠纹平。小舟从此逝，江海寄余生。
  6 */
  7 #include <cstdio>
  8 #include <cstring>
  9 #include <algorithm>
 10 #include <cmath>
 11 #include <queue>
 12 #include <vector>
 13 #include <iostream>
 14 #include <string>
 15 #include <ctime>
 16 #define LOCAL
 17 const double Pi = acos(-1.0);
 18 const int MAXN = 200000 + 10;
 19 using namespace std;
 20 typedef long long ll;
 21 struct Num {
 22    double a , b;
 23     //构造函数
 24    Num(double x = 0,double y = 0) {a = x; b = y;}
 25     //复数的三种运算
 26    Num operator + (const Num &c) {return Num(a + c.a, b + c.b);}
 27    Num operator - (const Num &c) {return Num(a - c.a, b - c.b);}
 28    Num operator * (const Num &c) {return Num(a * c.a - b * c.b, a * c.b + b * c.a);}
 29 }x1[MAXN] , x2[MAXN];
 30
 31 //注意loglen为log后的长度
 32 void change(Num *t, int len, int loglen){
 33     //蝶形变换后的序列编号
 34     for (int i = 0; i < len; i++){
 35         int x = i, k = 0;
 36         for (int j = 0; j < loglen; j++, x >>= 1) k = (k << 1) | (x & 1);
 37         if (k < i) swap(t[k], t[i]);
 38     }
 39 }
 40 //基2-FFT
 41 void FFT(Num *x, int len, int loglen){
 42     change(x, len, loglen);
 43     int t = 1;
 44     //t为长度
 45     for (int i = 0; i < loglen; i++, (t <<= 1)){
 46         int l = 0, r = l + t;
 47         while (l < len){
 48             //初始化
 49             Num a, b, wo(cos(Pi / t), sin(Pi / t)), wn(1, 0);
 50             for (int j = l; j < l + t; j++){
 51                 a = x[j];
 52                 b = x[j + t] * wn;
 53                 //蝶形计算
 54                 x[j] = a + b;
 55                 x[j + t] = a - b;
 56                 wn = wn * wo;
 57             }
 58             //注意是合并所以要走2t步
 59             l = r + t;
 60             r = l + t;
 61         }
 62     }
 63 }
 64 //离散傅里叶变换
 65 void DFT(Num *x, int len, int loglen){
 66     int t = (1<<loglen);
 67     for (int i = 0; i < loglen; i++){
 68         t >>= 1;
 69         int l = 0, r = l + t;
 70         while (l < len){
 71             Num a, b, wn(1, 0), wo(cos(Pi / t), -sin(Pi / t));
 72             for (int j = l; j < l + t; j++){
 73                 a = x[j] + x[j + t];
 74                 b = (x[j] - x[j + t]) * wn;
 75                 x[j] = a;
 76                 x[j + t] = b;
 77                 wn = wn * wo;
 78             }
 79             l = r + t;
 80             r = l + t;
 81         }
 82     }
 83     change(x, len, loglen);
 84     for (int i= 0; i < len; i++) x[i].a /= len;
 85 }
 86 int solve(char *a, char *b){
 87     int len1, len2, len, loglen;
 88     int t, over;
 89     len1 = strlen(a) << 1;
 90     len2 = strlen(b) << 1;
 91     len = 1;
 92     loglen = 0;
 93     while (len < len1) len <<= 1, loglen++;
 94     while (len < len2) len <<= 1, loglen++;
 95     //处理len1
 96     for (int i = 0; i < len; i++){
 97         if (a[i]) x1[i].a = a[i] - ‘0‘, x1[i].b = 0;
 98         else x1[i].a = x1[i].b = 0;
 99     }
100     for (int i = 0; i < len; i++){
101         if (b[i]) x2[i].a = b[i] - ‘0‘, x2[i].b = 0;
102         else x2[i].a = x2[i].b = 0;
103     }
104     FFT(x1, len, loglen);
105     FFT(x2, len, loglen);
106     for (int i = 0; i < len; i++) x1[i] = x1[i] * x2[i];
107
108     DFT(x1, len, loglen);
109     over = len = 0;
110     //转换成十进制的整数
111     for (int i = ((len1 + len2) / 2) - 2; i >= 0; i--){
112         t = (int)((double)x1[i].a + (double)over + 0.5);
113         a[len++] = t % 10;
114         over = t / 10;
115     }
116     while (over){
117         a[len++] = over % 10;
118         over /= 10;
119     }
120     return len;
121 }
122 //输出
123 void print(char *str, int len){
124      for(len--; len>=0 && !str[len];len--);
125     if(len < 0) putchar(‘0‘);
126     else for(;len>=0;len--) putchar(str[len]+‘0‘);
127     printf("\n");
128 }
129 char a[MAXN] , b[MAXN];
130
131 int main() {
132     int n;
133
134     scanf("%d", &n);
135     for (int i = 0; i < n; i++) scanf("%lf%lf", &x1[i].a, &x1[i].b);
136     int t = 0;
137     while ((1<<t) < n) t++;
138     FFT(x1, n, t);
139     for (int i = 0; i < n; i++) printf("%.2lf %.2lf\n", x1[i].a / n, x1[i].b / n);
140    return 0;
141 }

时间： 2024-10-13 18:57:58

【清橙A1084】【FFT】快速傅里叶变换的相关文章

准零基础搞懂FFT快速傅里叶变换及其实现程序（二）

上一篇文章我们了解了DFT的原理,FFT是基于DFT的更适合计算机运算的算法,本文我们就正式开始学习FFT的原理. 首先我么先来宏观的看一下FFT.如果我们把整个FFT的算法看成一个黑盒子的话,那么它的输入就是时间波形信号,比如声音波形(横轴为时间,纵轴为振幅).外什么FFT要比DFT速度更快呢?下面(图1)解释了FFT和DFT的(对于计算机的)算法复杂度 图1 从上面的数学表达式可以看出,一个1024采样点的FFT比DFT块了102.4倍.如果傅里叶变换的数量级更大,FFT的速度优势会更明显.

FFT —— 快速傅里叶变换

问题: 已知A[], B[], 求C[],使: 定义C是A,B的卷积,例如多项式乘法等. 朴素做法是按照定义枚举i和j,但这样时间复杂度是O(n2). 能不能使时间复杂度降下来呢? 点值表示法: 我们把A,B,C看作表达式. 即: A(x)=a0 + a1* x + a2 * x2 +... 将A={(x1,A(x1)), (x2,A(x2)), (x3,A(x3))...}叫做A的点值表示法. 那么使用点值表示法做多项式乘法就很简单了:对应项相乘. 那么,如何将A和B转换成点值表示法,再将C转

【Delphi】如何在三轴加速器的频谱分析中使用FFT(快速傅里叶变换)算法

关于傅里叶变换的作用,网上说的太过学术化,且都在说原理,已经如何编码实现,可能很多人有个模糊影响,在人工智能,图像识别,运动分析,机器学习等中,频谱分析成为了必备的手段,可将离散信号量转换为数字信息进行归类分析. 今天这里将的不是如何实现,而是如何使用傅里叶变换 但频谱分析中,涉及到的信号处理知识对大部分软件开发的人来说,太过于晦涩难懂,傅里叶变换,拉普拉斯,卷积,模相,实数,虚数,复数,三角函数等等,已经能让软件工程师望而却步,造成懂知识的人无法开发,懂开发的人无法分析,而同时具备两种技能的人

FFT快速傅里叶变换

摘自:https://www.cnblogs.com/RabbitHu/p/FFT.html 快速傅里叶变换(FFT)是一种能在O(nlogn)O(nlog?n)的时间内将一个多项式转换成它的点值表示的算法. 点值表示:设A(x)是一个n−1次多项式,那么把n个不同的x代入,会得到n个y.这n对(x,y)唯一确定了该多项.由多项式可以求出其点值表示,而由点值表示也可以求出多项式. 设有两个n−1次多项式A(x)和B(x)),我们的目标是——把它们乘起来.普通的多项式乘法是O(n^2),但有趣的是

【bzoj2179】FFT快速傅里叶变换（优化高精度乘法）

#include<bits/stdc++.h> using namespace std; #define pi acos(-1) typedef complex<double> C; const int N=201100; int n,m,l,r[N],ans[N]; C a[N],b[N]; char s[N],t[N]; void fft(C *a,int f){ for(int i=0;i<n;++i) if(r[i]>i) swap(a[i],a[r[i]]);

【模板】FFT快速傅里叶变换

1 struct Complex{ 2 double x, y; 3 inline Complex(double xx=0, double yy=0){ 4 x=xx; y=yy; 5 } 6 inline Complex operator + (Complex a){ 7 return Complex(x+a.x, y+a.y); 8 } 9 inline Complex operator - (Complex a){ 10 return Complex(x-a.x, y-a.y); 11 }

浅谈FFT(快速傅里叶变换)

本文主要简单写写自己学习FFT的经历以及一些自己的理解和想法. FFT的介绍以及入门就不赘述了,网上有许多相关的资料,入门的话推荐这篇博客:FFT(最详细最通俗的入门手册),里面介绍得很详细. 为什么要学习FFT呢?因为FFT能将多项式乘法的时间复杂度由朴素的$O(n^2)$降到$O(nlogn)$,这相当于能将任意形如$f[k]=\sum\limits _{i+j=k}f[i]*f[j]$的转移方程的计算在$O(nlogn)$的时间内完成.因此对于想要进阶dp的同学来说,FFT是必须掌握的技能

FFT 快速傅里叶变换

这个东西很神奇,看了半天网上的解释和课件,研究了很长时间,算是大概明白了它的原理. 话不多说先上图. 我们要求的h(x)=f(x)*g(x),f(x)=Σai*x^i,g(x)=Σbi*x^i. 朴素求复杂度是n2的,但一个x次多项式在平面上可以由x+1个点唯一插值表示,所以我们可以先用求出x+1个点(xi,f(xi))和(xi,g(xi)),再求出(xi,f(xi)*g(xi)),就可以反解出    h(x)的表达式. 那么我们需要在nlogn的时间内干完这两步,首先xi的取值需要特殊取,令x

快速傅里叶变换FFT

快速傅里叶变换FFT DFT是信号分析与处理中的一种重要变换.但直接计算DFT的计算量与变换区间长度N的平方成正比,当N较大时,计算量太大,直接用DFT算法进行谱分析和信号的实时处理是不切实际的. 1.直接计算DFT 长度为N的有限长序列x(n)的DFT为: 2.减少运算量的思路和方法 思路:N点DFT的复乘次数等于N2.把N点DFT分解为几个较短的DFT,可使乘法次数大大减少.另外,旋转因子WmN具有周期性和对称性. (考虑x(n)为复数序列的一般情况,对某一个k值,直接按上式计算X(k)值需