dct-iv.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);
		i1=i1-n1;
		BTRVS(C,i1+1,n1);
	}
}

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

void 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 FWT4(int m, int n)        /*  DWT-IV  */
{
	double t;
	int i,n1,n2;
	COEF(n);
	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];
	}
}

void DCTIV(int m,int n)
{
	int k;
	for (k=n+1;k<=2*n;k++)
		F[k]=0;
	FWT4(m+1,2*n);
	for (k=1;k<=n;k++)
		F[k]=(F[k]-F[2*n+1-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 );
	DCTIV(M,N);
	_ftime(&timebuffer2 );
	ComMs=(timebuffer2.time-timebuffer.time)+float(timebuffer2.millitm-timebuffer.millitm)/1000;
	for (i=1;i<=N;i++)
		printf("x[%d]=%f\n",i,F[i]);
	printf("运算时间为%f秒。\n",ComMs);
	getch();
}