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

algorithm, dct, mathematical-optimization, performance

Solution

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

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

Problem

Original source

Related problems

• How to compute Discrete Fourier Transform?