我在找一个简单的算法，快速DCT和矩阵的IDCT [N×M的] [英] I am looking for a simple algorithm for fast DCT and IDCT of matrix [NxM]

问题描述

我在找一个简单的算法来进行快速 DCT （类型2 ）任何尺寸[N×M个]的矩阵，并且也是一个算法的逆变换 IDCT <的/一>（也称为DCT类型3）。

I am looking for a simple algorithm to perform fast DCT (type 2) of a matrix of any size [NxM], and also an algorithm for the inverse transformation IDCT (also called DCT type 3).

我需要一个DCT-2D算法，但即使是DCT-1D算法是不够好，因为我可以使用DCT-1D实现DCT-2D（和IDCT-1D实现IDCT-2D）。

I need a DCT-2D algorithm, but even a DCT-1D algorithm is good enough because I can use DCT-1D to implement DCT-2D (and IDCT-1D to implement IDCT-2D ).

PHP code是preferable，但任何算法显然不够就行了。

PHP code is preferable, but any algorithm that is clear enough will do.

实现DCT / IDCT我目前的PHP脚本是很慢的，每当矩阵大小超过[200×200]。

My current PHP script for implementing DCT/IDCT is very slow whenever matrix size is more than [200x200].

我跳槽找到一种方式来达到[4000x4000]的preform DCT不到20秒之内。有谁知道该怎么办呢？

I was hopping to find a way to preform DCT of up to [4000x4000] within less than 20 seconds. Does anyone know how to do it?

推荐答案

下面是一维FDCT和IFDCT矿计算通过FFT长度相同：

Here is mine computation of 1D FDCT and IFDCT by FFT with the same length:

//---------------------------------------------------------------------------
void DFCTrr(double *dst,double *src,double *tmp,int n)
    {
    // exact normalized DCT II by N DFFT
    int i,j;
    double nn=n,a,da=(M_PI*(nn-0.5))/nn,a0,b0,a1,b1,m;
    for (j=  0,i=n-1;i>=0;i-=2,j++) dst[j]=src[i];
    for (j=n-1,i=n-2;i>=0;i-=2,j--) dst[j]=src[i];
    DFFTcr(tmp,dst,n);
    m=2.0*sqrt(2.0);
    for (a=0.0,j=0,i=0;i<n;i++,j+=2,a+=da)
        {
        a0=tmp[j+0]; a1= cos(a);
        b0=tmp[j+1]; b1=-sin(a);
        a0=(a0*a1)-(b0*b1);
        if (i) a0*=m; else a0*=2.0;
        dst[i]=a0;
        }
    }
//---------------------------------------------------------------------------
void iDFCTrr(double *dst,double *src,double *tmp,int n)
    {
    // exact normalized DCT III = iDCT II by N iDFFT
    int i,j;
    double nn=n,a,da=(M_PI*(nn-0.5))/nn,a0,m,aa,bb;
    m=1.0/sqrt(2.0);
    for (a=0.0,j=0,i=0;i<n;i++,j+=2,a+=da)
        {
        a0=src[i];
        if (i) a0*=m;
        aa= cos(a)*a0;
        bb=+sin(a)*a0;
        tmp[j+0]=aa;
        tmp[j+1]=bb;
        }
    m=src[0]*0.25;
    iDFFTrc(src,tmp,n);
    for (j=  0,i=n-1;i>=0;i-=2,j++) dst[i]=src[j]-m;
    for (j=n-1,i=n-2;i>=0;i-=2,j--) dst[i]=src[j]-m;
    }
//---------------------------------------------------------------------------

