java里的fft2和fft:)

fft是改网上c代码的,fft2是自己写的,谁能不能给一个求任意长度(即不需要将输入补成2的指数倍)的fft程序

   public class FftList{    //提供fft的输出
        double[] fr;
        double[] fi;
    }

    public class Fft2List{    //提供fft2的输出
        double[][] outr;
        double[][] outi;
    }

    //inr和ini为输入图像的实部和虚部,没有虚部只要令ini=null就可,w和h为输
//入图像的宽和高,ncols和nrows为输出图像的宽和高,注意ncols和nrows一定要是2的
//指数,ifft=true为逆fft2运算
    //此函数相当于matlab下的fft2(X,nrows,ncols)
    public Fft2List fft2(double[][] inr,double[][] ini,int w,int h,int n
cols,int nrows,boolean ifft)
    {

        int i,j,k;
        double[] pr=new double[ncols];      //暂存每行实部
        double[] pi=new double[ncols];      //暂存每行虚部
       k=(int)(Math.log(ncols)/Math.log(2.0));               //k=log2(ncols)

        Fft2List fft2list = new Fft2List();
        fft2list.outr = new double[ncols][nrows];
        fft2list.outi = new double[ncols][nrows];

        ////////////先对每行求fft/////////////
         for(i=0;i<h;i++)
        {
            for(j=0;j

    //inr和ini为输入图像的实部和虚部,没有虚部只要令ini=null就可,w和h为输
//入图像的宽和高,ncols和nrows为输出图像的宽和高,注意ncols和nrows一定要是2的
//指数,ifft=true为逆fft2运算
    //此函数相当于matlab下的fft2(X,nrows,ncols)
    public Fft2List fft2(double[][] inr,double[][] ini,int w,int h,int n
cols,int nrows,boolean ifft)
    {

        int i,j,k;
        double[] pr=new double[ncols];      //暂存每行实部
        double[] pi=new double[ncols];      //暂存每行虚部
       k=(int)(Math.log(ncols)/Math.log(2.0));               //k=log2(ncols)

        Fft2List fft2list = new Fft2List();
        fft2list.outr = new double[ncols][nrows];
        fft2list.outi = new double[ncols][nrows];

        ////////////先对每行求fft/////////////
         for(i=0;i<h;i++)
        {
            for(j=0;j<ncols;j++)               //清0
            {
                pr[j]=0;
                pi[j]=0;
            }
            for(j=0;j<w;j++)            //将每行实部和虚部赋给pr和pi
            {
                pr[j]=inr[j][i];
                if(ini != null)           //输入可能没有虚部
                    pi[j]=ini[j][i];
            }
            FftList fftlist = this.fft(pr,pi,ncols,k,ifft);

            for(j=0;j<ncols;j++)      //将每行结果保存到outr和outi中去
            {
                fft2list.outr[j][i]=fftlist.fr[j];
                fft2list.outi[j][i]=fftlist.fi[j];
            }
        }

        /////////////////////////////////////////

         pr=new double[nrows];
         pi=new double[nr
                pr[j]=0;
                pi[j]=0;
            }
            for(j=0;j<w;j++)            //将每行实部和虚部赋给pr和pi
            {
                pr[j]=inr[j][i];
                if(ini != null)           //输入可能没有虚部
                    pi[j]=ini[j][i];
            }
            FftList fftlist = this.fft(pr,pi,ncols,k,ifft);

            for(j=0;j<ncols;j++)      //将每行结果保存到outr和outi中去
            {
                fft2list.outr[j][i]=fftlist.fr[j];
                fft2list.outi[j][i]=fftlist.fi[j];
            }
        }

        /////////////////////////////////////////

         pr=new double[nrows];
         pi=new double[nrows];
         k=(int)(Math.log(nrows)/Math.log(2.0));               //k=log2(nrows)

        /////////////再对每列求fft/////////////////

        for(i=0;i<ncols;i++)
        {
            for(j=0;j<nrows;j++)          //清0
            {
                pr[j]=0;
                pi[j]=0;
            }
            for(j=0;j<h;j++)
            {
                pr[j]=fft2list.outr[i][j];        //将上一部分结果的每一
                                                  // 列赋给pr和pi
                pi[j]=fft2list.outi[i][j];
            }
            FftList fftlist = this.fft(pr,pi,nrows,k,ifft);
            for(j=0;j<nrows;j++)
            {
                fft2list.outr[i][j]=fftlist.fr[j];
                fft2list.outi[i][j]=fftlist.fi[j];
            }nt)(Math.log(nrows)/Math.log(2.0));               //k=log2(nrows)

        /////////////再对每列求fft/////////////////

        for(i=0;i<ncols;i++)
        {
            for(j=0;j<nrows;j++)          //清0
            {
                pr[j]=0;
                pi[j]=0;
            }
            for(j=0;j<h;j++)
            {
                pr[j]=fft2list.outr[i][j];        //将上一部分结果的每一
                                                  // 列赋给pr和pi
                pi[j]=fft2list.outi[i][j];
            }
            FftList fftlist = this.fft(pr,pi,nrows,k,ifft);
            for(j=0;j<nrows;j++)
            {
                fft2list.outr[i][j]=fftlist.fr[j];
                fft2list.outi[i][j]=fftlist.fi[j];
            }
        }
        return fft2list;
    }


     //n=2^k,其他的计算会出错
    public FftList fft(double[] pr,double[] pi,int n,int k,boolean ifft)
    {
        int it,m,is,i,j,nv,l0;
        double p,q,s,vr,vi,poddr,poddi;
        FftList cl = new FftList();
        cl.fr = new double[n];
        cl.fi = new double[n];

        for (it=0; it<=n-1; it++)
        {
            m=it; is=0;
            for (i=0; i<=k-1; i++)
            { j=m/2; is=2*is+(m-2*j); m=j;}
            cl.fr[it]=pr[is]; cl.fi[it]=pi[is];
        }

        pr[0]=1.0; pi[0
        return fft2list;
    }


     //n=2^k,其他的计算会出错
    public FftList fft(double[] pr,double[] pi,int n,int k,boolean ifft)
    {
        int it,m,is,i,j,nv,l0;
        double p,q,s,vr,vi,poddr,poddi;
        FftList cl = new FftList();
        cl.fr = new double[n];
        cl.fi = new double[n];

        for (it=0; it<=n-1; it++)
        {
            m=it; is=0;
            for (i=0; i<=k-1; i++)
            { j=m/2; is=2*is+(m-2*j); m=j;}
            cl.fr[it]=pr[is]; cl.fi[it]=pi[is];
        }

        pr[0]=1.0; pi[0]=0.0;
        p=2*Math.PI/(1.0*n);

        pr[1]=Math.cos(p); pi[1]=-Math.sin(p);

        if (ifft) pi[1]=-pi[1];    //求逆fft

        for (i=2; i<=n-1; i++)
        {
            p=pr[i-1]*pr[1]; q=pi[i-1]*pi[1];
            s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);
            pr[i]=p-q; pi[i]=s-p-q;
        }

        for (it=0; it<=n-2; it=it+2)
        {
            vr=cl.fr[it]; vi=cl.fi[it];
            cl.fr[it]=vr+cl.fr[it+1]; cl.fi[it]=vi+cl.fi[it+1];
            cl.fr[it+1]=vr-cl.fr[it+1]; cl.fi[it+1]=vi-cl.fi[it+1];
        }

        m=n/2; nv=2;
        for (l0=k-2; l0>=0; l0--)
        {
            m=m/2; nv=2*nv;
        pr[1]=Math.cos(p); pi[1]=-Math.sin(p);

        if (ifft) pi[1]=-pi[1];    //求逆fft

        for (i=2; i<=n-1; i++)
        {
            p=pr[i-1]*pr[1]; q=pi[i-1]*pi[1];
            s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);
            pr[i]=p-q; pi[i]=s-p-q;
        }

        for (it=0; it<=n-2; it=it+2)
        {
            vr=cl.fr[it]; vi=cl.fi[it];
            cl.fr[it]=vr+cl.fr[it+1]; cl.fi[it]=vi+cl.fi[it+1];
            cl.fr[it+1]=vr-cl.fr[it+1]; cl.fi[it+1]=vi-cl.fi[it+1];
        }

        m=n/2; nv=2;
        for (l0=k-2; l0>=0; l0--)
        {
            m=m/2; nv=2*nv;
            for (it=0; it<=(m-1)*nv; it=it+nv)
                for (j=0; j<=(nv/2)-1; j++)
                {
                    p=pr[m*j]*cl.fr[it+j+nv/2];
                    q=pi[m*j]*cl.fi[it+j+nv/2];
                    s=pr[m*j]+pi[m*j];
                    s=s*(cl.fr[it+j+nv/2]+cl.fi[it+j+nv/2]);
                    poddr=p-q; poddi=s-p-q;
                    cl.fr[it+j+nv/2]=cl.fr[it+j]-poddr;
                    cl.fi[it+j+nv/2]=cl.fi[it+j]-poddi;
                    cl.fr[it+j]=cl.fr[it+j]+poddr;
                    cl.fi[it+j]=cl.fi[it+j]+poddi;
                }
        }

        if (ifft)
            for (i=0;

        m=n/2; nv=2;
        for (l0=k-2; l0>=0; l0--)
        {
            m=m/2; nv=2*nv;
            for (it=0; it<=(m-1)*nv; it=it+nv)
                for (j=0; j<=(nv/2)-1; j++)
                {
                    p=pr[m*j]*cl.fr[it+j+nv/2];
                    q=pi[m*j]*cl.fi[it+j+nv/2];
                    s=pr[m*j]+pi[m*j];
                    s=s*(cl.fr[it+j+nv/2]+cl.fi[it+j+nv/2]);
                    poddr=p-q; poddi=s-p-q;
                    cl.fr[it+j+nv/2]=cl.fr[it+j]-poddr;
                    cl.fi[it+j+nv/2]=cl.fi[it+j]-poddi;
                    cl.fr[it+j]=cl.fr[it+j]+poddr;
                    cl.fi[it+j]=cl.fi[it+j]+poddi;
                }
        }

        if (ifft)
            for (i=0; i<=n-1; i++)
            {
                cl.fr[i]=cl.fr[i]/(1.0*n);
                cl.fi[i]=cl.fi[i]/(1.0*n);
            }

        return cl;
    }

阅读更多
个人分类: java
上一篇vb.net下对注册表的各种操作!(谢neu2000指路:))
想对作者说点什么? 我来说一句

没有更多推荐了,返回首页

关闭
关闭
关闭