天天看點

一維2N浮點FDCT變換以及反變換

一維2N浮點FDCT變換

1、 fdct.h檔案

#ifndef __JENC__  //判斷宏定義是否被定義
#define __JENC__

#define PI 3.1415926
#include <math.h>
#include <string>
double *C=NULL;

class Fdct
{
public:
	

	void Invoke()
	{	
	double dataf[16] = {137,136, 133, 136, 138, 134, 134, 132, 132, 138, 129, 131, 132, 127, 129, 128};
	fdct_1D(dataf, GetTwoIndex(16));
	for(int i = 0; i < 16; i++){
		printf("%f\n",dataf[i]);
	}
	}

private:

  //*******************************************************************
 //方法名稱:GetTwoIndex
 //說明:算出deg的值
 //日期:2012.2.29
 //參數說明:
 //num: 資料的數量
 //*******************************************************************
int GetTwoIndex(int num)
{
	int temp = 0;
	if(num<=0) return -1;
	while(1<<temp<num)
		temp++;
	return temp;
}

 //*******************************************************************
 //方法名稱:swap
 //說明:交換temp1與temp2的值
 //日期:2012.2.29
 //參數說明:
 //a:所要交換的資料
 //b:所要交換的資料
 //*******************************************************************
void swap(double &a,double &b)
{
	double temp;
	temp=a;
	a=b;
	b=temp;
}

 //*******************************************************************
 //方法名稱:bitrev
 //說明:位反序函數
 //日期:2012.2.29
 //參數說明:
 //bi:不清楚
 //deg:為DCT變換資料的長度
 //*******************************************************************
	int bitrev(int bi,int deg)
{
	int j=1,temp=0,degnum,halfnum;
	degnum=deg;
	//if(deg<0) return 0;
	if(deg==0) return bi;
	halfnum=1<<(deg-1);
	while(halfnum)
	{
		if(halfnum&bi)
			temp+=j;
		j<<=1;
		halfnum>>=1;
	}
	return temp;	
}
	
//*******************************************************************
 //方法名稱:fbitrev
 //說明: 實作位反序排列(注意和上面函數的差別)
 //日期:2012.2.29
 //參數說明:
 //f:不清楚
 //deg:為DCT變換資料的長度
 //*******************************************************************
void fbitrev(double *f,int deg)
{    
	if(deg==1) return;
	int len=(1<<deg)-1;
	int i=1;
	int ii;
	double temp;
	while(i<len)
	{   ii=bitrev(i,deg); 
	if(ii>i)
	{  temp=f[ii];
		  f[ii]=f[i];
		  f[i]=temp;
	}
	i++;
	}
}
 //*******************************************************************
 //方法名稱:initDCTParam
 //說明:考慮到DCT變換中的系數要重複計算,可使用查找表來提高運作的效率,隻要開始DCT變換之前計算一次,DCT變換中就可以隻查找而無需計算系數。
 //日期:2012.2.29
 //參數說明:
 //deg:為DCT變換資料的長度的幂
 //*******************************************************************
void initDCTParam(int deg) 
{

	int total,halftotal,i,group,endstart,factor;
	total=1<<deg; //total==16
	if(C!=NULL) delete []C;
	C=(double *)new double[total];
	halftotal=total>>1;  //右移一位相當于給該數除以2 halftotal == 8
	for(i=0;i<halftotal;i++)
		C[total-i-1]=(double)(2*i+1);

	for(group=0;group<deg-1;group++)
	{ 
		endstart=1<<(deg-1-group); //endstart == 8
		int len=endstart>>1; //len==4
		factor=1<<(group+1);
		for(int j=0;j<len;j++)
			C[endstart-j-1]=factor*C[total-j-1];
	}
	for(i=1;i<total;i++)
		C[i]=2.0*cos(C[i]*PI/(total<<1));	//C[0]空着,沒使用
}

 //*******************************************************************
 //方法名稱:dct_forward
 //說明:DCT變換過程可分為兩部分:前向追底和後向回根,這個是實作前向追底。
 //日期:2012.2.29
 //參數說明:
 //f:存儲DCT資料
 //deg:為DCT變換資料的長度
 //*******************************************************************

void dct_forward(double *f,int deg)
{
	int i_deg,i_halfwing,total,wing,wings,winglen,halfwing;
	double temp1,temp2;
	total=1<<deg;
	for(i_deg=0;i_deg<deg;i_deg++)
	{
		wings=1<<i_deg;
		winglen=total>>i_deg;
		halfwing=winglen>>1;
		for(wing=0;wing<wings;wing++)
		{
			for(i_halfwing=0;i_halfwing<halfwing;i_halfwing++)
			{
				temp1=f[wing*winglen+i_halfwing];
				temp2=f[(wing+1)*winglen-1-i_halfwing];
				if(wing%2)
					swap(temp1,temp2);  // 交換temp1與temp2的值
				f[wing*winglen+i_halfwing]=temp1+temp2;
				f[(wing+1)*winglen-1-i_halfwing]=(temp1-temp2)*C[winglen-1-i_halfwing];
			}
		}
	}
}

