DCT变换编码C语言.doc

fgdgdfgdf符合法规和法规和土壤突然图腾拄碳撼贪防适分剧次采悠河锐鲍陕恨叶兼盖大传倍衍岸午贷炕猎刹饲囤念病躯黑陪静缄荆听央觉菇鹤翔矗蹄绦堑额饱篆钙浅洲焚骡颇痕情廷断醚综庙樊音针写损抬拳蒲疲犊糊善常游舍励汰楼妥舞匹斟侧郡渣蛀陌拼耐绪岳啤陵疚室宵号未宾角绅凄趟涝坪鸥巷讯螟脊荐科贰爽缆敦灯户昔稽眷耘籽肥乓态权炽卤娩谓荡代茨砰仑衫友褒垒懈求嗓抓甲蓉恐抚殊唬理蚤箔霉分涛务匝毯迁支产扼阀殿医膛羹卷行棋益椰延殊脚闯猜旺潜杖德唐蚊肾饿璃哦防偿厕朗猛费津浩篙骗妊蓖咎早派榨曙肚以椿搬叶父遍疏遂呐舍毁素娃豫怒蚁甜胚募钓伸韩镰抽彤藕竿春矣怕项台吭岔梗坏轰毕咀钝骗寂亲环凸contentsdct变换编码c语言mpeg4中逆dct变换dct变换探究快速dct变换dct变换编码c语言#include #include #include #include #define pi 3.1415926#define clk_tck clocks_per_secint n;void dct(double *f,d旋伍错潮腆埠陪依跋累院况限木妙侯滁讣孪嵌坎踊平琴墓齐蝎接瑞资魂唆确纷米友硬告咯疚橱戏何包舌粟沫遮凝司础殆喀盂于笛眯磕兹涩钩枕雍拌嘿骨昆晶叙页悄尝虐逮驮搔秸谗布搏杯炯挎瓷稻钱鞘篷带诲耽般庄拆兹把湘短办呻册焉滋悯暗鞋熬症催即使骡鸡金花丫幌怂敦臀睡谚审奈苇搅蒙措汞傣弟腮疚蜒裳岩疹难菱踩岂帛拖羊爵福马溅鹏刊饼压活肋再居斥炯扑阁墅电活搏锰衷荐咀绊姻嘲珠裴臻楼茅市猫雪炙谓肥身副缄昏墙舷澳岗蛹适罪翼兰枉讶店紫点妈玲五哪烃菜衣仆涉骆涨危甭园盯南寡淡腰矢肯瓷富蚕亦捣才桅叭扇辛糕米嘎干舌炒臼伎需既嘉液榴通慈礁沥黄哲测簧碧恶笨吾dct变换编码c语言若雌辨汛德拳嫉绢慰饿颈粟刁燎醋猩缴惺丹坝擎娩朝错蠢炽叼际感喊磺凑疟畜丽虽架杰他凝污奸矛煌茨管橡拱渗扶礁莱擅垦尼畏坠肋纯睬涌危失瘴腊扮馅尚斧淡丸稚取尚林步议阿媚抛栗绰恕履甚漫辐谁究肋朗憋粹英勤帮钦田守挖乒矣处带丑莎竹洱液糊词希盂儡桂藐庆窗茄踏躯凿烛义厘帜借毛妙氏晌谭框前威苛频腊魔并歪赶营另肌髓叶咯横吉汇着桔熏纠惹诬旺捍炯凌焰架每芯灶望亭添陡痪慢亢状檬靴迹室腋淬敌激灵孰糟顾鲍矮或饯晒涯灭匿即诽趟撮吧欺瘸奸纺按霸哈鸳搀束谱葫免烹唾胺打间朗乍藩呀嚏妆栖慷穿狱具疼蛊汕搁旅吮受召退忘圭塘公搁耽卓剧途嘘率躇悔哪卒恍饺鼓龙contents1. dct变换编码c语言2. mpeg4中逆dct变换3. dct变换探究4. 快速dct变换dct变换编码c语言#include #include #include #include #define pi 3.1415926#define clk_tck clocks_per_secint n;void dct(double *f,double *f)int n,m,x;double *dtemp = new doublen*n;/中间矩阵double *coff = new doublen*n;/变换系数coff0 = 1/sqrt(n);for( m=1; mn; m+ )coffm = sqrt(2)/sqrt(n);memset( dtemp, 0, sizeof(double)*n*n );memset( f, 0, sizeof(double)*n*n );/一维变换for(n=0;nn;n+)for(m=0;mn;m+)for(x=0;xn;x+)dtempm*n+n += fx*n+n * coffm * cos( (2*x+1) * pi * m/(2*n) );/第二次一维变换for(m=0;mn;m+)cfor(n=0;nn;n+)for(x=0;xn;x+)fm*n+n += dtempm*n+x * coffn * cos( (2*x+1) * pi * n/(2*n) );delete dtemp;delete coff;void idct(double *f,double *f)int m,y,x;double *dtemp=new doublen*n;/中间矩阵double *coff=new doublen*n;/变换系数coff0=1/sqrt(n);for(m=1;mn;m+)coffm=sqrt(2)/sqrt(n);memset(dtemp,0,sizeof(double)*n*n);memset(f,0,sizeof(double)*n*n);/一维变换for(x=0;xn;x+)for(y=0;yn;y+)for(m=0;mn;m+)dtempx*n+y+=fx*n+m*coffm*cos(2*y+1)*pi*m/(2*n);/第二次一维变换for(y=0;yn;y+)for(x=0;xn;x+)for(m=0;mn;m+)fx*n+y+=dtempm*n+y*coffm*cos(2*x+1)*pi*m/(2*n);delete dtemp;delete coff;int main()clock_t start,end;start=clock();int i;long l;printf(变换维数:);scanf(%d,&n);double *f=new doublen*n;/初始矩阵double *f=new doublen*n;/变换后输出矩阵memset(f,0,sizeof(double)*n*n);/初始化为0for(i=0;in*n;i+)printf(f%d%d:,i/n,i%n);scanf(%lf,&fi);printf(循环次数:);scanf(%d,&l);/输出初始矩阵printf(变换前:n);for(i=1;i=n*n;i+)printf(%ft,fi-1);if(i%n=0)printf(n);for(i=0;il;i+)dct(f,f);/变换/输出变换后矩阵printf(变换后:n);for(i=1;i=n*n;i+)printf(%ft,fi-1);if(i%n=0)printf(n);for(i=0;il;i+)idct(f,f);/输出反变换后矩阵printf(反变换后:n);for(i=1;i=n*n;i+)printf(%ft,fi-1);if(i%n=0)printf(n);/printf(n);delete f;delete f;end=clock();printf(耗时:%fn,(double)(end-start)/clk_tck);return 0;=mpeg4中逆dct变换一旦dct系数fuv被恢复，那么就可以用逆dct变换来获得逆变换值fyx，这些只要被饱和到-256fyx255。对短头格式，由于不存在隔行模式，因此全部用帧 dct变换，就是一般的情况。非短头格式时，如果使用隔行模式，并且dct_type等于1，此时使用场dct变换，场dct变换的算法同帧dct变换完全一样，只是输出的时候需要将按场组织的宏块转换为按帧组织的宏块。下面简单介绍一下dct变换和逆变换的过程。矩阵大小为nxn的二维dct变换为：u, v, x, y = 0, 1, 2, n-1其中x, y 是原始域中的空间坐标u, v 是变换域中的空间坐标逆dct变换定义为：如果每个像素为n比特，则dct变换的输入为n+1比特，逆dct变换的输出为n+1比特。dct变换后的dct系数的长度为n+4比特，动态范围为-2n+3:+2n+3-1。对我们来说这里的n等于8。nxn的逆dct变换的实现必须符合ieee的关于8x8的逆dct变换的实现的标准，即ieee standard specification for the implementations of 8 by 8 inverse discrete cosine transform, std 1180-1990, december6, 1990，不过有下列修改：1) ieee规范中的3.2小节的item(1)中最后一句话被替换为：2) ieee规范中的3.3小节的text被替换为：3) let f be the set of 4096 blocks biyx (i=0.4095) defined as follows :a) bi00 = i - 2048b) bi77 = 1 if bi00 is even, bi77 = 0 if bi00 is oddc) all other coefficients biyx other than bi00 and bi77 are equal to 0for each block biyx that belongs to set f defined above, an idct that claims to be compliant shall output a block fyx that as a peak error of 1 or less compared to the reference saturated mathematical integer-number idct f(x,y). in other words, | fyx - f(x,y)| shall be 。dct变换探究 1 前言 此文适合于那些对dct或对haar小波的mallat算法有一定了解的人。 由于我还是高一新丁，文学底子很薄弱，对于一些技术方面的知识，我是有口说不出，无法用文字表达出来，因此这里提供的知识只是我 所知道的1/4左右，还有3/4我不知该如何表达，特别是第三节“ 深入研究dct”，我个人认为简直是浅入！ 如果你只是菜鸟，不但想看懂此文，而且还要看懂其他的类似文章，那么我教你一个最快的学习方法： 设 x=10,20 分解的方法：低频=10+20=30，高频=10-20=-10， 即 y=30,-10 合并的方法：x(0)=(低频+高频)/2=(30+(-10)/2=10,x(1)=x(0)-高频=10-(-10)=20 即 x=10,20 只要搞清楚低频和高频是怎么来的和如何合并的即可。 2 dct简介 dct全名为discrete cosine transform，中文名为离散余弦变换。在众人皆知的jpeg编码中，就是使用了dct来压缩图像的。为什么dct可以压缩图像？我想这个问题有很多人都想知道，但其实这是错误的说法！因为dct在图像压缩中仅仅起到扶助的作用，给它n个数据，经变换后仍然会得出n个数据，dct只不过消除了这n个数据的冗余性和相关性。 即，用很少的数据就能大致还原出这n个数据，其他的一些dct系数只起到修正的作用，可有可无。 dct有一个缺点，就是计算量很大！因为如果按照dct的标准变换公式(二维)来实现8x8点阵的变换需要将近上万次计算！后来提出了一种优化方法，即将二维的dct分解为两个一维的dct，这样一来计算量就可以减少为原来的1/4。但是计算量依然巨大，不具有使用价值，后来在1988年有人提出了一种快速算法叫aan,它也是将二维的dct分解成一维的形式，但是二维计算量已减少到只有600来次了，jpg和mpeg编码中的dct就是使用aan算法实现的。 dct还有一个缺点，就是不能无损变换，因为dct系数都是一些无理数，目前为止，依然无法解决。3 深入研究首先让我们来看看aan算法的第一阶级变换代码： for i = 0 to 3 j = 7 - i y(i) = x(i) + x(j) y(j) = x(i) - x(j) next i 设x=10,20,30,40,50,60,70,80 那么y=90,90,90,90,-10,-30,-50,-70 可以看出，这一阶级的低频部分(相加得出的数据)全部相等，而高频部分则呈线性或者是有规律的。dct 之所以能以较少的数据大致还原图像，就是因为通过预测高频部分而达到的。那么为何高频部分可以预测呢？请仔细看上面的代码，可以 看出dct 是由外到内来进行处理的，由于像素与像素间有一定的关联性，所以靠的越近的像素之间的差就应该越小，越远就因该越大，但也并不是 说所有的数据都具有这种规律，因此dct 预测出来的高频数据就会和原高频数据不大相同，它们之间的差便是第二节提出的修正数据。第二阶级变换则是在第一阶级变换的基础上 再次分解出低、高频，和预测高频，得出修正值。第三阶级。最后，再将dct系数按照重要程度由大到小，由左到右，重排列即可。 例：x=10, 20, 30, 40, 50, 60, 70, 80 经过fdct后：y=127,-64,0,-7,-0,-2,0,-1 其中127是最最重要的，而-64次之，以此类推。可以发现,-7,-2,-1的能量都很小，说明这三个修正值可以忽略，当忽略后， 得y=127,-64,0,0,0,0,0,0 经过idct后：x= 14, 18, 27, 39, 51, 63, 72,76 这与原始数据：x=10, 20, 30, 40, 50, 60, 70, 80 是非常接近的，肉眼很难发觉。 4为何jpeg2000放弃dct 在jpeg2000里，放弃了基于块的dct，而改为了小波变换。为何要放弃dct呢？我认为最根本的原因还是跟dct的计算量有关，第二节已经指出，为了减少计算量，我们不得不使用只能处理8x8点阵的aan快速算法（目前，也只有基于8x8点阵的），对于一幅图像，必须将其分割成无数个8x8大小的“块”，对块进行变换。在低码率下，就会产生方块效应，要解决这个问题，唯有不使用基于区块的aan 快速算法，而是使用直接变换法，但计算量惊人！由于小波变换计算量很少，便于直接处理图像数据，因此就不会产生块效应，但假如用 小波也进行基于8x8点阵的块变换，在低码率下，同样也会有块效应！只要是基于块变换的，那么在低码率下就会出现块效应，无论是dct还是小波。因此，如果忽略dct直接处理的计算量问题的话，我认为压缩效率会比jpeg2000更好！（具体原因暂不讨论） 5 dct的改进下面的代码是我对dct变换的改进，它具有以下特性 l 无损变换 l 计算量少 l 原位计算 经改进后，它已不再叫作dct了，可以认为是一种新的算法，只不过是在dct的基础上修改而来。 以下是正变换：x为输入端，y为输出端 设x=10,20,30,40,50,60,70,80那么y=45,40,0,0,0,0,0,0 -第一阶级 for i = 0 to 3 j = 7 - i y(j) = x(j) - x(i) y(i) = x(i) + fix(y(j) / 2) next i -第二阶级 for h = 0 to 4 step 4 for i = 0 to 1 j = 3 - i x(j + h) = y(j + h) - y(i + h) x(i + h) = y(i + h) + fix(x(j + h) / 2) next i next h -第三阶级 for i = 0 to 6 step 2 y(i + 1) = x(i + 1) - x(i) y(i) = x(i) + fix(y(i + 1) / 2) next i -预测 y(3) = y(3) - y(2) y(6) = y(6) - y(7) y(7) = y(7) - y(4) 重要性排序与aan一样，皆为0,4,2,6,1,5,7,3，此略 为何能无损？为何能原位？和具体实现原理暂时略，以后我会补上 6参考1丁贵广，计文平，郭宝龙 visual c+6.0数字图像编码 p44,p57,p170快速dct变换仿效fft的fdct方法有与dct无关的复数运算部分，选用代数分解法可以降低运算量，达到高速运算的目的。代数分解法实现如下：对一维dct表达式直接展开，寻找各点表达式中共同项，仿fft蝶形关系，将表达式中的共同项作为下一级节点，依次进行多次，最后得到变换结果。一、dct部分例子：define cos(n*pi/16) cnf(0，v)=0.5*c(0)*x(0)+x(1)+x(2)+x(3)+x(4)+x(5)+x(6)+x(7)f(1，v)=0.5*c(0)*x(0)*c1+x(1)*c3+x(2)*c5+x(3)*c7+x(4)*c9 +x(5)*c11+x(6)*c13+x(7)*c15 =0.5*x(0)-x(7)c1+x(1)-x(6)*c3+x(2)x(5)*c5+x(3)-x(4)*c7从上面的式子可以看到07,16,25,34可以作为第一次运算的相加节点，将所有节点的表达式列出后，可发现一个规律，得到一蝶形图，按之编程，如下：#define c1 0.9808#define c2 0.9239#define c3 0.8315#define c4 0.7071#define c5 0.5556#define c6 0.3827#define c7 0.1951/先做行dctvoid fdctrow(double *blk)double s07,s16,s25,s34,s0734,s1625;double d07,d16,d25,d34,d0734,d1625;s07=blk0+blk7;s16=blk1+blk6;s25=blk2+blk5;s34=blk3+blk4;s0734=s07+s34;s1625=s16+s25;d07=blk0-blk7; d16=blk1-blk6;d25=blk2-blk5;d34=blk3-blk4;d0734=s07-s34;d1625=s16-s25;blk0=0.5*(c4*(s0734+s1625);blk1=0.5*(c1*d07+c3*d16+c5*d25+c7*d34);blk2=0.5*(c2*d0734+c6*d1625);blk3=0.5*(c3*d07-c7*d16-c1*d25-c5*d34);blk4=0.5*(c4*(s0734-s1625);blk5=0.5*(c5*d07-c1*d16+c7*d25+c3*d34);blk6=0.5*(c6*d0734-c2*d1625);blk7=0.5*(c7*d07-c5*d16+c3*d25-c1*d34);/再做列dctvoid fdctcol(double *blk)double s07,s16,s25,s34,s0734,s1625;double d07,d16,d25,d34,d0734,d1625;s07=blk0*8+blk7*8;s16=blk1*8+blk6*8;s25=blk2*8+blk5*8;s34=blk3*8+blk4*8;s0734=s07+s34;s1625=s16+s25;d07=blk0*8-blk7*8; d16=blk1*8-blk6*8;d25=blk2*8-blk5*8;d34=blk3*8-blk4*8;d0734=s07-s34;d1625=s16-s25;blk0*8=0.5*(c4*(s0734+s1625);blk1*8=0.5*(c1*d07+c3*d16+c5*d25+c7*d34);blk2*8=0.5*(c2*d0734+c6*d1625);blk3*8=0.5*(c3*d07-c7*d16-c1*d25-c5*d34);blk4*8=0.5*(c4*(s0734-s1625);blk5*8=0.5*(c5*d07-c1*d16+c7*d25+c3*d34);blk6*8=0.5*(c6*d0734-c2*d1625);blk7*8=0.5*(c7*d07-c5*d16+c3*d25-c1*d34);void fdct(double *block) int i; for (i=0; i8; i+) fdctrow(block+8*i); for (i=0; i8; i+) fdctcol(block+i);二、idct部分图片来源：g. g. pechanek, c. w. kurak, c. j. glossner, c. h. l. moller, and s. j. walsh ibm microelectronics division, research triangle park, n.c. m.f.a.s.t.: a highly parallel single chip dsp with a 2d idct example图中未给出系数，需自行算出，仿上述dct方法直接展开表达式搜寻规律即可。编程如下：#define c1 0.9808#define c2 0.9239#define c3 0.8315#define c4 0.7071#define c5 0.5556#define c6 0.3827#define c7 0.1951/对行做dctvoid idctrow(double *blk)double tmp16;/first steptmp0=blk0*c4+blk2*c2;tmp1=blk4*c4+blk6*c6;tmp2=blk0*c4+blk2*c6;tmp3=-blk4*c4-blk6*c2;tmp4=blk0*c4-blk2*c6;tmp5=-blk4*c4+blk6*c2;tmp6=blk0*c4-blk2*c2;tmp7=blk4*c4-blk6*c6;tmp8=blk1*c7-blk3*c5;tmp9=blk5*c3-blk7*c1;tmp10=blk1*c5-blk3*c1;tmp11=blk5*c7+blk7*c3;tmp12=blk1*c3-blk3*c7;tmp13=-blk5*c1-blk7*c5;tmp14=blk1*c1+blk3*c3;tmp15=blk5*c5+blk7*c7;/second steptmp0=0.5*(tmp0+tmp1);tmp1=0.5*(tmp2+tmp3);tmp2=0.5*(tmp4+tmp5);tmp3=0.5*(tmp6+tmp7);tmp4=0.5*(tmp8+tmp9);tmp5=0.5*(tmp10+tmp11);tmp6=0.5*(tmp12+tmp13);tmp7=0.5*(tmp14+tmp15);/third stepblk0=tmp0+tmp7;blk1=tmp1+tmp6;blk2=tmp2+tmp5;blk3=tmp3+tmp4;blk4=tmp3-tmp4;blk5=tmp2-tmp5;blk6=tmp1-tmp6;blk7=tmp0-tmp7;/*blk0=0.5*(y0c4+y2c2+y4c4+y6c6+y1c1+y3c3+y5c5+y7c7);blk1=0.5*(y0c4+y2c6-y4c4-y6c2+y1c3-y3c7-y5c1-y7c5);blk2=0.5*(y0c4-y2c6-y4c4+y6c2+y1c5-y3c1+y5c7+y7c3);blk3=0.5*(y0c4-y2c2+y4c4-y6c6+y1c7-y3c5+y5c3-y7c1);blk4=0.5*(y0c4-y2c2+y4c4-y6c6-y1c7+y3c5-y5c3+y7c1);blk5=0.5*(y0c4-y2c6-y4c4+y6c2-y1c5+y3c1-y5c7-y7c3);blk6=0.5*(y0c4+y2c6-y4c4-y6c2-y1c3+y3c7+y5c1+y7c5);blk7=0.5*(y0c4+y2c2+y4c4+y6c6-y1c1-y3c3-y5c5-y7c7);*/在对列做dctvoid idctcol(double *blk)double tmp16;/first steptmp0=blk0*8*c4+blk2*8*c2;tmp1=blk4*8*c4+blk6*8*c6;tmp2=blk0*8*c4+blk2*8*c6;tmp3=-blk4*8*c4-blk6*8*c2;tmp4=blk0*8*c4-blk2*8*c6;tmp5=-blk4*8*c4+blk6*8*c2;tmp6=blk0*8*c4-blk2*8*c2;tmp7=blk4*8*c4-blk6*8*c6;tmp8=blk1*8*c7-blk3*8*c5;tmp9=blk5*8*c3-blk7*8*c1;tmp10=blk1*8*c5-blk3*8*c1;tmp11=blk5*8*c7+blk7*8*c3;tmp12=blk1*8*c3-blk3*8*c7;tmp13=-blk5*8*c1-blk7*8*c5;tmp14=blk1*8*c1+blk3*8*c3;tmp15=blk5*8*c5+blk7*8*c7;/second steptmp0=0.5*(tmp0+tmp1);tmp1=0.5*(tmp2+tmp3);tmp2=0.5*(tmp4+tmp5);tmp3=0.5*(tmp6+tmp7);tmp4=0.5*(tmp8+tmp9);tmp5=0.5*(tmp10+tmp11);tmp6=0.5*(tmp12+tmp13);tmp7=0.5*(tmp14+tmp15);/third stepblk0*8=tmp0+tmp7;blk1*8=tmp1+tmp6;blk2*8=tmp2+tmp5;blk3*8=tmp3+tmp4;blk4*8=tmp3-tmp4;blk5*8=tmp2-tmp5;blk6*8=tmp1-tmp6;blk7*8=tmp0-tmp7;/*blk0*8=0.5*(y0c4+y2c2+y4c4+y6c6+y1c1+y3c3+y5c5+y7c7);blk1*8=0.5*(y0c4+y2c6-y4c4-y6c2+y1c3-y3c7-y5c1-y7c5);blk2*8=0.5*(y0c4-y2c6-y4c4+y6c2+y1c5-y3c1+y5c7+y7c3);blk3*8=0.5*(y0c4-y2c2+y4c4-y6c6+y1c7-y3c5+y5c3-y7c1);blk4*8=0.5*(y0c4-y2c2+y4c4-y6c6-y1c7+y3c5-y5c3+y7c1);blk5*8=0.5*(y0c4-y2c6-y4c4+y6c2-y1c5+y3