FDCT快速算法

1、浮点FDCT快速变换c程序代码最近在做FDCT快速算法的项目，网上找了很多资料，发现基本上都有谬错，几经周折，终于自己写出了一个基于FFT的二维行-列快速变换的FDCT，经过反复测试，该程序运行良好，输出结果与标准值无二。下面是源代码以及详细算法推导过程。#define DCTSIZE 8typedef float FAST_FLOAT;void jpeg_fdct_float (FAST_FLOAT * lpBuff) FAST_FLOAT tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7; FAST_FLOAT tmp10, tmp11, tmp

2、12, tmp13; FAST_FLOAT * dataptr; FAST_FLOAT A,B,C,D; FAST_FLOAT V4=0,0,0,0; int ctr,i; /* 第一部分，对行进行计算 */ dataptr = lpBuff; for (ctr = DCTSIZE-1; ctr >= 0; ctr-) tmp0 = dataptr0 + dataptr7; tmp7 = dataptr0 - dataptr7; tmp1 = dataptr1 + dataptr6; tmp6 = dataptr1 - dataptr6; tmp2 = dataptr2 + datapt

3、r5; tmp5 = dataptr2 - dataptr5; tmp3 = dataptr3 + dataptr4; tmp4 = dataptr3 - dataptr4; /* 对偶数项进行运算 */ tmp10 = tmp0 + tmp3; /* phase 2 */ tmp13 = tmp0 - tmp3; tmp11 = tmp1 + tmp2; tmp12 = tmp1 - tmp2; dataptr0 = tmp10 + tmp11; /* phase 3 */ dataptr4 = (tmp10 - tmp11)*(FAST_FLOAT)0.707106781); datapt

4、r2=tmp13*(FAST_FLOAT)0.923879000)+ tmp12*(FAST_FLOAT)0.382683000); /* phase 5 */ dataptr6=tmp13*(FAST_FLOAT)0.382683000) - tmp12*(FAST_FLOAT)0.923879000); /* 对奇数项进行计算 */*构造过渡表达式，v(n)=(x(n)-x(7-n)*2cos(/16)(2n+1),先计算v(n)的DCT，再通过过渡变换因子转换为X2k+1的表达式。对偶数项X1,X5进行运算，按照之前的推论，可以得出*/ A=2*tmp7*(FAST_FLOAT)0.98

5、0785283)+2*tmp4*(FAST_FLOAT)0.195090322); B=2*tmp6*(FAST_FLOAT)0.831469612)+2*tmp5*(FAST_FLOAT)0.555570233);V0=A+B;V2= (FAST_FLOAT)0.707106781)*(A-B);/*对奇数项X3,X7进行计算，按照2点第四类DCT算法，可以得出*/C=2*tmp7*(FAST_FLOAT)0.980785283)-2* tmp4*(FAST_FLOAT) 0.195090322);D=2*tmp6*(FAST_FLOAT)0.831469612)-2* tmp5*(FAST

6、_FLOAT)0.555570233);V1=C*(FAST_FLOAT)0.923879532)+D*(FAST_FLOAT)0.382683432);V3=C*(FAST_FLOAT)0.382683432)-D*(FAST_FLOAT)0.923879532);dataptr1=(FAST_FLOAT)0.500000000)*V0;dataptr3=V1-dataptr1;dataptr5=V2- dataptr3;dataptr7=V3- dataptr5; for(i=0;i<8;i+)if (i=0) dataptri=( (FAST_FLOAT)0.353553390)*

7、dataptri;else dataptri=( (FAST_FLOAT)0.500000000)*dataptri; dataptr += DCTSIZE; /* 将指针指向下一行 */ /* 第二部分，对列进行计算 */ dataptr = lpBuff; for (ctr = DCTSIZE-1; ctr >= 0; ctr-) tmp0 = dataptrDCTSIZE*0 + dataptrDCTSIZE*7; tmp7 = dataptrDCTSIZE*0 - dataptrDCTSIZE*7; tmp1 = dataptrDCTSIZE*1 + dataptrDCTSIZE

8、*6; tmp6 = dataptrDCTSIZE*1 - dataptrDCTSIZE*6; tmp2 = dataptrDCTSIZE*2 + dataptrDCTSIZE*5; tmp5 = dataptrDCTSIZE*2 - dataptrDCTSIZE*5; tmp3 = dataptrDCTSIZE*3 + dataptrDCTSIZE*4; tmp4 = dataptrDCTSIZE*3 - dataptrDCTSIZE*4; /* 对偶数项进行运算 */ tmp10 = tmp0 + tmp3; /* phase 2 */ tmp13 = tmp0 - tmp3; tmp11

9、 = tmp1 + tmp2; tmp12 = tmp1 - tmp2; dataptrDCTSIZE*0 = tmp10 + tmp11; /* phase 3 */ dataptrDCTSIZE*4 = (tmp10 - tmp11)*(FAST_FLOAT)0.707106781); dataptrDCTSIZE*2=tmp13*(FAST_FLOAT)0.923879000)+ tmp12*(FAST_FLOAT)0.382683000); /* phase 5 */ dataptrDCTSIZE*6=tmp13*(FAST_FLOAT)0.382683000)- - tmp12*(F

10、AST_FLOAT)0.923879000); /* 对奇数项进行计算 */*构造过渡表达式，v(n)=(x(n)-x(7-n)*2cos(/16)(2n+1),先计算v(n)的DCT，再通过过渡变换因子转换为X2k+1的表达式。对偶数项X1,X5进行运算，按照之前的推论，可以得出*/ A=2*tmp7*(FAST_FLOAT)0.980785283)+2*tmp4*(FAST_FLOAT)0.195090322); B=2*tmp6*(FAST_FLOAT)0.831469612)+2*tmp5*(FAST_FLOAT)0.555570233);V0=A+B;V2=( (FAST_FLOAT

11、) 0.707106781)*(A-B);/*对奇数项X3,X7进行计算，按照2点第四类DCT算法，可以得出*/C=2* tmp7*(FAST_FLOAT)0.980785283)-2* tmp4*(FAST_FLOAT)0.195090322);D=2*tmp6*(FAST_FLOAT)0.831469612)-2* tmp5*(FAST_FLOAT)0.555570233);V1=C*(FAST_FLOAT)0.923879532)+D*(FAST_FLOAT)0.382683432);V3=C*(FAST_FLOAT)0.382683432)-D*(FAST_FLOAT)0.923879532);dataptr1* DCTSIZE=( (FAST_FLOAT)0.500000000)*V0;dataptr3* DCTSIZE=V1-dataptr1* DCTSIZE;dataptr5* DCTSIZE=V2- dataptr3* DCTSIZE;dataptr7* DCTSIZE=V3- dataptr5* DCTSIZE; fo