我正在寻找一个简单的快速algorithmDCT和IDCTmatrix

我正在寻找一种简单的algorithm来执行任意大小的matrix[NxM]的快速DCT （types2）以及用于逆变换IDCT （也称为DCTtypes3）的algorithm。

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

PHP代码是可取的，但任何足够清晰的algorithm都可以。

当前matrix大小大于[200×200]的我用于实现DCT / IDCT的PHP脚本非常慢。

我正在跳跃寻找一种方法来在不到20秒的时间内完成[4000×4000]的DCT。 有谁知道如何做到这一点？

Answer 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是目的地vector[n]
• src是源vector[n]
• tmp是临时向量[2n]

这些数组不应该重叠！ 它是从我的转换类，所以我希望不要忘了复制一些东西。

• XXXrr表示目的地是真实的，源也是真实的域
• XXXrc表示目的地是真实的，源是复杂的域
• XXXcr意味着目的地是复杂的，源是真实的域名

所有的数据都是double数组，对于复数域第一个数是Real和第二个虚数部分，所以数组是2N大小。 这两个函数都使用FFT和iFFT，如果你需要代码也为他们评论我。 只是为了确定我在下面添加了不快的实现。 复制它比较容易，因为快速的使用太多的变换类层次结构

慢DFT，用于testing的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; } } //--------------------------------------------------------------------------- 

因此，为了testing，只需将代码重写为DFFTcr和iDFFTrc （或者用它们来比较您的FFT,iFFT ），然后执行您自己的FFT，iFFT有关更多信息，请参阅：

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

2D DFCT

1. 将srcmatrix调整为2幂

   通过添加零，使用快速algorithm的大小必须始终为2 ！

2. 分配NxN实matrixtmp,dst和1 NxN复vectort

3. 由DFCTrr转换行

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

4. 转置tmpmatrix

5. 由DFCTrr转换行

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

6. 转置dstmatrix

7. 通过将matrix乘以0.0625归一化dst

2D iDFCT

与上面相同，但使用iDFCTrr ，而不是乘以16.0 。

[笔记]

确保在实施自己的FFT和iFFT之前，他们给我的结果是一样的，否则DCT / iDCT将无法正常工作！