I am looking for a simple algorithm for fast DCT and IDCT of matrix [NxM]

Question

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?

Answer 1

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 domain
• XXXrc means destination is real and source is complex domain
• XXXcr means destination is complex and source is real domain

All 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:

• How to compute Discrete Fourier Transform?

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

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 !!!