一维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==16if(C!=NULL) delete []C;C=(double *)new double[total];halftotal=total>>1;  //右移一位相当于给该数除以2 halftotal == 8for(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 == 8int len=endstart>>1; //len==4factor=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);//111F[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;}