a_dct.cpp

一个不错

/* Copyright is licensed under GNU LGPL.                 by I.J.Wang 2006

   Fast DCT transformation

   Build: make a_dct
*/
#include <math.h>
#include "../src/wy_array.h"
#include "../src/wy_uty.h"
#include <stdlib.h>

// DCT of N^2 implement
//
//  F(s)= c(s)* Zigma(m=0,M-1){ f(m)*cos(PI*s*(2*m+1)/(2*M)) }
//    c(s)=sqrt(2.0/M) ,s!=0
//        =sqrt(1.0/M) ,s==0
//
void dct_basic(const double* src, const int nsize, Wy_Array<double>& dest)
{
const double r=M_PI/(2.0*nsize);
const double a=sqrt(2.0/nsize);

 if(nsize<=4) {
   switch(nsize) {
     case 1: dest.push_back(src[0]);
             break;
     case 2: dest.push_back(M_SQRT1_2*(src[0]+src[1]));
             dest.push_back(src[0]*cos(M_PI/4.0)+src[1]*cos(M_PI*3.0/4.0));
             break;
     case 3: dest.push_back(M_SQRT1_2*a*(src[0]+src[1]+src[2]));
             dest.push_back(a*(src[0]*cos(M_PI    /6.0)
                              +src[1]*cos(M_PI*3.0/6.0)
                              +src[2]*cos(M_PI*5.0/6.0)));
             dest.push_back(a*(src[0]*cos(M_PI*2.0/6.0)
                              +src[1]*cos(M_PI)
                              +src[2]*cos(M_PI*10.0/6.0)));
             break;
     case 4: dest.push_back(M_SQRT1_2*a*(src[0]+src[1]+src[2]+src[3]));
             dest.push_back(a*(src[0]*cos(M_PI     /8.0)
                              +src[1]*cos(M_PI* 3.0/8.0)
                              +src[2]*cos(M_PI* 5.0/8.0)
                              +src[3]*cos(M_PI* 7.0/8.0)));
             dest.push_back(a*(src[0]*cos(M_PI* 2.0/8.0)
                              +src[1]*cos(M_PI* 6.0/8.0)
                              +src[2]*cos(M_PI*10.0/8.0)
                              +src[3]*cos(M_PI*14.0/8.0)));
             dest.push_back(a*(src[0]*cos(M_PI* 3.0/8.0)
                              +src[1]*cos(M_PI* 9.0/8.0)
                              +src[2]*cos(M_PI*15.0/8.0)
                              +src[3]*cos(M_PI*21.0/8.0)));
             break;
   }
   return;
 }

int s,m;
 for(s=0; s<nsize; ++s) {
   const double rs=r*s;
   double ftmp=0;
   for(m=0; m<nsize; ++m) {
     ftmp+=src[m]*cos(rs*((m<<1)+1));
   }
   if(s==0) {
     dest.push_back(ftmp*a*M_SQRT1_2);
   } else {
     dest.push_back(ftmp*a);
   }
 }
}

// Inverse DCT of N^2 implement
//
void idct_basic(const double* src, const int nsize, Wy_Array<double>& dest)
{

 if(nsize<=0) {
   return;
 }

const double r=M_PI/(2.0*nsize);
const double a=sqrt(2.0/nsize);
double ftmp;
int s,m;

 for(m=0; m<nsize; ++m) {
   ftmp=0;
   for(s=0; s<nsize; ++s) {
     if(s==0) {
       ftmp+=M_SQRT1_2*src[s]*cos(r*s*((m<<1)+1));
     } else {
       ftmp+=src[s]*cos(r*s*((m<<1)+1));
     }
   }
   dest.push_back(ftmp*a);
 }
}

// Main subroutine of NlogN algorithm
//  src size must be power of 2
//
static void dct_fast2n_1(const double* src, const int N, Wy_Array<double>& dest)
{
int i,j;

 if((N<=7)||(N&1)) {
   switch(N) {
     case 1: dest.push_back(src[0]);
             break;
     case 2: dest.push_back(src[0]+src[1]);
             dest.push_back(src[0]*cos(M_PI/4.0)+src[1]*cos(M_PI*5.0/4.0));
             break;
     case 3: dest.push_back(src[0]+src[1]+src[2]);
             dest.push_back(src[0]*cos(M_PI     /6.0)
                           +src[1]*cos(M_PI* 5.0/6.0)
                           +src[2]*cos(M_PI* 9.0/6.0));
             dest.push_back(src[0]*cos(M_PI* 2.0/6.0)
                           +src[1]*cos(M_PI*10.0/6.0)
                           +src[2]*cos(M_PI*18.0/6.0));
             break; 
     case 4: dest.push_back(src[0]+src[1]+src[2]+src[3]);
             dest.push_back(src[0]*cos(M_PI     /8.0)
                           +src[1]*cos(M_PI* 5.0/8.0)
                           +src[2]*cos(M_PI* 9.0/8.0)
                           +src[3]*cos(M_PI*13.0/8.0));
             dest.push_back(src[0]*cos(M_PI* 2.0/8.0)
                           +src[1]*cos(M_PI*10.0/8.0)
                           +src[2]*cos(M_PI*18.0/8.0)
                           +src[3]*cos(M_PI*26.0/8.0));
             dest.push_back(src[0]*cos(M_PI* 3.0/8.0)
                           +src[1]*cos(M_PI*15.0/8.0)
                           +src[2]*cos(M_PI*27.0/8.0)
                           +src[3]*cos(M_PI*39.0/8.0));
             break;
     default:
          const double r=M_PI/(2.0*N);
          for(i=0; i<N; ++i) {
            const double rs=r*i;
            double ftmp=0;
            for(j=0; j<N; ++j) {
              ftmp+=src[j]*cos(rs*((j<<2)+1));
            }
            dest.push_back(ftmp);
          }
   }
   return;
 }
 
const int N2=N>>1;
Wy_Array<double> buf,Y0,Y1;

 buf.resize(N2);
 for(i=0,j=N2; j<N; ++i,++j) {
   buf[i]=src[i]+src[j];
 }
 dct_fast2n_1(buf.begin(),N2,Y0);

 for(i=0,j=N2; j<N; ++i,++j) {
   buf[i]=(src[i]-src[j])*2.0*cos((M_PI*2.0*i+M_PI_2)/N);
 }
 dct_fast2n_1(buf.begin(),N2,Y1);

 buf[0]=Y1[0]/2.0;
 for(i=1; i<N2; ++i) {
   buf[i]=Y1[i]-buf[i-1];
 }

 for(i=0; i<N2; ++i) {
   dest.push_back(Y0[i]);
   dest.push_back(buf[i]);
 }
}

