logo标签列表

寻找一个简单的算法，用于矩阵[NxM]的快速离散余弦变换（DCT）和逆离散余弦变换（IDCT）。

performancealgorithmmathematical-optimizationdct
6
我正在寻找一个简单的算法，可以快速执行任意大小[NxM]矩阵的DCT（类型2），以及逆变换IDCT（也称为DCT类型3）的算法。需要一个DCT-2D算法，但是即使是DCT-1D算法也足够好，因为我可以使用DCT-1D来实现DCT-2D（并使用IDCT-1D来实现IDCT-2D）。PHP代码最好，但任何清晰的算法都可以。我的当前PHP脚本在矩阵大小超过[200x200]时非常缓慢。我希望找到一种在不到20秒的时间内执行高达[4000x4000]的DCT的方法。有人知道如何做吗？
1
如果代码不是太多，你能否给我们展示一下你目前已经完成的部分呢？ - Teepeemm
到目前为止，我或多或少是按照它们的正式定义计算DCT和IDCT，这非常慢。 - algo-rithm
rwong，谢谢，我会看一下的。你知道将DFT转换为DCT或者反过来的快速/简单方法吗？ - algo-rithm
1
你想用哪个DCT（1、2、3、4）？我建议使用FFT计算DCT，因为快速DCT非常复杂，往往比通过FFT计算更慢。你还想要直接的还是归一化的？此外，DFCT算法的文档非常差！（从未见过功能实现，只有不完整的方程式） - Spektre
1
你可以在网络上找到许多快速余弦变换，但它们都是针对像8x8这样的固定大小，完全没有帮助...要实现直接离散余弦变换(DFCT)或反向离散余弦变换(iDFCT)，你还需要实现正弦变换(DST)、反向正弦变换(iDST)，而且递归与DFT和iDFT中那样容易分离，你将得到三个项，而不是两个... - Spektre
显示剩余2条评论
1个回答
2

这里是我使用相同长度的FFT计算1D 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;
        }
    }
//---------------------------------------------------------------------------

因此，为了测试，请将名称重写为DFFTcr和iDFFTrc（或使用它们与您的FFT，iFFT进行比较）。当代码正常工作时，再实现自己的FFT，iFFT。有关更多信息，请参见：

2D DFCT

  1. resize src matrix to power of 2

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

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

  3. transform lines by DFCTrr

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

  4. transpose tmp matrix

  5. transform lines by DFCTrr

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

  6. transpose dst matrix

  7. normalize dst by multiply matrix by 0.0625

2D iDFCT

2D iDFCT与上述相同，但使用iDFCTrr并乘以16.0。

[注]

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

Spektre，非常感谢，但我真的需要[NxM]。你知道我怎么修改你的代码[NxN]来得到[NxM]吗？ - algo-rithm
应用任何类型的快速算法，N和M都必须是2的幂次方。您可以像在NxN上一样对NxM进行2D变换，只需调整转置函数即可...(唯一可能改变的是归一化常数...大多数算法使用NM、1/NM，但我更喜欢16、1/16，因为它不会改变信号的幅度)。 - Spektre
我尝试编译和运行DFTcr函数，但似乎当dst数组大小超过限制时，在循环内部的某个点上dst[j+j+1]=b*_n会引发异常。 - algo-rithm
你是不是忘记将复杂数组分配为实际数组的两倍了？ - Spektre
当您的DCT / iDCT正常工作时，请不要忘记实现快速FT，iFT，以使您的DCT，iDCT也变得更快。 - Spektre
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接