dct-iii.cpp

在matlab中实现DCT变换

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

#define PI 3.14159265359
#define N 16    /*  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 DCTIII(int m,int n)
{	double t;
	int i0,i,i1,j,j1,n1,k1,k2,n2,n3,n4,i2,i3,n5,n6,k,n0,j2,p;
	for (k=n+1;k<=2*n;k++)
		F[k]=0;
	p=n<<1;
	COEF(p);
	if (p!=4)
	{
		i0=1;
		for (i=3;i<=m+1;i++)
		{
			i1=i0;
			i0=i1<<1;
			for (j=1;j<=i1;j++)
			{
				j1=1-j;
				n1=p+j1;
				k1=i0+j1;
				k2=k1+i0;
				t=F[n1];
				for (k=n1;k>=k2;k=k-i0)
					F[k]+=F[k-i0];
				F[k1]-=t;
			}
		}
	}
	n1=p>>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[i3]*=C[1];
	}
	if (p!=4)
	{
		i1=2;
		n0=1;
		for (i=3;i<=m+1;i++)
		{
			n0=n0<<1;
			n5=n0<<1;
			n3=p/(n5<<1);      
			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[k2]*=(-C[j2]);
				}
			}
			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,p);
	for (k=1;k<=n;k++)
		F[k]=(F[k]+F[2*n+1-k])/2;
}

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

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 );
	DCTIII(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();
}