zoukankan      html  css  js  c++  java
  • FFT

    A fast Fourier transform (FFT) is an algorithm that computes the discrete Fourier transform (DFT) of a sequence, or its inverse (IDFT). Fourier analysis converts a signal from its original domain (often time or space) to a representation in the frequency domain and vice versa. The DFT is obtained by decomposing a sequence of values into components of different frequencies.[1] This operation is useful in many fields, but computing it directly from the definition is often too slow to be practical. An FFT rapidly computes such transformations by factorizing the DFT matrix into a product of sparse (mostly zero) factors.[2] As a result, it manages to reduce the complexity of computing the DFT from O ( n 2 ) {displaystyle O(n^{2})} , which arises if one simply applies the definition of DFT, to O ( n log ⁡ n ) {displaystyle O(nlog n)} , where n {displaystyle n} is the data size. The difference in speed can be enormous, especially for long data sets where N may be in the thousands or millions. In the presence of round-off error, many FFT algorithms are much more accurate than evaluating the DFT definition directly. There are many different FFT algorithms based on a wide range of published theories, from simple complex-number arithmetic to group theory and number theory.

    Fast Fourier transforms are widely used for applications in engineering, science, and mathematics. The basic ideas were popularized in 1965, but some algorithms had been derived as early as 1805.[1] In 1994, Gilbert Strang described the FFT as "the most important numerical algorithm of our lifetime",[3][4] and it was included in Top 10 Algorithms of 20th Century by the IEEE journal Computing in Science & Engineering.[5] In practice, the computation time can be reduced by several orders of magnitude in such cases, and the improvement is roughly proportional to N / {displaystyle /} log N. This huge improvement made the calculation of the DFT practical; FFTs are of great importance to a wide variety of applications, from digital signal processing and solving partial differential equations to algorithms for quick multiplication of large integers.

    The best-known FFT algorithms depend upon the factorization of N, but there are FFTs with O(N log N) complexity for all N, even for prime N. Many FFT algorithms only depend on the fact that e − 2 π i / N {displaystyle e^{-2pi i/N}} is an N-th primitive root of unity, and thus can be applied to analogous transforms over any finite field, such as number-theoretic transforms. Since the inverse DFT is the same as the DFT, but with the opposite sign in the exponent and a 1/N factor, any FFT algorithm can easily be adapted for it.
    FFT is to solve the following problem.

    [P(x)=sum_{i=0}^{n}a_ix_i=a_0+a_1x+a_2x_2+⋯+a_nx_n ]

    Complex number

    A=a+ib,(a,b∈ℝ,A∈ℂ)
    $ i=sqrt{-1} ( ) A_1+A_2=(a_1+i* b_1)+(a_2+i* b_2)=(a_1+a_2)+i* (b_1+b_2) $$ A_1−A_2=(a_1+i* b_1)−(a_2+i* b_2)=(a_1−a_2)+i* (b_1−b_2) $

    $ A_1×A_2=(a_1+i* b_1)(a_2+i* b_2)=a_1a_2+i* 2b_1b_2+i* a_1b_2+i* a_2b_1=(a_1a_2−b_1b_2)+i* (a_1b_2+a_2b_1) $

    单位根

    记 $ ω_n = cos(2πn)+i* sin(2πn) $,则 n 次单位根就是

    $ ω_0n,ω_1n,⋯,ω{n−1}n $ 。 通项公式, $ ω^i_n=cos(2iπn)+i* sin(2iπn) i∈[1,n] $

    it has the following fact

    $ ω{dk}_{dn}=ωk_n $

    $ ω{k+1}_n=ωk_n* ω_n $

    $ ω^{frac{n}{2}}_n=-1 $

    $ ω{k+frac{n}{2}}_n=-ωk_n $
    DFT

    [F(x)=sum_{i=0}^{n-1}a_i* x^i ]

    [F_0(x)=sum_{i=0}^{frac{n}{2}-1}a_{2i}* x^i ]

    [F_1(x)=sum_{i=0}^{frac{n}{2}-1}a_{2i+1}* x^i ]

    [F(x)=F_0(x)+x* F_1(x) ]

    [F(ω^k_n)=F_0(ω^{2* k}n)+ω^k_n* F_1(ω^{2* k} n)=F_0(ω^k_{frac{n}{2}})+ω^k_n* F_1(ω^k_{frac{n}{2}}) ]

    [F(ω^{k+frac{n}{2}}n)=(-ω^k_n)=F_0(ω^{2* k}n)-ω^k_n* F_1(ω^{2* k} n)=F_0(ω^k{frac{n}{2}})-ω^k_n* F_1(ω^k_{frac{n}{2}}) ]

    分治!!!
    IDFT

    so as we make our $ ω^k_n $ to $ -ω^k_n $,Ans to a,we can solve it.
    Better

    first we consider a (n) is the smallest number that n>len1+len2 and lowerbit(n)=n 我们发现分治到边界之后的下标是原下标的二进制位翻转 so

    f[0]=0;
    for(int i=1;i<=n;i++)
    if(i&1) f[i]=f[i-1]|(n>>1);
    else f[i]=f[i>>1]>>1;

    as we calculate FFT we can use the following steps

    $ tmp=ω^k_n* a_{frac{n}{2}+k} $

    $ a_{frac{n}{2}+k}=a_{k}-tmp ( ) a_{k}=a_{k}+tmp $

    //zzd's code not writer's
    #include <bits/stdc++.h>
    using namespace std;
    const int N=1<<20;
    const double PI=acos(-1.0);
    struct C{
        double r,i;
        C(){}
        C(double a,double b){r=a,i=b;}
        C operator + (C x){return C(r+x.r,i+x.i);}
        C operator - (C x){return C(r-x.r,i-x.i);}
        C operator * (C x){return C(rx.r-ix.i,rx.i+ix.r);}
    }a[N],b[N],w[N];
    int A,B,n,L,R[N];
    void FFT(C a[],int n){
        for (int i=0;i<n;i++)
    	if (R[i]>i)
    	    swap(a[R[i]],a[i]);
        for (int t=n>>1,d=1;d<n;d<<=1,t>>=1)
    	for (int i=0;i<n;i+=(d<<1))
    	    for (int j=0;j<d;j++){
    		C tmp=w[tj]a[i+j+d];
    		a[i+j+d]=a[i+j]-tmp;
    		a[i+j]=a[i+j]+tmp;
    	    }
    }
    int main(){
        scanf("%d",&A);A++;
        scanf("%d",&B);B++;
        for (int i=0;i<A;i++)
    	scanf("%lf",&a[i].r);
        for (int i=0;i<B;i++)
    	scanf("%lf",&b[i].r);
        for (n=1,L=0;n<=A+B;n<<=1,L++);
        for (int i=0;i<n;i++){
    	R[i]=(R[i>>1]>>1)|((i&1)<<(L-1));
    	w[i]=C(cos(2.0iPI/n),sin(2.0iPI/n));
        }
        FFT(a,n),FFT(b,n);
        for (int i=0;i<n;i++)
    	a[i]=a[i]*b[i],w[i].i=-w[i].i;
        FFT(a,n);
        A--,B--;
        for (int i=0;i<=A+B;i++)
    	printf("%d ",int(a[i].r/n+0.5));
        return 0;
    }
    
    //writer's
    #include<cstdio>
    #include<cmath>
    #include<algorithm>
    #include<iostream>
    using namespace std;
    const int N=1<<20;
    const double pai=acos(-1.0);
    int pos[N],n,m,li,L;
    struct C{
        double x,y;
        C(){}
        C(double a,double b){x=a,y=b;}
        C operator + (C xx){ return C(x+xx.x,y+xx.y);}
        C operator - (C xx){ return C(x-xx.x,y-xx.y);}
        C operator * (C xx){ return C(xxx.x-yxx.y,xxx.y+yxx.x);}
    }a[N],b[N],c[N],w[N];
    void FFT(C a[],int n){
        for(int i=0;i<n;i++)
    	if(pos[i]>i)
    	    swap(a[pos[i]],a[i]);
        for(int t=n>>1,d=1;d<n;d<<=1,t>>=1)
    	for(int i=0;i<n;i+=(d<<1))
    	    for(int j=0;j<d;j++){
    		C tmp=w[tj]a[i+j+d];
    		a[i+j+d]=a[i+j]-tmp;
    		a[i+j]=a[i+j]+tmp;
    	    }
    }
    int main(){
        scanf("%d%d",&n,&m);
        n++;m++;
        for(int i=0;i<n;i++) scanf("%lf",&a[i].x);
        for(int i=0;i<m;i++) scanf("%lf",&b[i].x);
        for(li=1,L=0;li<=n+m;li<<=1,L++);
        for(int i=0;i<li;i++){
    	if(i&1) pos[i]=pos[i-1]|(li>>1);
    	else pos[i]=pos[i>>1]>>1;
    	w[i]=C(cos(2.0ipai/li),sin(2.0ipai/li));
        }
        //for(int i=0;i<li;i++) printf("%d ",pos[i]);
        FFT(a,li);
        //for(int i=0;i<li;i++) printf("%lf %lf
    ",a[i].x,a[i].y);
        FFT(b,li);
        for(int i=0;i<li;i++){
    	c[i]=a[i]*b[i];
    	w[i].y=-w[i].y;
        }
        FFT(c,li);
        for(int i=0;i<=n+m-2;i++)
    	printf("%d ",int(c[i].x/li+0.5));
        return 0;
    }
    
  • 相关阅读:
    mysql数据向Redis快速导入
    jquery.cookie.js使用
    怎么才能在职场中如鱼得水(转)
    内部类(编程思想)
    main方法原来只要放在public static类中就能跑,涨知识了
    匿名内部类--工厂
    Java通过继承外部类来建立该外部类的protected内部类的实例(转)
    监听器的使用例子 ServletContextListener
    Class.getResource()方法的使用
    maven打包资源文件(转)
  • 原文地址:https://www.cnblogs.com/ac_forever/p/10390099.html
Copyright © 2011-2022 走看看