// DCT of NlogN algorithm
//
void dct(const double* src, const int N, Wy_Array<double>& dest)
{
 
 if((N<8)||(N&1)) {
   dct_basic(src,N,dest);
   return;
 }

Wy_Array<double> g;
int i,j;

 g.resize(N);
 for(i=0,j=0; j<N; ++i,j+=2) {
   g[i]=src[j];
 }
 for(  j=N-1; i<N; ++i,j-=2) {
   g[i]=src[j];
 }
 dct_fast2n_1(g.begin(),N,dest);

 // Multiply coefficient
 // As dest is defined to be pushed (size may not be N), index dest items
 // from the back.
 //
const double a=sqrt(2.0/N);
 for(i=N-1,j=dest.size()-1; i>=0; --i,--j) {
   dest[j]*=a;
 }
 dest[j+1]*=M_SQRT1_2;
}

//---------------------
// Test
//---------------------

void dump(const Wy_Array<double>& f)
{
 for(size_t i=0; i<f.size(); ++i) {
   Wy::cout << f[i] << ", ";
 }
 Wy::cout << '\n';
}

// Return 0  : f1==f2
//        >0 : index where f1[i]!=f2[i] (differs by thresh)
//        <0 : size not identical
int cmp(const Wy_Array<double>& f1, const Wy_Array<double>& f2)
{
//const double thresh=1e-10;   // fail on size=398
const double thresh=1e-9;    // fail on size=
//const double thresh=1e-8;      
const size_t N=f1.size();
 if(N!=f2.size()) {
   return(-1);
 }
 for(size_t i=0; i<N; ++i) {
   if(fabs(f1[i]-f2[i])>thresh) {
     return(i);
   }
 }
 return(0);
}

Wy_Array<double> f,F1,G1;

static void t1(void)
{
const int tsiz=256;
const int tfreq=1000;
int i;
 f.reset();
 for(i=0; i<tsiz; ++i) {
   f.push_back((random()-(RAND_MAX>>1))    // integer [-RAND_MAX/2, RAND_MAX/2]
            *100.0/(double)(RAND_MAX>>1)); // scale to [-100,+100]
 }
 for(i=0; i<tfreq; ++i) {
   F1.reset();
   //dct_basic(f.begin(),f.size(),F1);
   dct(f.begin(),f.size(),F1);
 }
 f.reset();
}

static void t2(void)
{
int tsiz;
int i,k;
 
 // make samples of size 1-99
 //
 for(tsiz=1; tsiz<99; ++tsiz) {
   f.reset();
   for(i=0; i<tsiz; ++i) {
     f.push_back((random()-(RAND_MAX>>1))  // integer [-RAND_MAX/2, RAND_MAX/2]
            *100.0/(double)(RAND_MAX>>1)); // scale to [-100,+100]
   }
   // check dct_basic/idct consistency
   F1.reset();
   G1.reset();
   dct_basic(f.begin(),f.size(),F1);
   idct_basic(F1.begin(),F1.size(),G1);
   k=cmp(f,G1);
   if(k!=0) {
     Wy::cout << "sample size=" << tsiz << '\n';
     if(k<0) {
       Wy::cout << "dct_basic()/idct_basic() output different sizes\n";
     } else {
       Wy::cout << "dct_basic[i]!=idct_basic[i], i==" << k << '\n';
     }
     Wy::cout << "sample[]=\n";
     dump(f);
     Wy::cout << "dct_basic[]=\n";
     dump(F1);
     Wy::cout << "idct_basic[]=\n";
     dump(G1);
     if(k>0) {
       Wy::cout << "sample[i]=" << f[k] << "  dct_basic[i]-idct_basic[i]=" << G1[k] << '\n';
     }
     return;
   }
   G1.reset();
   dct(f.begin(),f.size(),G1);
   k=cmp(F1,G1);
   if(k!=0) {
     Wy::cout << "sample size=" << tsiz << '\n';
     if(k<0) {
       Wy::cout << "dct_basic()/dct() output different sizes\n";
     } else {
       Wy::cout << "dct_basic[i]!=dct[i], i==" << k << '\n';
     }
     Wy::cout << "sample[]=\n";
     dump(f);
     Wy::cout << "dct_basic[]=\n";
     dump(F1);
     Wy::cout << "dct[]=\n";
     dump(G1);
     if(k>0) {
       Wy::cout << "sample[i]=" << f[k] << "  dct_basic[i]=" << F1[k] << "  dct[i]=" << G1[k] << '\n';
     }
     return;
   }
 }
}

int main(void)
try {
 t1();
 t2();
 Wy::cout << "Ok\n";
}
catch(const WyRet& e) {
 Wy::cerr << Wy::wrd(e) << "\n";
 return(-1);
}
catch(...) {
 Wy::cerr << "unknown stack unwind\n";
 return(-1);
};