一尘不染

我正在寻找用于矩阵[NxM]的快速DCT和IDCT的简单算法

algorithm

我正在寻找一种简单的算法来执行任意大小[NxM]的矩阵的快速DCT(类型2),以及一种用于逆变换IDCT(也称为DCT类型3)的算法。

我需要DCT-2D算法,但是即使DCT-1D算法也足够好,因为我可以使用DCT-1D来实现DCT-2D(和IDCT-1D来实现IDCT-2D)。

最好使用PHP代码,但是任何足够清晰的算法都可以。

每当矩阵大小大于[200x200]时,我当前用于实现DCT / IDCT的PHP脚本就会非常慢。

我希望找到一种方法,在不到20秒的时间内完成高达[4000x4000]的DCT。有人知道怎么做吗?


阅读 213

收藏
2020-07-28

共1个答案

一尘不染

这是使用相同长度的FFT对一维FDCT和IFDCT进行的矿山计算:

//---------------------------------------------------------------------------
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 意味着目的地是复杂的,而来源是真实的领域

所有数据均为double数组,对于复杂域,第一个数字为实数,第二个虚部为数组2N。如果您还需要代码,则这两个函数都使用 FFT
iFFT 。只是为了确保我在下面没有添加快速实现它们。复制起来容易得多,因为快速的代码占用了太多的转换类层次结构

用于测试的慢速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;
        }
    }
//---------------------------------------------------------------------------

因此,对于测试,只需在代码正常运行时将名称重写为DFFTcriDFFTrc(或使用它们与进行比较FFT,iFFT),然后实现自己的

二维DFCT

  1. 调整src矩阵的幂2

通过添加零,要使用快速算法,大小必须始终为幂2

  1. 分配NxN实矩阵tmp,dst1xN复矢量t

  2. 通过变换线DFCTrr

     DFCT(tmp.line(i),src.line(i),t,N)
    
  3. 转置tmp矩阵

  4. 通过变换线DFCTrr

     DFCT(dst.line(i),tmp.line(i),t,N)
    
  5. 转置dst矩阵

  6. dst 通过乘以矩阵归一化0.0625

二维iDFCT

与上述相同,但改用iDFCTrr乘以16.0

[笔记]

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

2020-07-28