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

Question

我正在寻找一种简单的算法来执行任意大小 [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。有谁知道怎么做？

Answer 1

这是我对相同长度的 1D FDCT 和 IFDCT 通过 FFT 的计算：

//---------------------------------------------------------------------------
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;
        }
    }
//---------------------------------------------------------------------------

因此，为了进行测试，只需在代码正常工作时将名称重写为DFFTcr和iDFFTrc（或使用它们与您的 进行比较），然后实现您自己的FFT，iFFT有关详细信息，请参阅：FFT,iFFT

• 如何计算离散傅里叶变换？

二维离散余弦变换

1. 将矩阵大小调整src为 的幂2

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

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

3. 将线变换为DFCTrr

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

4. 转置tmp矩阵

5. 将线变换为DFCTrr

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

6. 转置dst矩阵

7. dst通过将矩阵乘以标准化0.0625

二维离散离散傅里叶变换

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

[笔记]

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