回 帖 发 新 帖 刷新版面

主题:[讨论]fft  算法

void four1(double data[65],int nn,int isign)
{
    int n,j,i,m,mmax,istep ;
    double tempr,tempi,theta,wpr,wpi,wr,wi,wtemp ;
    n=2*nn ;
    j=1 ;
    for(i=1;i<=n;i=i+2)
    {
        if(j>i) {
            tempr=data[j];
            tempi=data[j+1];
            data[j]=data[i];
            data[j+1]=data[i+1];
            data[i]=tempr ;
            data[i+1]=tempi ;
        }
        m=n/2 ;
        while(m>=2&&j>m) {
            j=j-m ;  m=m/2 ;
        }
        j=j+m ;
    }
    mmax=2 ;
    while(n>mmax)
    {
        istep=2*mmax ;
        theta=6.28318530717959/(isign*mmax);
        wpr=-2.0*sin(0.5*theta)*sin(0.5*theta);
        wpi=sin(theta);
        wr=1.0 ;
        wi=0.0 ;
        for(m=1;m<=mmax;m=m+2)
        {
            for(i=m;i<=n;i=i+istep)
            {
                j=i+mmax ;
                tempr=double(wr)*data[j]-double(wi)*data[j+1];
                tempi=double(wr)*data[j+1]+double(wi)*data[j];
                data[j]=data[i]-tempr ;
                data[j+1]=data[i+1]-tempi ;
                data[i]=data[i]+tempr ;
                data[i+1]=data[i+1]+tempi ;
            }
            wtemp=wr ;
            wr=wr*wpr-wi*wpi+wr ;
            wi=wi*wpr+wtemp*wpi+wi ;
        }
        mmax=istep ;
    }
}

void twofft(double data1[],double data2[],double fft1[],
double fft2[],int&n)
{
    int j,n2,j2 ;
    double c1r,c1i,c2r,c2i,conjr,conji,h1r,h1i,h2r,h2i ;
    c1r=0.5 ;
    c1i=0.0 ;
    c2r=0.0 ;
    c2i=-0.5 ;
    for(j=1;j<=n;j++) {
        fft1[2*j-1]=data1[j];
        fft1[2*j]=data2[j];
    }
    four1(fft1,n,1);
    fft2[1]=fft1[2];
    fft2[2]=0.0 ;
    fft1[2]=0.0 ;
    n2=2*(n+2);
    for(j=2;j<=n/2+1;j++)
    {
        j2=2*j ;
        conjr=fft1[n2-j2-1];
        conji=-fft1[n2-j2];
        h1r=c1r*(fft1[j2-1]+conjr)-c1i*(fft1[j2]+conji);
        h1i=c1i*(fft1[j2-1]+conjr)+c1r*(fft1[j2]+conji);
        h2r=c2r*(fft1[j2-1]-conjr)-c2i*(fft1[j2]-conji);
        h2i=c2i*(fft1[j2-1]-conjr)+c2r*(fft1[j2]-conji);
        fft1[j2-1]=h1r ;          fft1[j2]=h1i ;
        fft1[n2-j2-1]=h1r ;     fft1[n2-j2]=-h1i ;
        fft2[j2-1]=h2r ;          fft2[j2]=h2i ;
        fft2[n2-j2-1]=h2r ;      fft2[n2-j2]=-h2i ;
    }
}
int cint(double x)
{
    int temp ;
    double iprt ;
    if(x>0) {7758521
        x=modf(x,&iprt);
        if(fabs(x)<0.5) temp=int(iprt);
        else  temp=int(iprt+1);
    }
    else if(x==0) temp=0 ;
    else  {
        x=modf(x,&iprt);
        if(fabs(x)<0.5) temp=int(iprt);
        else  temp=int(iprt)-1 ;
    }
    return temp ;
}

void prntft(double data[],double nn2)
{
    int n,m,mm ;
    cout<<"n     real(n)     imag.(n)     real(N-n)     imag.(N-n)"<<endl ;
    cout<<setw(1)<<"0" ;
    cout<<setw(14)<<data[1]; cout<<setw(14)<<data[2];
    cout<<setw(14)<<data[1]; cout<<setw(14)<<data[2]<<endl ;
    for(n=3;n<=(nn2/2)+1;n=n+2)
    {
        m=(n-1)/2 ;
        mm=nn2+2-n ;
        cout<<setiosflags(ios :: fixed);
        cout<<setprecision(0)<<setw(1)<<m ;
        cout<<setprecision(6)<<setw(14)<<data[n];
        cout<<setprecision(6)<<setw(14)<<data[n+1];
        cout<<setprecision(6)<<setw(14)<<data[mm];
        cout<<setprecision(6)<<setw(14)<<data[mm+1]<<endl ;
    }
}
void main()
{
    //program   d12r2  driver   for   routine   twofft
    int n,i,n2,isign ;
    double data1[33],data2[33],fft1[65],fft2[65],per,x ;
    n=32 ;
    n2=2*n ;
    per=8.0 ;
    const double pi=3.1415926 ;
    for(i=1;i<=n;i++) {
        x=2.0*pi*i/per ;
        data1[i]=cint(cos(x));
        data2[i]=cint(sin(x));
    }
    twofft(data1,data2,fft1,fft2,n);
    cout<<setw(1)<<"Fourier   transform   of   first   function:"<<endl ;
    prntft(fft1,n2);
    cout<<setw(1)<<"Fourier   transform   of   second   function:"<<endl ;
    prntft(fft2,n2);
    //invert   transform
    isign=-1 ;
    four1(fft1,n,isign);
    cout<<setw(1)<<"Inverted   transform   =   first   function:"<<endl ;
    prntft(fft1,n2);
    four1(fft2,n,isign);
    cout<<setw(1)<<"Inverted   transform   =   second   function:"<<endl ;
    prntft(fft2,n2);
}
/* 希望和和大家讨论 ,我对这个程序不是很明白,但是去已经是正确的 */

