www.pudn.com > AudioWave.rar > DCT.cpp


// DCT.cpp: implementation of the CDCT class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "Audio.h"
#include "DCT.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif
#include "conio.h"
#include "stdio.h"
#include "math.h"
#include "stdlib.h"
#include "sys/timeb.h"

#define PI 3.14159265359

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CDCT::CDCT()
{

}

CDCT::~CDCT()
{

}

void CDCT::BTRVS(double a[],int W,int n)
{
double t;
int n1,i,j,k;
j=1;
n1=n-1;
for (i=1;i<=n1;i++)
{
if (i<j)
{
t=a[j+W-1];
a[j+W-1]=a[i+W-1];
a[i+W-1]=t;
}
k=n>>1;
while (k<j)
{
j=j-k;
k=k>>1;
}
j=j+k;
}
}

void CDCT::COEF(int n) /* 求变换系数　*/
{
double p,A;
int n3,i1,n1,i,i2,n2;
C[1]=sqrt(2.0);
if (n>2)
{
n3=n>>1;
i1=1;
n1=1;
do
{
n1=n1<<1;
n2=n1<<2;
p=PI/n2;
for (i=1;i<=n1;i++)
{
i2=(i<<1)-1;
A=cos(i2*p);
A=A*2;
C[i1+i]=1.0/A;
}
BTRVS(C,i1+1,n1);
i1=i1+n1;
}
while (i1< n3);
i1=i1-n1;
BTRVS(C,i1+1,n1);
}
}

void CDCT::FWT3(int m,int n) /* DWT-III */
{
double t;
int i0,i,i1,j,j1,n1,k1,k2,n2,n3,n4,i2,i3,n5,n6,k,n0,j2;
if (n!=4)
{
i0=1;
for (i=3;i<=m;i++)
{
i1=i0;
i0=i1<<1;
for (j=1;j<=i1;j++)
{
j1=1-j;
n1=n+j1;
k1=i0+j1;
k2=k1+i0;
t=F[n1];
for (k=n1;k>=k2;k=k-i0)
F[k]=F[k]+F[k-i0];
F[k1]=F[k1]-t;
}
}
}
n1=n>>2;
n2=n1<<1;
n3=n1+n2;
n4=n2;
for (i=1;i<=n1;i++)
{
i1=n1+i;
i2=n2+i;
i3=n3+i;
t=F[i];
F[i]=t+F[i2];
F[i2]=t-F[i2];
F[i1]=C[1]*F[i1];
F[i3]=C[1]*F[i3];
}
if (n!=4)
{
i1=2;
n0=1;
for (i=3;i<=m;i++)
{
n0=n0<<1;
n5=n0<<1;
n3=n/(n5+n5);
n2=n3<<1;
n1=n2<<1;
n5=-n1;
for (j=1;j<=n0;j++)
{
n5=n1+n5;
n6=n5+n2;
for (k=1;k<=n2;k++)
{
k1=n5+k;
k2=n6+k;
t=F[k1];
F[k1]=t+F[k2];
F[k2]=t-F[k2];
}
}
i2=i1+n0;
i3=i1-1;
n5=n3-n1;
for (j=1;j<=n0;j++)
{
n5=n5+n1;
n6=n5+n2;
j1=i3+j;
j2=i2-j;
for (k=1;k<=n3;k++)
{
k1=n5+k;
k2=n6+k;
F[k1]=C[j1]*F[k1];
F[k2]=-C[j2]*F[k2];
}
}
i1=i2;
}
}
for (i=1;i<=n4;i++)
{
i1=i<<1;
t=F[i1-1];
F[i1-1]=t+F[i1];
F[i1]=t-F[i1];
}
BTRVS(F,1,n);
}

void CDCT::FWT4(int m, int n) /* DWT-IV */
{
double t;
int i,n1,n2;
t=F[n];
for (i=n;i>=2;i--)
F[i]=F[i]+F[i-1];
F[1]=F[1]-t;
FWT3(m,n);
n1=n>>1;
n2=n1-1;
for (i=1;i<=n1;i++)
{
F[i]=C[n2+i]*F[i];
F[i+n1]=-C[n-i]*F[i+n1];
}
}

int CDCT::log2(int number) /* 求N的幂 */
{
int i;
i=-1;
while(1)
{
if(number==0)break;
number=number>>1;
i++;
}
return i;
}

void CDCT::DCTIV(float *fData, int nPower)
{
int i;
int M;
N = 1<<nPower;
double sqrtN = sqrt(N);
F = new double[N+1];
C = new double[N+1];
for (i=1;i<=N;i++)
{
F[i] = fData[i-1];
}
M=log2(N);
COEF(N);
FWT4(M,N);
for (i=1;i<=N;i++)
{
F[i]=F[i]/sqrtN;
fData[i-1] = (float)F[i];
}
delete[] F;
delete[] C;
}

void CDCT::IDCTIV(float *fData, int nPower)
{
int i;
int M;
N = 1<<nPower;
double sqrtN = sqrt(N);
F = new double[N+1];
C = new double[N+1];
for (i=1;i<=N;i++)
{
F[i] = fData[i-1];
}
M=log2(N);
COEF(N);
FWT4(M,N);
for (i=1;i<=N;i++)
{
F[i]=F[i]/sqrtN;
fData[i-1] = (float)F[i];
}
delete[] F;
delete[] C;
}