我正在寻找一种用于矩阵 [NxM] 的快速 DCT 和 IDCT 的简单算法 [英] I am looking for a simple algorithm for fast DCT and IDCT of matrix [NxM]

问题描述

我正在寻找一种简单的算法来执行快速 DCT(类型 2) 的任意大小 [NxM] 矩阵，以及逆变换算法 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 代码更可取，但任何足够清晰的算法都可以.

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

当矩阵大小超过 [200x200] 时，我当前用于实现 DCT/IDCT 的 PHP 脚本非常慢.

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

我一直在寻找一种方法，可以在不到 20 秒的时间内完成高达 [4000x4000] 的 DCT.有人知道怎么做吗?

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?

推荐答案

这里是我用相同长度的 FFT 计算 1D FDCT 和 IFDCT:

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]
• dst is destination vector [n]
• src is source vector [n]
• tmp is temp vector [2n]

这些数组不应该重叠！！！它取自我的变换类所以我希望不要忘记复制一些东西.

These arrays should not overlap !!! It is taken from mine transform class so I hope did not forget to copy something.

• XXXrr 表示目的地是真实的，源也是真实的域
• XXXrc 表示目的地是真实的，源是复杂的域
• XXXcr 表示目的地是复杂的，源是真实的域
• 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

所有数据都是double数组，对于复数域，第一个数字是实部，第二个是虚部，所以数组是2N大小.这两个函数都使用 FFT 和 iFFT 如果您还需要它们的代码，请评论我.为了确保我在下面添加了它们的快速实现.复制它要容易得多，因为快速的使用了太多的转换类层次结构

All data are double arrays, for complex 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 added 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 实现:

//---------------------------------------------------------------------------
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 进行比较)工作正常然后实现你自己的 FFT,iFFT 有关更多信息，请参阅:

So for testing just rewrite the names to DFFTcr and iDFFTrc (or use them to compare to your FFT,iFFT) when the code works properly then implement your own FFT,iFFT For more info on that see:

• 如何计算离散傅立叶变换?

二维 DFCT

1. 将 src 矩阵调整为 2

1. resize src matrix to power of 2

通过添加零，要使用快速算法，大小必须始终是 2 的幂！！！

by adding zeros, to use fast algorithm the size must be always power of 2 !!!

分配NxN实矩阵tmp,dst和1xN复向量t

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

通过DFCTrr

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

• 转置tmp矩阵

• transpose tmp matrix

  通过DFCTrr

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

• 转置dst矩阵

• transpose dst matrix

  通过将矩阵乘以0.0625

  二维 iDFCT

  与上面相同，但使用 iDFCTrr 并乘以 16.0 代替.

  Is the same as above but use iDFCTrr and multiply by 16.0 instead.

  [注释]

  在实施您自己的 FFT 和 iFFT 之前，请确保它们给出与我的相同的结果，否则 DCT/iDCT 将无法正常工作！！！

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