⭐ 欢迎来到虫虫下载站! | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站

📄 dct_64.cpp

📁 DCT快速算法64位的C语言实现
💻 CPP
📖 第 1 页 / 共 4 页
字号:
/*********************************************************************
*利用zhao等人论文原理,手动展开迭代的DCT变换程序                     *
*以下程序中N=64                                                     *
*其中f为64位的输入流,F为经过DCT变换后的64位输出流                   *                                
**********************************************************************/
#include"stdio.h"
#include"stdlib.h"

#define c_pi_1_128   0.99969881869620     //表示cos(PI*(奇数)/128)的值,以下类似
#define c_pi_3_128   0.99729045667869
#define c_pi_5_128   0.99247953459871
#define c_pi_7_128   0.98527764238894
#define c_pi_9_128   0.97570213003853
#define c_pi_11_128  0.96377606579544
#define c_pi_13_128  0.94952818059304
#define c_pi_15_128  0.93299279883474
#define c_pi_17_128  0.91420975570353
#define c_pi_19_128  0.89322430119552
#define c_pi_21_128  0.87008699110871
#define c_pi_23_128  0.84485356524971
#define c_pi_25_128  0.81758481315158
#define c_pi_27_128  0.78834642762661
#define c_pi_29_128  0.75720884650648
#define c_pi_31_128  0.72424708295147
#define c_pi_33_128  0.68954054473707
#define c_pi_35_128  0.65317284295378
#define c_pi_37_128  0.61523159058063
#define c_pi_39_128  0.57580819141785
#define c_pi_41_128  0.53499761988710
#define c_pi_43_128  0.49289819222978
#define c_pi_45_128  0.44961132965461
#define c_pi_47_128  0.40524131400499
#define c_pi_49_128  0.35989503653499
#define c_pi_51_128  0.31368174039889
#define c_pi_53_128  0.26671275747490
#define c_pi_55_128  0.21910124015687
#define c_pi_57_128  0.17096188876030
#define c_pi_59_128  0.12241067519922
#define c_pi_61_128  0.07356456359967
#define c_pi_63_128  0.02454122852291

#define c_pi_1_64    0.99879545620517     //表示cos(PI*(奇数)/64)的值,以下类似
#define c_pi_3_64    0.98917650996478
#define c_pi_5_64    0.97003125319454
#define c_pi_7_64	 0.94154406518302
#define c_pi_9_64	 0.90398929312344
#define c_pi_11_64	 0.85772861000027
#define c_pi_13_64   0.80320753148064
#define c_pi_15_64   0.74095112535496
#define c_pi_17_64   0.67155895484702
#define c_pi_19_64   0.59569930449243
#define c_pi_21_64   0.51410274419322
#define c_pi_23_64   0.42755509343028
#define c_pi_25_64   0.33688985339222
#define c_pi_27_64   0.24298017990326
#define c_pi_29_64   0.14673047445536
#define c_pi_31_64   0.04906767432742

#define c_pi_1_32    0.99518472667220    //表示cos(PI*(奇数)/32)的值,以下类似
#define c_pi_3_32    0.95694033573221
#define c_pi_5_32    0.88192126434836
#define c_pi_7_32    0.77301045336274
#define c_pi_9_32    0.63439328416365
#define c_pi_11_32   0.47139673682600
#define c_pi_13_32   0.29028467725446
#define c_pi_15_32   0.09801714032956

#define c_pi_1_16    0.98078528040323     //表示cos(PI*(奇数)/16)的值,以下类似
#define c_pi_3_16    0.83146961230255
#define c_pi_5_16    0.55557023301960
#define c_pi_7_16    0.19509032201613

#define c_pi_1_8     0.92387953251129	  //表示cos(PI*(奇数)/8)的值,以下类似
#define c_pi_3_8	 0.38268343236509

#define c_pi_1_4     0.70710678118655      //表示cos(PI*(奇数)/4)的值,以下类似