回复列表 (共12个回复)

沙发

不能写一些注释吗 让别人能快点看啊?!

板凳

fft?? 快速傅立叶变化么~~  呵呵 先把理论的东西弄清楚了吧

3 楼


同意,确实应该了解理论知识 
不然不好理解 碟行算法

4 楼

我只做了基2的,时间抽选和频率抽选,不用这么长程序吧

5 楼

// 按时间抽选法 带地址变换的 基2 fft
    
    void fft(complex *seq,uinteger l)
    {
        //printf("fft...\n");
        uinteger n(1<<l),i,j;
        complex temp;
        //倒位序地址变换
        
        for(i=0;i<n;++i)
        {
            uinteger iv=inorder(i,l);
            if(i<iv)
            {
                temp=seq[i];
                seq[i]=seq[iv];
                seq[iv]=temp;
            }
        }
        
        uinteger len(1),l2(2);
        //l轮蝶形运算
        for(i=1;i<=l;++i,len=l2,l2<<=1)
        {
            real q=pi/len;
            complex u((real)cos(q),(real)(-sin(q)));
            complex w(1);
            for(uinteger k=0;k<len;++k)
            {
                for(j=0;j<n;j+=l2)
                {
                    //蝶形迭代
                    uinteger s1=j+k;
                    uinteger s2=s1+len;
                    temp=seq[s2]*w;
                    seq[s2]=seq[s1]-temp;
                    seq[s1]+=temp;
                }
                w*=u;
            }
        }
    }

这个蛮简单的

6 楼

能不能 谈谈 倒位序地址变换
        
        for(i=0;i<n;++i)
        {
            uinteger iv=inorder(i,l);
            if(i<iv)
            {
                temp=seq[i];
                seq[i]=seq[iv];
                seq[iv]=temp;
            }
        }

这段程序的思路呀 ?
麻烦了 

我上面的程序 能实现两种数据输入格式的数据变换 ,而且还有反变换 
所以看起来比较长 

7 楼

void butterfly4t(struct complex*xin,int l)
{
    int i,j,r,s,a,b,m,m1,m2,m3,m4,N ;
    static double u,v,ar,rr[4],ii[4];
    N=(int)(pow(4,l));
    for(i=1;i<=l;i++)
    {
        m1=pow(4,i-1);
        m2=4*m1 ;
        m3=N/m2 ;
        for(j=0;j<m3;j++)
        {
            m4=m2*j ;
            for(r=0;r<m1;r++)
            {
                for(a=0;a<4;a++)
                {
                    s=a*m1+m4+r ;
                    rr[s]=xin[s].x ;
                    ii[s]=xin[s].y ;
                }
                
                for(b=0;b<4;b++)
                {
                    u=0 ;
                    v=0 ;
                    for(m=0;m<4;m++)
                    {
                        s=m4+r+m*m1 ;
                        ar=2*M_PI*m*(b*m1+r)/m2 ;
                        u+=rr[s]*cos(ar)+ii[s]*sin(ar);
                        v+=ii[s]*cos(ar)-rr[s]*sin(ar);
                    }
                    s=m4+r+b*m1 ;
                    xin[s].x=u ;
                    xin[s].y=v ;
                }
         }
        }
    }
}这个是 FFT4中的碟行算法 自己加个 倒位序地址变换
就是 FFT4了

8 楼

倒位序我这么变的


inline uinteger inorder(uinteger addr, uinteger l)
    {
        uinteger ioaddr(InOrder::map[addr & 0xFF]);
        if(l>8)
        {
            addr>>=8;
            while(l>8)
            {
                ioaddr=(ioaddr<<8) | InOrder::map[addr & 0xFF];
                addr>>=8;
                l-=8;
            }
        }
        return ioaddr >> (8-l);
        //printf("%u -> %u l=%u\n",addr,ioaddr,l);
        //return ioaddr;
    }


unsigned char InOrder::map[256]=
{0,128,64,192,32,160,96,224,16,144,80,208,48,176,112,240,8,136,72,200,40,
168,104,232,24,152,88,216,56,184,120,248,4,132,68,196,36,164,100,228,20,
148,84,212,52,180,116,244,12,140,76,204,44,172,108,236,28,156,92,220,60,
188,124,252,2,130,66,194,34,162,98,226,18,146,82,210,50,178,114,242,10,138,
74,202,42,170,106,234,26,154,90,218,58,186,122,250,6,134,70,198,38,166,102,
230,22,150,86,214,54,182,118,246,14,142,78,206,46,174,110,238,30,158,94,222,
62,190,126,254,1,129,65,193,33,161,97,225,17,145,81,209,49,177,113,241,9,
137,73,201,41,169,105,233,25,153,89,217,57,185,121,249,5,133,69,197,37,165,
101,229,21,149,85,213,53,181,117,245,13,141,77,205,45,173,109,237,29,157,
93,221,61,189,125,253,3,131,67,195,35,163,99,227,19,147,83,211,51,179,115,
243,11,139,75,203,43,171,107,235,27,155,91,219,59,187,123,251,7,135,71,199,
39,167,103,231,23,151,87,215,55,183,119,247,15,143,79,207,47,175,111,239,31,
159,95,223,63,191,127,255};

9 楼

基4的倒位序跟基2的一样吗,我没弄过基4的

10 楼


你 看上面的代码 就知道不一样了 

其实 基二是 二进制的反序
     基四是 四进制的反序

我来回复

您尚未登录,请登录后再回复。点此登录或注册