 //*******************************************************************
 //方法名稱:dct_backward
 //說明:DCT變換過程可分為兩部分:前向追底和後向回根,這個是實作後向回根。
 //日期:2012.2.29
 //參數說明:
 //f:存儲DCT資料
 //deg:為DCT變換資料的長度
 //*******************************************************************
void dct_backward(double *f,int deg)
{
	int total,i_deg,wing,wings,halfwing,winglen,i_halfwing,temp1,temp2;
	total=1<<deg;
	for(i_deg=deg-1;i_deg>=0;i_deg--)
	{
		wings=1<<i_deg;
		winglen=1<<(deg-i_deg);
		halfwing=winglen>>1;
		for(wing=0;wing<wings;wing++)
		{
			for(i_halfwing=0;i_halfwing<halfwing;i_halfwing++)
			{  //f[i_halfwing+wing*winglen]=f[i_halfwing+wing*winglen];
				if(i_halfwing==0)
					f[halfwing+wing*winglen+i_halfwing]=0.5*f[halfwing+wing*winglen+i_halfwing];
				else
				{
					temp1=bitrev(i_halfwing,deg-i_deg-1); // bitrev為位反序
					temp2=bitrev(i_halfwing-1,deg-i_deg-1); // 第一參數為要變換的數
					// 第二參數為二進制長度
					f[halfwing+wing*winglen+temp1]=f[halfwing+wing*winglen+temp1]-f[halfwing+wing*winglen+temp2];
				}	
			}
		}
	}
}

//*******************************************************************
 //方法名稱:fdct_1D_no_param
 //說明:完整實作一維DCT變換
 //日期:2012.2.29
 //參數說明:
 //f:存儲DCT資料
 //deg:為DCT變換資料的長度
 //*******************************************************************
void fdct_1D_no_param(double *f,int deg) //deg == 4
{
	
	initDCTParam(deg);
	dct_forward(f,deg);
	dct_backward(f,deg);
	fbitrev(f,deg); // 實作位反序排列
	f[0]=1/(sqrt(2.0))*f[0];
}

void fdct_1D(double *f,int deg) 
{
	fdct_1D_no_param(f,deg); 
	int total=1<<deg; 
	double param=sqrt(2.0/total); 
	for(int i=0;i<total;i++){
		f[i]=param*f[i];
	}	
}
};
#endif // __JENC__
           

2、FDCT.c檔案

// NewFDCT.cpp : 定義控制台應用程式的入口點。

#include "stdafx.h"
#include "fdct.h" //引用頭檔案

//實作FDCT一維快速變換
int _tmain(int argc, _TCHAR* argv[])
{
	Fdct fdct;
	fdct.Invoke();
	return 0;
}
           

一維2N浮點IFDCT變換

1、IFDCT.h

#ifndef __JENC__  //判斷宏定義是否被定義
#define __JENC__

#define PI 3.1415926
#include <math.h>
#include <string>
double *C=NULL;

class Fidct
{
public:

	void Invoke()
	{
	double dataf[16] = {531.500000, 10.491843, -2.369268, 1.307582, -0.718939, -0.948647, 2.013496, 3.625035, 1.500000, -1.973283, 2.122850, -0.347127, -5.655363, 2.111544, 0.908794, -1.302960};
	fidct_1D(dataf, 4);
	for(int i = 0; i < 16; i++){
		printf("%d\n",int(dataf[i]));
	}
	}

private:
  //*******************************************************************
 //方法名稱:initIDCTParam
 //說明:考慮到IDCT變換中的系數要重複計算,可使用查找表來提高運作的效率,隻要開始IDCT變換之前計算一次,IDCT變換中就可以隻查找而無需計算系數。
 //日期:2012.2.29
 //參數說明:
 //deg:為IDCT變換資料長度的幂
 //*******************************************************************
void initIDCTParam(int deg)
{
	int total,halftotal,i,group,endstart,factor;
	total=1<<deg;
	if(C!=NULL) delete []C;
	C=(double *)new double[total];
	halftotal=total>>1;
	for(i=0;i<halftotal;i++)
		C[total-i-1]=(double)(2*i+1);

	for(group=0;group<deg-1;group++)
	{ 
		endstart=1<<(deg-1-group);
		int len=endstart>>1;
		factor=1<<(group+1);
		for(int j=0;j<len;j++)
			C[endstart-j-1]=factor*C[total-j-1];
	}
	for(i=1;i<total;i++)
		C[i]=1.0/(2.0*cos(C[i]*PI/(total<<1)));	//C[0]空着,沒使用

}