• dst为目标向量[N]
• 在src为源矢量[N]
• 在TMP是临时矢量[2N]
• 在这些不应重叠!!!
• 在它转换类，所以我希望也没有忘记复制一些取自煤矿
• XXXrr意味着目的地是真实的，来源也是真正的域
• XXXrc意味着目的地是真实的，来源是复杂的域
• XXXcr表示目标是复杂的，来源是真实的域
• 所有的数据都是双阵列
• 的COMPEX域第一个数字是真实的，第二虚部分，因此数组为2N尺寸
• 在这两个函数使用FFT和IFFT，如果你需要code也为他们评我
• 只是要确定我加入他们的不是快速实现以下
• 这是很容易复制，因为快的人使用太多的改变类层次结构的
• dst is destination vector [n]
• src is source vector [n]
• tmp is temp vector [2n]
• these should not overlap !!!
• it is taken from mine transform class so I hope did not forget to copy something
• XXXrr means destination is real and source is also real domain
• XXXrc means destination is real and source is complex domain
• XXXcr means destination is complex and source is real domain
• all data are double arrays
• for compex domain first number is Real and second Imaginary part so array is 2N size
• both functions use FFT and iFFT if you need code also for them comment me
• just to be sure i add not fast implementation of them below
• it is much easier to copy that because fast ones use too much of the transform class hierarchy

慢DFT和IDFT的实现来进行测试：

slow DFT,iDFT implementations for testing:

//---------------------------------------------------------------------------
void transform::DFTcr(double *dst,double *src,int n)
    {
    int i,j;
    double a,b,a0,_n,q,qq,dq;
    dq=+2.0*M_PI/double(n); _n=2.0/double(n);
    for (q=0.0,j=0;j<n;j++,q+=dq)
        {
        a=0.0; b=0.0;
        for (qq=0.0,i=0;i<n;i++,qq+=q)
            {
            a0=src[i];
            a+=a0*cos(qq);
            b+=a0*sin(qq);
            }
        dst[j+j  ]=a*_n;
        dst[j+j+1]=b*_n;
        }
    }
//---------------------------------------------------------------------------
void transform::iDFTrc(double *dst,double *src,int n)
    {
    int i,j;
    double a,a0,a1,b0,b1,q,qq,dq;
    dq=+2.0*M_PI/double(n);
    for (q=0.0,j=0;j<n;j++,q+=dq)
        {
        a=0.0;
        for (qq=0.0,i=0;i<n;i++,qq+=q)
            {
            a0=src[i+i  ]; a1=+cos(qq);
            b0=src[i+i+1]; b1=-sin(qq);
            a+=(a0*a1)-(b0*b1);
            }
        dst[j]=a*0.5;
        }
    }
//---------------------------------------------------------------------------

• 在这样的测试只改写名DFFTcr和iDFFTrc（或用它们来比较你FFT，IFFT）
• 当code正常工作，然后实现自己的FFT，IFFT

2D DFCT

1.resize SRC矩阵2的功率

1.resize src matrix to power of 2

• 加入零
• 要使用快速算法的大小必须始终功率2！

2.allocate的N×N实矩阵TMP，DST和的1×N复杂的向量t

2.allocate NxN real matrices tmp,dst and 1xN complex vector t

3.transform线由DFCTrr

3.transform lines by DFCTrr

DFCT(tmp.line(i),src.line(i),t,N)

4.transpose TMP矩阵

4.transpose tmp matrix

5.transform线由DFCTrr

5.transform lines by DFCTrr

DFCT(dst.line(i),tmp.line(i),t,N)

6.transpose DST矩阵

6.transpose dst matrix

7.normalize DST由乘矩阵0.0625;

7.normalize dst by multiply matrix by 0.0625;

2D iDFCT

• 是与上述相同，但使用iDFCTrr和16.0乘

[注意事项]

• 请务必实现自己的FFT和IFFT，他们给予同样的结果之前，为我的
• ，否则的DCT / IDCT将无法正常工作！
• be sure before implementing your own FFT and iFFT that they give the same result as mine
• otherwise the DCT/iDCT will not work properly !!!

这篇关于我在找一个简单的算法，快速DCT和矩阵的IDCT [N×M的]的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！