zoukankan      html  css  js  c++  java
  • 【清橙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=w2ii
      所以当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("
    ");
    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
    ", x1[i].a / n, x1[i].b / n);
    140    return 0;
    141 }
    View Code
  • 相关阅读:
    java中native的用法
    用uWSGI和Nginx部署Flask项目
    elasticsearch之使用Python批量写入数据
    mysql 远程访问
    Chrome扩展及应用开发
    jQuery ajax
    Chrome扩展及应用开发-储存数据
    Chrome扩展及应用开发-扩展页面间的通信
    Python3 将本地时间转换成指定时区时间
    python如何编译py文件生成pyc、pyo、pyd以及如何和C语言结合使用
  • 原文地址:https://www.cnblogs.com/hoskey/p/4377054.html
Copyright © 2011-2022 走看看