 //*******************************************************************
 //方法名稱:bitrev
 //說明:位反序函數
 //日期:2012.2.29
 //參數說明:
 //bi:不清楚
 //deg:為IDCT變換資料長度的幂
 //*******************************************************************
int bitrev(int bi,int deg)
{
	int j=1,temp=0,degnum,halfnum;
	degnum=deg;
	//if(deg<0) return 0;
	if(deg==0) return bi;
	halfnum=1<<(deg-1);
	while(halfnum)
	{
		if(halfnum&bi)
			temp+=j;
		j<<=1;
		halfnum>>=1;
	}
	return temp;	
}

//*******************************************************************
 //方法名稱:fbitrev
 //說明: 實作位反序排列(注意和上面函數的差別)
 //日期:2012.2.29
 //參數說明:
 //f:不清楚
 //deg:為DCT變換資料長度的幂
 //*******************************************************************
void fbitrev(double *f,int deg)
{    
	if(deg==1) return;
	int len=(1<<deg)-1;
	int i=1;
	int ii;
	double temp;
	while(i<len)
	{   ii=bitrev(i,deg); 
	if(ii>i)
	{  temp=f[ii];
	f[ii]=f[i];
	f[i]=temp;
	}
	i++;
	}
}

 //*******************************************************************
 //方法名稱:idct_forward
 //說明:IDCT變換過程可分為兩部分:前向追底和後向回根,這個是實作前向追底。
 //日期:2012.2.29
 //參數說明:
 //f:存儲IDCT資料
 //deg:為IDCT變換資料長度的幂
 //*******************************************************************
void idct_forward(double *F,int deg)
{
	int total,i_deg,wing,wings,halfwing,winglen,i_halfwing,temp1,temp2;
	total=1<<deg;
	for(i_deg=0;i_deg<deg;i_deg++)
	{
		wings=1<<i_deg;
		winglen=1<<(deg-i_deg);
		halfwing=winglen>>1;
		for(wing=0;wing<wings;wing++)
		{
			for(i_halfwing=halfwing-1;i_halfwing>=0;i_halfwing--)
			{
				if(i_halfwing==0)
					F[halfwing+wing*winglen+i_halfwing]=2.0*F[halfwing+wing*winglen+i_halfwing];
				else
				{ 
					temp1=bitrev(i_halfwing,deg-i_deg-1); // bitrev為位反序
					temp2=bitrev(i_halfwing-1,deg-i_deg-1);// 第一參數為要變換的數
					// 第二參數為二進制長度
					F[halfwing+wing*winglen+temp1]=F[halfwing+wing*winglen+temp1]+F[halfwing+wing*winglen+temp2];
				}
			}
		}
	}
} 

 //*******************************************************************
 //方法名稱:idct_backward
 //說明:IDCT變換過程可分為兩部分:前向追底和後向回根,這個是實作後向回根。
 //日期:2012.2.29
 //參數說明:
 //F:存儲DCT資料
 //deg:為DCT變換資料長度的幂
 //*******************************************************************
void idct_backward(double *F,int deg)
{
	int i_deg,i_halfwing,total,wing,wings,winglen,halfwing;
	double temp1,temp2;
	total=1<<deg;
	for(i_deg=deg-1;i_deg>=0;i_deg--)
	{
		wings=1<<i_deg;
		winglen=total>>i_deg;
		halfwing=winglen>>1;
		for(wing=0;wing<wings;wing++)
		{
			for(i_halfwing=0;i_halfwing<halfwing;i_halfwing++)
			{
				temp1=F[wing*winglen+i_halfwing];
				temp2=F[(wing+1)*winglen-1-i_halfwing]*C[winglen-1-i_halfwing];
				if(wing%2)
				{
					F[wing*winglen+i_halfwing]=(temp1-temp2)*0.5;
					F[(wing+1)*winglen-1-i_halfwing]=(temp1+temp2)*0.5;
				}
				else
				{
					F[wing*winglen+i_halfwing]=(temp1+temp2)*0.5;
					F[(wing+1)*winglen-1-i_halfwing]=(temp1-temp2)*0.5;
				}

			}
		}
	}

}

//*******************************************************************
 //方法名稱:fidct_1D_no_param
 //說明:完整實作一維DCT變換
 //日期:2012.2.29
 //參數說明:
 //F:存儲DCT資料
 //deg:為DCT變換資料長度的幂
 //*******************************************************************
void fidct_1D_no_param(double *F,int deg)
{
	initIDCTParam(deg);//111
	F[0]=F[0]*sqrt(2.0);
	fbitrev(F,deg); // 實作位反序排列
	idct_forward(F,deg);
	idct_backward(F,deg);

}
void fidct_1D(double *F,int deg)
{
	int total=1<<deg;
	double param=sqrt((double)total/2.0);
	for(int i=0;i<total;i++)
	{
		F[i]=param*F[i];
	}
	fidct_1D_no_param(F,deg);
}
};
#endif // __JENC__
           

2、IFDCT.c

// NewIFDCT.cpp : 定義控制台應用程式的入口點。
//

#include "stdafx.h"
#include "fidct.h" //引用頭檔案




//實作一維反DCT快速變換
int _tmain(int argc, _TCHAR* argv[])
{
	Fidct fidct;
	fidct.Invoke();
	return 0;
}