I am looking for a simple algorithm to perform fast DCT (type 2) of a matrix of any size [NxM], and also an algorithm for the inverse transformation IDCT (also called DCT type 3).
I need a DCT-2D algorithm, but even a DCT-1D algorithm is good enough because I can use DCT-1D to implement DCT-2D (and IDCT-1D to implement IDCT-2D ).
PHP code is preferable, but any algorithm that is clear enough will do.
My current PHP script for implementing DCT/IDCT is very slow whenever matrix size is more than [200x200].
I was hopping to find a way to preform DCT of up to [4000x4000] within less than 20 seconds. Does anyone know how to do it?
Here is mine computation of 1D FDCT and IFDCT by FFT with the same length:
//---------------------------------------------------------------------------
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
is destination vector [n]
src
is source vector [n]
tmp
is temp vector [2n]
These arrays should not overlap !!! It is taken from mine transform class so I hope did not forget to copy something.
XXXrr
means destination is real and source is also real domainXXXrc
means destination is real and source is complex domainXXXcr
means destination is complex and source is real domainAll data are double
arrays, for complex domain first number is Real and second Imaginary part so array is 2N
size. Both functions use FFT and iFFT if you need code also for them comment me. Just to be sure I added not fast implementation of them below. It is much easier to copy that because fast ones use too much of the transform class hierarchy
slow DFT,iDFT implementations for testing:
//---------------------------------------------------------------------------
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;
}
}
//---------------------------------------------------------------------------
So for testing just rewrite the names to DFFTcr
and iDFFTrc
(or use them to compare to your FFT,iFFT
) when the code works properly then implement your own FFT,iFFT For more info on that see:
2D DFCT
resize src
matrix to power of 2
by adding zeros, to use fast algorithm the size must be always power of 2
!!!
allocate NxN
real matrices tmp,dst
and 1xN
complex vector t
transform lines by DFCTrr
DFCT(tmp.line(i),src.line(i),t,N)
transpose tmp
matrix
transform lines by DFCTrr
DFCT(dst.line(i),tmp.line(i),t,N)
transpose dst
matrix
normalize dst
by multiply matrix by 0.0625
2D iDFCT
Is the same as above but use iDFCTrr
and multiply by 16.0
instead.
[Notes]
Be sure before implementing your own FFT and iFFT that they give the same result as mine otherwise the DCT/iDCT will not work properly !!!
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With