dct-ii.cpp

一般的快速DCT算法只有8点和16点的

#include "conio.h"
#include "stdio.h"
#include "math.h"
#include "stdlib.h"
#include "sys/timeb.h"

#define PI 3.14159265359
#define N 8    /*  N 为2的幂  */

double F[2*N+1],C[2*N+1],y[2*N+1];         /*  定义全局变量  */

void 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 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);
	}
}

void FWT2(int m, int n)         /*  DWT-II  */
{
	double t; 
	int n1,n2,n0,n3,n4,n5,i,j,k,i1,i2,j2,j3,j4,k1,k2,k3,k4,j1;
	COEF(n);
	BTRVS(F,1,n);
	n1=1;
	for (i=2;i<=m;i++)
	{
		n0=n1;
		n1=n1<<1;
		n2=n1<<1;
		j2=-n2;
		n3=n/n2;      
		n4=n3<<1;
		n5=n3-1;
		for (j=1;j<=n3;j++)
		{
			j2=j2+n2;
			j3=n5+j;
			j4=n4-j;
			for (k=1;k<=n0;k++)
			{
				k1=j2+k;
				k2=k1+n0;
				k3=k1+n1;
				k4=k2+n1;
				t=F[k1];
				F[k1]=t+F[k2];
				F[k2]=C[j3]*(t-F[k2]);
				t=F[k3];
				F[k3]=t+F[k4];
				F[k4]=-C[j4]*(t-F[k4]);
			}
		}
	}
	n2=n1+n0;
	for (i=1;i<=n0;i++)
	{
		i1=n1+i;
		i2=n2+i;
		F[i2]=-F[i2];
		t=F[i];
		F[i]=t+F[i1];
		F[i1]=t-F[i1];
	}
	if (n!=4)
	{
		for (i=3;i<=m;i++) 
		{
			n1=n0;
			n0=n0>>1;
			n2=n-n0;
			for (j=1;j<=n0;j++)
			{
				j1=j+n0;
				j2=j+n2;
				j3=j2-n1;
				t=F[j1];
				for (k=j1;k<=j3;k=k+n1)
					F[k]=F[k]+F[k+n1];
				F[j2]=F[j2]-t;
			}
		}
	}
}

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

void DCTII(int m,int n)
{	int k;
	for (k=n+1;k<=2*n;k++)
	F[k]=0;
	FWT2(m+1,2*n);
	for (k=2;k<=n;k++)
		F[k]=(F[k]-F[2*n+2-k])/2;
}

void main()            /*  主函数　*/
{
	int i,M;
	double ComMs;
	struct _timeb timebuffer,timebuffer2;
	F[0]=0;
  
  /*  随机生成输入序列 */
	/*randomize();*/
	for (i=1;i<=N;i++)
		F[i]=i-1;
 
 
  /*  DWT 过程  */
	M=log2(N);
	_ftime(&timebuffer );
	DCTII(M,N);
	DCTII(M,N);
	_ftime(&timebuffer2 );
	ComMs=(timebuffer2.time-timebuffer.time)+float (timebuffer2.millitm-timebuffer.millitm)/1000;
	for (i=1;i<=N;i++)
		printf("F[%d]=%f\n",i,F[i]);
    printf("运算时间为%f秒。\n",ComMs);

	getch();
}