我正在寻找一种简单的算法来执行任何大小的矩阵[NxM]的快速
DCT(类型2),以及用于逆变换
IDCT(也称为DCT类型3)的算法.
我需要DCT-2D算法,但即使是DCT-1D算法也足够好,因为我可以使用DCT-1D来实现DCT-2D(和IDCT-1D来实现IDCT-2D).
PHP代码是首选,但任何足够清晰的算法都可以.
当矩阵大小超过[200×200]时,我当前用于实现DCT / IDCT的PHP脚本非常慢.
我正在寻找一种方法,在不到20秒的时间内完成最高[4000×4000]的DCT.有谁知道怎么做?
最佳答案 这是通过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表示目标是复杂的,源是真实域
所有数据都是双数组,对于复数域,第一个数字是实数,第二个虚数部分,因此数组是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有关详细信息,请参阅:
> How to compute Discrete Fourier Transform?
2D DFCT
>将src矩阵调整为2的幂
通过添加零,使用快速算法,大小必须始终为2!
>分配NxN实矩阵tmp,dst和1xN复矢量t
>通过DFCTrr转换行
DFCT(tmp.line(i),src.line(i),t,N)
>转置tmp矩阵
>通过DFCTrr转换行
DFCT(dst.line(i),tmp.line(i),t,N)
>转置dst矩阵
>通过乘以矩阵将dst归一化0.0625
2D iDFCT
与上面相同,但使用iDFCTrr并乘以16.0.
[笔记]
在实现自己的FFT和iFFT之前,请确保它们与我的结果相同,否则DCT / iDCT将无法正常工作!