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;
    }
    
  • 相关阅读:
    HTTP 错误 404.2
    SQL Server 2008 R2如何开启数据库的远程连接(转)
    CSS中font-family:中文字体对应的英文名称
    15/18位身份证号码正则表达式(详细版)
    C#获取系统时间及时间格式
    C#正则表达式判断输入日期格式是否正确
    Linux查看机器负载
    模拟HTTP请求超时时间设置
    MySQL show命令的用法
    innodb事务隔离级别
  • 原文地址:https://www.cnblogs.com/ac_forever/p/10390099.html
Copyright © 2011-2022 走看看