double * DCT_64(double * f)
{
   /*double f[]={1,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0, 
		        0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0,
				0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0,
				0,0,0,0, 0,0,0,0, 0,0,0,0, 0,0,0,0};
   */
	double * F,* t0,* t1,* t2,* t3,* t4,* t5,* t6,* t7;

	//01    f[64]乘以A64分解式的右侧 
	//     64次加法     0次乘法
	{
		t0=(double *)malloc(512);

		t0[0]=f[0]+f[63];		t0[32]=f[0]-f[63];
		t0[1]=f[1]+f[62];		t0[33]=f[1]-f[62];
		t0[2]=f[2]+f[61];		t0[34]=f[2]-f[61];
		t0[3]=f[3]+f[60];		t0[35]=f[3]-f[60];
		t0[4]=f[4]+f[59];		t0[36]=f[4]-f[59];
		t0[5]=f[5]+f[58];		t0[37]=f[5]-f[58];
		t0[6]=f[6]+f[57];		t0[38]=f[6]-f[57];
		t0[7]=f[7]+f[56];		t0[39]=f[7]-f[56];
		t0[8]=f[8]+f[55];		t0[40]=f[8]-f[55];
		t0[9]=f[9]+f[54];		t0[41]=f[9]-f[54];
		t0[10]=f[10]+f[53];		t0[42]=f[10]-f[53];
		t0[11]=f[11]+f[52];		t0[43]=f[11]-f[51];
		t0[12]=f[12]+f[51];		t0[44]=f[12]-f[51];
		t0[13]=f[13]+f[50];		t0[45]=f[13]-f[50];
		t0[14]=f[14]+f[49];		t0[46]=f[14]-f[49];
		t0[15]=f[15]+f[48];		t0[47]=f[15]-f[48];
		t0[16]=f[16]+f[47];		t0[48]=f[16]-f[47];
		t0[17]=f[17]+f[46];		t0[49]=f[17]-f[46];
		t0[18]=f[18]+f[45];		t0[50]=f[18]-f[45];
		t0[19]=f[19]+f[44];		t0[51]=f[19]-f[44];
		t0[20]=f[20]+f[43];		t0[52]=f[20]-f[43];
		t0[21]=f[21]+f[42];		t0[53]=f[21]-f[42];
		t0[22]=f[22]+f[41];		t0[54]=f[22]-f[41];
		t0[23]=f[23]+f[40];		t0[55]=f[23]-f[40];
		t0[24]=f[24]+f[39];		t0[56]=f[24]-f[39];
		t0[25]=f[25]+f[38];		t0[57]=f[25]-f[38];
		t0[26]=f[26]+f[37];		t0[58]=f[26]-f[37];
		t0[27]=f[27]+f[36];		t0[59]=f[27]-f[36];
		t0[28]=f[28]+f[35];		t0[60]=f[28]-f[35];
		t0[29]=f[29]+f[34];		t0[61]=f[29]-f[34];
		t0[30]=f[30]+f[33];		t0[62]=f[30]-f[33];
		t0[31]=f[31]+f[32];		t0[63]=f[31]-f[32];
	}


	//02   t0[0]--t0[31]乘以A32
	{
		//02_1  t0[0]--t0[31]乘以A32分解式右侧的矩阵
		//   32次加法  0 次乘法
		{
			t1=(double *)malloc(256);   //为t1  分配32 个double的空间

			t1[0]=t0[0]+t0[31];      t1[16]=t0[0]-t0[31];
			t1[1]=t0[1]+t0[30];      t1[17]=t0[1]-t0[30];
			t1[2]=t0[2]+t0[29];      t1[18]=t0[2]-t0[29];
			t1[3]=t0[3]+t0[28];      t1[19]=t0[3]-t0[28];
			t1[4]=t0[4]+t0[27];      t1[20]=t0[4]-t0[27];
			t1[5]=t0[5]+t0[26];      t1[21]=t0[5]-t0[26];
			t1[6]=t0[6]+t0[25];      t1[22]=t0[6]-t0[25];
			t1[7]=t0[7]+t0[24];      t1[23]=t0[7]-t0[24];
			t1[8]=t0[8]+t0[23];      t1[24]=t0[8]-t0[23];
			t1[9]=t0[9]+t0[22];      t1[25]=t0[9]-t0[22];
			t1[10]=t0[10]+t0[21];    t1[26]=t0[10]-t0[21];
			t1[11]=t0[11]+t0[20];    t1[27]=t0[11]-t0[20];
			t1[12]=t0[12]+t0[19];    t1[28]=t0[12]-t0[19];
			t1[13]=t0[13]+t0[18];    t1[29]=t0[13]-t0[18];
			t1[14]=t0[14]+t0[17];    t1[30]=t0[14]-t0[17];
			t1[15]=t0[15]+t0[16];    t1[31]=t0[15]-t0[16];
			
		}

		
		


		//02_2  t1[0]--t1[15]乘以A16    记录在f[0]--f[15]里
		//
		{

			//02_2_1   t1[16]乘以A16分解式右侧的矩阵
			//      16次加法     0 次乘法
			{
				t2=(double *)malloc(128);     //为t2 分配16 个double的空间
				
				t2[0]=t1[0]+t1[15];		t2[8]=t1[0]-t1[15];
				t2[1]=t1[1]+t1[14];		t2[9]=t1[1]-t1[14];
				t2[2]=t1[2]+t1[13];		t2[10]=t1[2]-t1[13];
				t2[3]=t1[3]+t1[12];		t2[11]=t1[3]-t1[12];
				t2[4]=t1[4]+t1[11];		t2[12]=t1[4]-t1[11];
				t2[5]=t1[5]+t1[10];		t2[13]=t1[5]-t1[10];
				t2[6]=t1[6]+t1[9];		t2[14]=t1[6]-t1[9];
				t2[7]=t1[7]+t1[8];		t2[15]=t1[7]-t1[8];
			}



			//02_2_2   t2[0]--t2[7]乘以A8       最后仍保存到t2[0]--t2[7]里
			//
			{
				//02_2_2_1   t2[0]--t2[7]乘以A8后分解式右侧的矩阵
				//        8次加法   0次乘法
				{
					t3=(double *)malloc(64);       //为t3 分配8 个double的空间

					t3[0]=t2[0]+t2[7];		t3[4]=t2[0]-t2[7];
					t3[1]=t2[1]+t2[6];		t3[5]=t2[1]-t2[6];
					t3[2]=t2[2]+t2[5];		t3[6]=t2[2]-t2[5];
					t3[3]=t2[3]+t2[4];		t3[7]=t2[3]-t2[4];
				}

				//02_2_2_2      t3[0]--t3[3]乘以A4
				//
				{
					//02_2_2_2_1       t3[0]--t3[3]乘以A4分解式右侧的矩阵
					//				4次加法     0次乘法
					{
						t4=(double *)malloc(32);		//为t4 分配4 个double的空间

						t4[0]=t3[0]+t3[3];		t4[2]=t3[0]-t3[3];
						t4[1]=t3[1]+t3[2];		t4[3]=t3[1]-t3[2];
					}
					//02_2_2_2_2      t4[0]--t4[1]乘以A2    这时可以直接乘了!!!!
					//             2 次加法    1 次乘法
					{
						t5=(double *)malloc(32);		//为t5 分配4 个double的空间
						t5[0]=t4[0]+t4[1];        
						t5[1]=(t4[0]-t4[1])*c_pi_1_4;
					}
					//02_2_2_2_3       t4[2]--t4[3]乘以B2
					//				 4次加法   3次乘法
					{
						t4[2]*=c_pi_1_8;		t4[3]*=c_pi_3_8;

						t5[2]=t4[2]+t4[3];		t5[3]=(t4[2]-t4[3])*c_pi_1_4;

						t5[3]=t5[3]+t5[3];

						t5[3]=t5[3]-t5[2];
					}//   此时 t5 记录的是P4之前的列向量
					//02_2_2_2_4    t5乘以P4  用t4记录
					{
						t4[0]=t5[0];
						t4[1]=t5[2];
						t4[2]=t5[1];
						t4[3]=t5[3];
					}

				}//此时 t4[0]--t4[3]记录的是乘以A4后的t3[0]--t3[3]元素  t5 已没用

				//02_2_2_3        t3[4]--t3[7]乘以B4
				//        
				{
					//02_2_2_3_1	  t3[4]--t3[7]乘以B4分解式最右侧的以 cos(PI*(2*i+1)/2/8)为对角元的对角阵
					//            0次加法   4次乘法
					{
						t3[4]*=c_pi_1_16;		t3[5]*=c_pi_3_16;
						t3[6]*=c_pi_5_16;		t3[7]*=c_pi_7_16;
					}

					//02_2_2_3_2   接下来 t3[4]--t3[7]乘以A4
					{
						//02_2_2_3_2_1         t3[4]--t3[7]乘以A4分解式右侧的矩阵 t5记录
						//				    4次加法     0次乘法
						{
							t5[0]=t3[4]+t3[7];		t5[2]=t3[4]-t3[7];
							t5[1]=t3[5]+t3[6];		t5[3]=t3[5]-t3[6];
						}
						//02_2_2_3_2_2    		t5[0]--t5[1]乘以乘以A2    这时可以直接乘了!!!!
						//                  2 次加法    1 次乘法
						{
							t6=(double *)malloc(32);		//为t6 分配4 个double的空间
							t6[0]=t5[0]+t5[1];        
							t6[1]=(t5[0]-t5[1])*c_pi_1_4;
						}
						//02_2_2_3_2_3       t5[2]--t5[3]乘以B2
						//				  4次加法   3次乘法
						{
							t5[2]*=c_pi_1_8;		t5[3]*=c_pi_3_8;

							t6[2]=t5[2]+t5[3];		t6[3]=(t5[2]-t5[3])*c_pi_1_4;

							t6[3]=t6[3]+t6[3];

							t6[3]=t6[3]-t6[2];
						}//   此时 t6 记录的是P4之前的列向量
						//02_2_2_3_2_4    t6乘以P4   用t5记录
						//
						{
							t5[0]=t6[0];
							t5[1]=t6[2];
							t5[2]=t6[1];
							t5[3]=t6[3];
						}//此时 t5[0]--t5[3]记录的是乘以A4后的t3[4]--t3[7]元素  t6 已没用
					}

					//02_2_2_3_3         t5[0]--t5[3]乘以 一二阵 和 正负一阵
					//                9次加法   0次乘法
					{
							t5[1]=t5[1]+t5[1];  t5[2]=t5[2]+t5[2];  t5[3]=t5[3]+t5[3];

							t6[0]=t5[0];
							t6[1]=t5[1]-t5[0];
							t6[2]=t5[0]-t5[1]+t5[2];
							t6[3]=t5[1]-t5[0]+t5[3]-t5[2];

					}
				}// 此时t6[0]--t6[3]记录的是t3[4]--t3[7]乘以B4后的元素
				//02_2_2_4         P8之前的t3[0]--t3[7]乘以P8, 仍记录在t2[0]--t2[7]里
				//
				{
					t2[0]=t4[0];
					t2[1]=t6[0];
					t2[2]=t4[1];
					t2[3]=t6[1];
					t2[4]=t4[2];
					t2[5]=t6[2];
					t2[6]=t4[3];
					t2[7]=t6[3];
				}
			}//此时 t2[0]--t2[7]记录的是乘完A8的t2[0]--t2[7]



			//02_2_3    t2[8]--t2[15]乘以B8       仍保存到t2[8]--t2[15]里
			{
				//02_2_3_1  t2[8]--t2[15]乘以B8分解式最右侧的以 cos(PI*(2*i+1)/2/16)为对角元的对角阵
				//        0次加法    8次乘法
				{
					t2[8]*=c_pi_1_32;		t2[9]*=c_pi_3_32;
					t2[10]*=c_pi_5_32;		t2[11]*=c_pi_7_32;
					t2[12]*=c_pi_9_32;		t2[13]*=c_pi_11_32;
					t2[14]*=c_pi_13_32;		t2[15]*=c_pi_15_32;
				}

				//02_2_3_2    接下来 t2[8]--t2[15]乘以A8
				{
					//02_2_3_2_1   tt2[8]--t2[15]乘以A8后分解式右侧的矩阵
					//          8次加法   0次乘法
					{
						free(t4);
						t4=(double *)malloc(64);       //重新为t4 分配8 个double的空间

						t4[0]=t2[8]+t2[15];		t4[4]=t2[8]-t2[15];
						t4[1]=t2[9]+t2[14];		t4[5]=t2[9]-t2[14];
						t4[2]=t2[10]+t2[13];	t4[6]=t2[10]-t2[13];
						t4[3]=t2[11]+t2[12];	t4[7]=t2[11]-t2[12];
					}
					//02_2_3_2_2      t4[0]--t4[3]乘以A4
					{
						//02_2_3_2_2_1    t4[0]--t4[3]乘以A4分解式右侧的矩阵
					   //				4次加法     0次乘法
						{
							t5[0]=t4[0]+t4[3];		t5[2]=t4[0]-t4[3];
							t5[1]=t4[1]+t4[2];		t5[3]=t4[1]-t4[2];

						}
						//02_2_3_2_2_2    t5[0]--t5[1]乘以A2    这时可以直接乘了!!!!
						//             2 次加法    1 次乘法
						{
							t6[0]=t5[0]+t5[1];
							t6[1]=(t5[0]-t5[1])*c_pi_1_4;
						}

						//02_2_3_2_2_3      t5[2]--t5[3]乘以B2
						//				 4次加法   3次乘法
						{
							t5[2]*=c_pi_1_8;		t5[3]*=c_pi_3_8;

							t6[2]=t5[2]+t5[3];		t6[3]=(t5[2]-t5[3])*c_pi_1_4;

							t6[3]=t6[3]+t6[3];

							t6[3]=t6[3]-t6[2];
						}//   此时 t6 记录的是P4之前的列向量

						//02_2_3_2_2_4     t6乘以P4  用t5记录
						{
							t5[0]=t6[0];
							t5[1]=t6[2];
							t5[2]=t6[1];
							t5[3]=t6[3];
						}

					}//此时 t5[0]--t5[3]记录的是乘以A4后的t4[0]--t4[3]元素  t6 已没用

					//02_2_3_2_3		t4[4]--t4[7]乘以B4
					//
					{
						//02_2_3_2_3_1        t4[4]--t4[7]乘以B4分解式最右侧的以 cos(PI*(2*i+1)/2/8)为对角元的对角阵
						//                 0次加法   4次乘法
						{
							t4[4]*=c_pi_1_16;		t4[5]*=c_pi_3_16;
							t4[6]*=c_pi_5_16;		t4[7]*=c_pi_7_16;
						}
						//02_2_3_2_3_2     接下来 t4[4]--t4[7]乘以A4
						{
							//02_2_3_2_3_2_1       t4[4]--t4[7]乘以A4分解式右侧的矩阵 t6记录
							//				    4次加法     0次乘法
							{
								t6[0]=t4[4]+t4[7];		t6[2]=t4[4]-t4[7];
								t6[1]=t4[5]+t4[6];		t6[3]=t4[5]-t4[6];

							}
							//02_2_3_2_3_2_2      t6[0]--t6[1]乘以乘以A2    这时可以直接乘了!!!!
							//                 2 次加法    1 次乘法
							{
								t7=(double *)malloc(32);		//为t7 分配4 个double的空间
								t7[0]=t6[0]+t6[1];        
								t7[1]=(t6[0]-t6[1])*c_pi_1_4;

							}
							//02_2_3_2_3_2_3     t6[2]--t6[3]乘以B2
							//				  4次加法   3次乘法
							{
								t6[2]*=c_pi_1_8;		t6[3]*=c_pi_3_8;

								t7[2]=t6[2]+t6[3];		t7[3]=(t6[2]-t6[3])*c_pi_1_4;

								t7[3]=t7[3]+t7[3];

								t7[3]=t7[3]-t7[2];

							}//  此时    t7 记录的是P4之前的列向量
							//02_2_3_2_3_2_4          t7乘以P4   用t6记录
							//
							{
								t6[0]=t7[0];
								t6[1]=t7[2];
								t6[2]=t7[1];
								t6[3]=t7[3];
							}//此时 t6[0]--t6[3]记录的是乘以A4后的t4[4]--t4[7]元素  t7 已没用
						}
						//02_2_3_2_3_3      t6[0]--t6[3]乘以 一二阵 和 正负一阵
						//                9次加法   0次乘法
						{
							t6[1]=t6[1]+t6[1];	t6[2]=t6[2]+t6[2];	t6[3]=t6[3]+t6[3];

							
							t7[0]=t6[0];
							t7[1]=t6[1]-t6[0];
							t7[2]=t6[0]-t6[1]+t6[2];
							t7[3]=t6[1]-t6[0]+t6[3]-t6[2];

						}// 此时t7[0]--t7[3]记录的是t4[4]--t4[7]乘以B4后的元素
					}
					//02_2_3_2_4      P8之前的43[0]--t4[7]乘以P8,记录在t4[0]--t4[7]里
					{
						t4[0]=t5[0];
						t4[1]=t7[0];
						t4[2]=t5[1];
						t4[3]=t7[1];
						t4[4]=t5[2];
						t4[5]=t7[2];
						t4[6]=t5[3];
						t4[7]=t7[3];
					}
				}//此时t4[0]--t4[7]记录的是乘完A8后的t2[8]--t2[15]

				//02_2_3_3    t4[0]--t4[7]乘以 一二阵 和 正负一阵,得到  ,记录在t2[8]--t2[15]里
				{
					t4[1]=t4[1]+t4[1];

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -