MATLAB 高級編程與工程應(yīng)用 人臉識別 實(shí)驗(yàn)報(bào)告+源代碼

1、MATLAB高級編程與工程應(yīng)用實(shí)驗(yàn)四 圖像處理第一章 基礎(chǔ)知識1、MATLAB 提供了圖像處理工具箱，在命令窗口輸入help images 可查看該工具箱內(nèi)的所有函數(shù)。請閱讀并大致了解這些函數(shù)的基本功能。大致了解。2、利用MATLAB 提供的Image file I/O 函數(shù)分別完成以下處理：(a)以測試圖像的中心為圓心，圖像的長和寬中較小值的一半為半徑畫一個(gè)紅顏色的圓；分析：直接利用半徑條件，滿足條件的點(diǎn)將紅色元素置為255，綠色和藍(lán)色元素置為0，于是得到如下圖像：源代碼：load('hall_color.mat');%首先獲得三原數(shù)組R = hall_color(:,:,1

2、);G = hall_color(:,:,2);B = hall_color(:,:,3);%將圓上的點(diǎn)改為紅色for i = 1:120 for j = 1:168 a = abs(i - 60.5); b = abs(j - 84.5); d = sqrt(a 2 + b 2); if(abs(d - 60) < 0.5) R(i,j) = 255; G(i,j) = 0; B(i,j) = 0; end endend%生成新的矩陣hall_color1(:,:,1) = R;hall_color1(:,:,2) = G;hall_color1(:,:,3) = B; imshow(

3、hall_color1);(b)將測試圖像涂成國際象棋狀的“黑白格”的樣子，其中“黑”即黑色，“白”則意味著保留原圖。用一種看圖軟件瀏覽上述兩個(gè)圖，看是否達(dá)到了目標(biāo)。分析：首先設(shè)置標(biāo)記flag在進(jìn)行循環(huán)，對不同方格實(shí)行顏色更改就行。效果：源代碼：clear all; load('hall_color.mat'); R = hall_color(:,:,1);G = hall_color(:,:,2);B = hall_color(:,:,3); flag = 1;for i = 1:8 flag = mod(flag + 1),2); for j = 1:8 if(flag =

4、 1) for m = 15*(i - 1) + 1:15*i for n = 21*(j - 1) + 1:21*j R(m,n) = 0; G(m,n) = 0; B(m,n) = 0; end end end flag = mod(flag + 1),2); endend hall_color1(:,:,1) = R;hall_color1(:,:,2) = G;hall_color1(:,:,3) = B;imshow(hall_color1);用看圖軟件打開成功：第二章 圖像壓縮編碼1、圖像的預(yù)處理是將每個(gè)像素灰度值減去128 ，這個(gè)步驟是否可以在變換域進(jìn)行？請?jiān)跍y試圖像中截取一塊驗(yàn)

5、證你的結(jié)論。分析：可以在變換域進(jìn)行，這個(gè)操作對應(yīng)于在變換域?qū)⒅绷鞣至繙p去128*8*8/8 = 1024，于是可以得到如下圖像：原圖為：直接將灰度值減去128得到圖像為：通過改變變換域數(shù)據(jù)的方式得到圖像為：對比發(fā)現(xiàn)，兩種變換方式得到的預(yù)處理后圖像都是一樣的。另外需要說明一下，為了讓預(yù)處理后的圖像顯示出來不全為黑色，我將原圖像的灰度像素均乘以了一個(gè)2，這樣可以更加方便的對比變換后的圖像。2、請編程實(shí)現(xiàn)二維DCT ，并和MATLAB 自帶的庫函數(shù)dct2 比較是否一致。分析：由PDF中的補(bǔ)充知識我們可以知道，任意二維矩陣做DCT變換的時(shí)候都可以由它和與之大小相同的余弦序列矩陣加權(quán)表示，也就是說C

6、= DPD¬T，那么我們先找出相應(yīng)的D矩陣，于是我寫出程序循環(huán)得到了相應(yīng)維度的D矩陣，然后可以直接使用矩陣乘法得到相應(yīng)的二維DCT變換，這樣與MATLAB自帶的庫函數(shù)dct2雖然在代碼實(shí)現(xiàn)上有一定的差異（dct2中調(diào)用了一維DCT函數(shù)實(shí)現(xiàn)二維DCT），但是結(jié)果實(shí)際上是一致的。得到DCT系數(shù)矩陣的代碼如下所示：%生成8*8 DCT系數(shù)矩陣clc;%參數(shù)定義及初始化N = 8;DCT = zeros(N,N);%計(jì)算DCT系數(shù)矩陣DCT(1,:) = sqrt(1/N);for i = 2:N; for j = 1:N DCT(i,j) = sqrt(2/N)*cos(i - 1)*(2

7、*j - 1)*pi/(2*N); endend3、如果將DCT 系數(shù)矩陣中右側(cè)四列的系數(shù)全部置零，逆變換后的圖像會發(fā)生什么變化？選取一塊圖驗(yàn)證你的結(jié)論。如果左側(cè)的四列置零呢？分析：檢測圖像如下圖所示：如果DCT系數(shù)矩陣中右側(cè)四列的系數(shù)全部置零，逆變后的圖像高頻成分將會被抹去，得到圖像如下圖所示：比較兩個(gè)圖，可以發(fā)現(xiàn)在黑灰、白灰以及黑白交界處有些許變化，右下角的變化較為明顯，之后我將兩個(gè)圖像矩陣作差，這樣可以更加直觀的看到除去高頻成分后的差異：很明顯的我們可以看到在上面所謂的交界處差值絕對值交大，而越往相鄰顏色差異不大的塊走，那個(gè)差值的絕對值越小。進(jìn)一步調(diào)整。如果DCT系數(shù)矩陣中左側(cè)四列的系數(shù)

8、全部置零，逆變后的圖像直流分量以及低頻分量都被抹去，得到的圖像如下圖所示：可以看到，完全就是一片黑的。查看逆變后的矩陣我發(fā)現(xiàn)，兩次逆變后的矩陣相加之和就是原來的矩陣了，也就是說他們是互補(bǔ)的。4、若對DCT 系數(shù)分別做轉(zhuǎn)置、旋轉(zhuǎn)90 度和旋轉(zhuǎn)180 度操作(rot90) ，逆變換后恢復(fù)的圖像有何變化？選取一塊圖驗(yàn)證你的結(jié)論。分析：我選擇如下圖像作驗(yàn)證：首先，將DCT系數(shù)作轉(zhuǎn)置，逆變換后恢復(fù)得到的圖像如下圖所示：明顯我們可以看到逆變換后恢復(fù)的圖像是原圖像逆時(shí)針旋轉(zhuǎn)了90度后的結(jié)果，原因是對DCT系數(shù)轉(zhuǎn)置，原來各橫行頻率信息變?yōu)楦髫Q列頻率信息，對應(yīng)于原圖像的橫豎調(diào)換。進(jìn)一步，將DCT系數(shù)做旋轉(zhuǎn)90度

9、操作，逆變換后恢復(fù)得到的圖像如下圖所示：可以看到圖像發(fā)生了較大變化，分析其原因主要是DCT系數(shù)旋轉(zhuǎn)90度后，頻域內(nèi)原來的豎列直流信息的行列最高頻分量變成了圖像的直流分量，這個(gè)分量較小，而原圖像直流分量變成了各列最高頻分量的直流分量，這個(gè)分量比較大，還原回去的新圖像各豎列的波動較大，而整體灰度較暗。下面， 將DCT系數(shù)做旋轉(zhuǎn)90度操作，逆變換后恢復(fù)得到的圖像如下圖所示：觀察發(fā)現(xiàn)圖像的變化更大了，黑白交錯得更厲害，分析其原因在于將DCT系數(shù)旋轉(zhuǎn)180度后，原圖像直流分量編程各橫行豎列的最高頻率分量，這個(gè)量是比較大的，而原圖像橫行豎列的最高頻率分量變?yōu)榱藞D像的直流分量，這個(gè)量是相當(dāng)小的，所以說最終逆

10、變后的圖像沿橫豎兩個(gè)方向的灰度波動都交大，并且圖像整體灰度較暗。源代碼如下：image_2_4.mclc;C = dct2(P);C1 = C'C2 = rot90(C);C3 = rot90(C2);imtool(idct2(C),0 255);imtool(idct2(C1),0 255);imtool(idct2(C2),0 255);imtool(idct2(C3),0 255);5、如果認(rèn)為差分編碼是一個(gè)系統(tǒng)，請繪出這個(gè)系統(tǒng)的頻率響應(yīng)，說明它是一個(gè)_高通_（低通、高通、帶通、帶阻）濾波器。DC 系數(shù)先進(jìn)行差分編碼再進(jìn)行熵編碼，說明DC 系數(shù)的_低_頻率分量更多。分析：首先寫出

11、系統(tǒng)的差分方程如下：y(n) = x(n-1) x(n)對方程進(jìn)行Z變換后得到如下傳遞函數(shù)：Hz= Y(z)X(z)=1-zz然后畫出頻率相應(yīng)得到：由此可知系統(tǒng)呈現(xiàn)高通特性，說明DC系數(shù)的低頻分量更多，通過差分編碼后可以盡可能多的減短碼長。源代碼如下：image_2_5.mclc,close all;a = 1;b = -1 1;tf(b,a);figure;freqz(b,a);6、DC預(yù)測誤差的取值和Category 值有何關(guān)系？如何利用預(yù)測誤差計(jì)算出其Category ？分析：Category值就是DC預(yù)測誤差的絕對值的二進(jìn)制位長，可以得到一個(gè)轉(zhuǎn)換公式：Category=floorlog

12、2DC+1,DC00,DC=07、你知道哪些實(shí)現(xiàn)Zig-Zag掃描的方法？請利用MATLAB的強(qiáng)大功能設(shè)計(jì)一種最佳方法。分析：Zigzag掃描實(shí)現(xiàn)將每一個(gè)8*8矩陣的DCT系數(shù)按照特定的順序存儲成一列數(shù)據(jù)，但是在實(shí)現(xiàn)掃描的Z字型移動時(shí)，會遇到較多麻煩，在編程實(shí)現(xiàn)時(shí)，可以建立一個(gè)映射關(guān)系，采取順序掃描，離散映射的方法快速獲得列矢量各元素。下表表示矩陣各系數(shù)對應(yīng)編號：126715162829358141727304349131826314244101219253241455411202433404653552123343947525661223538485157606236374950585963

13、64 按照二維順序掃描此8*8圖像，然后按照對應(yīng)編號將每個(gè)分量映射到矢量的對應(yīng)位置。 寫出函數(shù)zigzag()源代碼如下所示：zigzag.m:function Z = zigzag(C)order = 1 2 6 7 15 16 28 29; 3 5 8 14 17 27 30 43; 4 9 13 18 26 31 42 44; 10 12 19 25 32 41 45 54; 11 20 24 33 40 46 53 55; 21 23 34 39 47 52 56 61; 22 35 38 48 51 57 60 62; 36 37 49 50 58 59 63 63;Z = zero

14、s(64,1);for i = 1:8 for j = 1:8 Z(order(i,j) = C(i,j); endendend8、對測試圖像分塊、DCT和量化，將量化后的系數(shù)寫成矩陣的形式，其中每一列為一個(gè)塊的DCT系數(shù)Zig-Zag掃描后形成的列矢量，第一行為各個(gè)塊的DC系數(shù)。分析：首先將測試圖像分成互補(bǔ)重疊的8*8小塊，然后依次對每一個(gè)8*8小塊進(jìn)行DCT變換，得到8*8的DCT系數(shù)矩陣，系數(shù)矩陣除以量化矩陣得到根據(jù)權(quán)重量化的DCT系數(shù)，接著調(diào)用7問中的zigzag()函數(shù)將量化后的DCT系數(shù)按順序存儲成為64*1的列向量，最后將所有列向量按行并起來即可得到用來表征整個(gè)圖像的DCT系數(shù)矩

15、陣，這個(gè)整體圖像DCT系數(shù)矩陣的第一行即為DC系數(shù)，其余為AC系數(shù)。原代碼如下：image_2_8.mclc,close all,clear all;load JpegCoeff.mat;load hall.mat;%實(shí)現(xiàn)分塊a b = size(hall_gray);M = 8*ones(1,a/8);N = 8*ones(1,b/8);Block = mat2cell(hall_gray,M,N);%DCT變換、量化、zigzag掃描J = zeros(64,a*b/64); %a*b/64為分塊數(shù)目for i = 1:a*b/64 Blocki = double(Blocki) - 12

16、8; %預(yù)處理 C = dct2(Blocki); %DCT變換 C_qual = round(C./QTAB); %量化 J(:,i) = zigzag(C_qual); %zigzag掃描end9、請實(shí)現(xiàn)本章介紹的JPEG編碼（不包括寫JFIF文件），輸出為DC系數(shù)的碼流、AC系數(shù)的碼流、圖像高度和圖像寬度，將這四個(gè)變量寫入Jpegcodes.mat文件。分析：按照pdf中所給的編碼方式分別寫出DC編碼和AC編碼的函數(shù)，然后將使用8問方法處理過后的矢量調(diào)用這兩個(gè)函數(shù)分別得到了DC系數(shù)的碼流和AC系數(shù)的碼流。源代碼如下：clc,close all,clear all;load JpegCoe

17、ff.mat;load hall.mat;%實(shí)現(xiàn)分塊a b = size(hall_gray);M = 8*ones(1,a/8);N = 8*ones(1,b/8);Block = mat2cell(hall_gray,M,N);%DCT變換、量化、zigzag掃描 J = zeros(64,a*b/64); %a*b/64為分塊的數(shù)目for i = 1:a*b/64 Blocki = double(Blocki) - 128; %預(yù)處理 C = dct2(Blocki); %DCT變換 C_qual = round(C./QTAB); %量化 J(:,i) = zigzag(C_qual)

18、; %zigzag掃描end%熵編碼DC = DCencode(J(1,:);%調(diào)用函數(shù)DCencode實(shí)現(xiàn)直流編碼AC = ACencode(J); %調(diào)用函數(shù)Cencode實(shí)現(xiàn)交流編碼%保存save Jpegcodes.mat DC AC a b;DCencode.mfunction DC = DCencode(dc)%讀入DCTABload JpegCoeff.mat;%進(jìn)行差分運(yùn)算a b = size(dc);dc_dif = zeros(1,b);dc_dif(1) = dc(1);dc_dif(2:end) = dc(1:end - 1) - dc(2:end);DC = zero

19、s(0);for i = 1:b %計(jì)算各category，用于定位DCTAB中對應(yīng)行數(shù) if dc_dif(i) = 0; category = floor(log2(abs(dc_dif(i) + 1; else category = 0; end %提取category的huffman碼 len = DCTAB(category + 1,1); dc_huff = DCTAB(category + 1,2:len + 1); %生成相應(yīng)二進(jìn)制碼 dc_bin = de2bi(abs(dc_dif(i),'left-msb'); if(dc_dif(i) < 0) d

20、c_bin = dc_bin; end %生成完整編碼 if(category = 0) DC = DC,dc_huff; else DC = DC,dc_huff,dc_bin; endendendACencode.mfunction AC = ACencode(J)%讀入ACTABload JpegCoeff.mat;%初始化ZRL = 1 1 1 1 1 1 1 1 0 0 1;EOB = 1 0 1 0;a b = size(J);AC = zeros(0);for i = 1:b if(any(J(2:end,i) %如果此子向量中存在非0系數(shù) place = find(J(2:en

21、d,i); %定位此序列中向量非0系數(shù)所在下標(biāo) last = place(end) + 1; %校正下標(biāo)值，應(yīng)加1 Run = 0; for j = 2:last if(J(j,i) = 0) Run = Run + 1; if(Run = 16) AC = AC,ZRL; Run = 0; end else Size = floor(log2(abs(J(j,i) + 1; len = ACTAB(Run*10 + Size,3); ac_huff = ACTAB(Run*10 + Size,4:len + 3);%得到對應(yīng)huffman編碼 ac_bin = de2bi(abs(J(j,i

22、),'left-msb'); %轉(zhuǎn)化為二進(jìn)制數(shù) if(J(j,i) < 0) %對負(fù)數(shù)取反 ac_bin = ac_bin; end AC = AC,ac_huff,ac_bin; Run = 0; end end end AC = AC,EOB; %插入塊結(jié)束符endend10、計(jì)算壓縮比（輸入文件長度/輸出碼流長度），注意轉(zhuǎn)換為相同進(jìn)制。分析：輸出碼流長度 = DC碼流長度 + AC碼流長度 + 高的二進(jìn)制位長 + 寬的二進(jìn)制位長，由此可以得到輸出碼流長度 = 2173 + 23072 + 8 + 8 = 25261，而輸入文件長度 = 120*168*8 = 161

23、280，最終我們可以得到壓縮比 = 161280/25261 6.385。11、請實(shí)現(xiàn)本章介紹的JPEG解碼，輸入是你生成的Jpegcodes.mat文件。分別用客觀(PSNR)和主觀方式評價(jià)編解碼效果如何。分析：按照編碼的逆順序處理了自己生成的Jpegcodes.mat文件，生成如下圖像：按照課件中的方法編程計(jì)算出了PSNR系數(shù)的大?。河纱瞬还苁菑闹鞴苡^察還是客觀的用PSNR系數(shù)來判斷，都可以看到解碼出來的圖像與原圖像相比的效果是不太好的。源代碼如下：image_2_11.mclc,close all,clear all;%載入相關(guān)數(shù)據(jù)load JPEGCodes.mat;load Jpeg

24、Coeff.mat;%解碼J_dc = DCdecode(DC,a,b);J_ac = ACdecode(AC,a,b,ACTAB);J = cell2mat(J_dc;J_ac);%第一行是直流分量，之后為交流分量%反zigzag掃描，反量化，反DCT變換Block = cell(a/8,b/8);for i = 1:a*b/64 Q = Izigzag(J(:,i);%調(diào)用自定義反zigzag函數(shù)，恢復(fù)量化系數(shù)矩陣 C = Q.*QTAB; %反量化 Blocki = idct2(C);%反DCT變換 Blocki = Blocki + 128; %反預(yù)處理end%拼接恢復(fù)圖像hall_g

25、ray2 = cell2mat(Block);hall_gray2 = uint8(hall_gray2);save hall_gray2.mat hall_gray2;%畫出圖像進(jìn)行對比load hall.matfigure;subplot(1,2,1),imshow(hall_gray,0 255);subplot(1,2,2),imshow(hall_gray2);izigzag.mfunction C = Izigzag(J)order = 1 2 6 7 15 16 28 29; 3 5 8 14 17 27 30 43; 4 9 13 18 26 31 42 44; 10 12 1

26、9 25 32 41 45 54; 11 20 24 33 40 46 53 55; 21 23 34 39 47 52 56 61; 22 35 38 48 51 57 60 62; 36 37 49 50 58 59 63 64;C = zeros(8,8);for i = 1:64 C(i) = J(order(i);endendDCdecode.mfunction J = DCdecode(DC,a,b)%解碼恢復(fù)預(yù)測誤差序列J_dif = ones(1,a*b/64);i = 1;j = 1;while i <= length(DC); %計(jì)算category if(DC(i)

27、 = 0) if(DC(i + 1) = 0) category = 0; i = i + 2; elseif(DC(i + 2) = 0) category = 1; i = i + 3; else category = 2; i = i + 3; end elseif(DC(i + 1) = 0) if(DC(i + 2) = 0) category = 3; i = i + 3; else category = 4; i = i + 3; end elseif(DC(i + 2) = 0); category = 5; i = i + 3; elseif(DC(i + 3) = 0) c

28、ategory = 6; i = i + 4; elseif(DC(i + 4) = 0) category = 7; i = i + 5; elseif(DC(i + 5) = 0) category = 8; i = i + 6; elseif(DC(i + 6) = 0) category = 9; i = i + 7; elseif(DC(i + 7) = 0) category = 10; i = i + 8; elseif(DC(i + 8) = 0) category = 11; i = i + 9; end %二進(jìn)制轉(zhuǎn)化為十進(jìn)制 if category = 0 J_dif(j)

29、 = 0; j = j + 1; else DC_bin = DC(i:i + category - 1); if(DC_bin(1) = 0) DC_bin = DC_bin; J_dif(j) = - bi2de(DC_bin,'left-msb'); j = j + 1; else J_dif(j) = bi2de(DC_bin,'left-msb'); j = j + 1; end i = i + category; endend%恢復(fù)DC系數(shù)J = zeros(1,a*b/64);J(1) = J_dif(1);for i = 2:a*b/64 J(i

30、) = J(i - 1) - J_dif(i);endendACdecode.mfunction J = ACdecode(AC,a,b,ACTAB)%初始化ZRL = 1 1 1 1 1 1 1 1 0 0 1;EOB = 1 0 1 0;J = zeros(63,a*b/64);blocknum = 1;%塊位置編號i = 1; %列位置編號j = 1; %熵編碼位置編號while(j < length(AC) if(all(AC(j:j + 3) = EOB) blocknum = blocknum + 1; %進(jìn)入下一個(gè)塊 i = 1; %重置列位置編號 j = j + 4; %

31、熵編碼位置往后推進(jìn) elseif(all(AC(j:j + 10) = ZRL) i = i + 16;%相當(dāng)于填16個(gè)零 j = j + 11; else for x = 1:160 len = ACTAB(x,3); %將符合條件的列找出來，從而確定原數(shù)據(jù)信息 if(j + len <= length(AC)&&all(AC(j:j + len - 1) = ACTAB(x,4:len + 3) Run = ACTAB(x,1); Size = ACTAB(x,2); i = i + Run;%填充相應(yīng)個(gè)數(shù)的零 j = j + len; AC_bin = AC(j:j

32、 + Size - 1); if(AC_bin(1) = 0)%對負(fù)數(shù)的情況單獨(dú)處理 AC_bin = AC_bin; J(i,blocknum) = - bi2de(AC_bin,'left-msb'); i = i + 1; else J(j,blocknum) = bi2de(AC_bin,'left-msb'); i = i + 1; end j = j + Size; break;%找到對應(yīng)列后可以直接跳出循環(huán)，提高轉(zhuǎn)換效率 end end endendend12、將量化步長減小為原來的一半，重做編解碼。同標(biāo)準(zhǔn)量化步長的情況比較壓縮比和圖像質(zhì)量。分析：

33、只需要將代碼中的QTAB整體除以2即可。代碼在文件image_2_12.m中。最終恢復(fù)出來的圖像為：壓縮比 = 120*168*8/(2576 + 34164 + 16) 4.388直觀上看，兩種步長的畫質(zhì)都差不多，但是PSNR有細(xì)微差別：13、看電視時(shí)偶爾能看到美麗的雪花圖像（見snow.mat），請對其編解碼。和測試圖像的壓縮比和圖像質(zhì)量進(jìn)行比較，并解釋比較結(jié)果。分析：只需要將12問中的hall_gray變量賦值為snow，并且將步長改回QTAB即可，得到文件image_2_13.m，運(yùn)行文件可得到如下圖像:壓縮比 = 120*168*8/(43546 + 1426 + 16) 3.585

34、觀察數(shù)據(jù)發(fā)現(xiàn)AC系數(shù)的編碼數(shù)較多，這是因?yàn)閳D像塊內(nèi)黑白交替過快，高頻分量較多，導(dǎo)致無法有效進(jìn)行壓縮，這是導(dǎo)致壓縮比偏小的主要原因。第三章 信息隱藏1、 實(shí)現(xiàn)本章介紹的空域隱藏方法和提取方法。驗(yàn)證其抗JPEG編碼能力。分析：假設(shè)要隱藏的數(shù)據(jù)為：“Matlab is an interesting and useful tool”。隱藏后圖像效果如下所示：可以看到圖像是幾乎沒有發(fā)生變化的，所以說空域隱藏信息的隱蔽性比較好。但是經(jīng)過JPEG編碼之后再解碼，就無法從解碼后的圖像中還原得到信息，最后恢復(fù)得到的信息如下所示：words = ?i?4?t?%j4?%¨Ual?j? K?ìc

35、?$yü?f +wR&/?ü?f ?n?6?r?<?ò?üLk4?j ?êU.?顯然，這樣的隱藏方法根本無法抵抗JPEG編碼。原代碼如下：clc;load hall.mat;words = 'Matlab is an interesting useful tool'%生成隱藏信息w_bin = de2bi(uint8(words),8); %將隱藏信息依次轉(zhuǎn)換成8位二進(jìn)制碼w_bin = reshape(w_bin.',1); %生成二進(jìn)制碼流word_len = length(words); %記錄隱藏信

36、息字節(jié)數(shù)hidden = de2bi(uint8(word_len),8)'w_bin;%信息前八位用于記錄信息字節(jié)數(shù)%將隱藏信息添加到像素值最低位hall_bin = de2bi(hall_gray);%由于函數(shù)默認(rèn)將二進(jìn)制方向調(diào)換，所以最低位被放在了首列hall_bin(1:length(hidden),1) = hidden;hall_hidden = bi2de(hall_bin);hall_hidden = reshape(hall_hidden,size(hall_gray);figure;subplot(1,2,1),imshow(hall_gray);subplot(1

37、,2,2),imshow(hall_hidden); a b = size(hall_hidden);M = 8*ones(1,a/8);N = 8*ones(1,b/8);Block = mat2cell(hall_hidden,M,N); %DCT變換、量化、zigzag掃描 J = zeros(64,a*b/64); %a*b/64為分塊的數(shù)目for i = 1:a*b/64 Blocki = double(Blocki) - 128; %預(yù)處理 C = dct2(Blocki); %DCT變換 C_qual = round(C./QTAB); %量化 J(:,i) = zigzag(C

38、_qual); %zigzag掃描end%熵編碼DC = DCencode(J(1,:);%調(diào)用函數(shù)DCencode實(shí)現(xiàn)直流編碼AC = ACencode(J); %調(diào)用函數(shù)ACencode實(shí)現(xiàn)交流編碼 %解碼J_dc = DCdecode(DC,a,b);J_ac = ACdecode(AC,a,b,ACTAB);J = cell2mat(J_dc;J_ac);%第一行是直流分量，之后為交流分量%反zigzag掃描，反量化，反DCT變換Block = cell(a/8,b/8);for i = 1:a*b/64 Q = Izigzag(J(:,i);%調(diào)用自定義反zigzag函數(shù)，恢復(fù)量化系

39、數(shù)矩陣 C = Q.*(QTAB/2); %反量化 Blocki = idct2(C);%反DCT變換 Blocki = Blocki + 128; %反預(yù)處理end%拼接恢復(fù)出圖像hall_hidden2 = cell2mat(Block);hall_hidden2 = uint8(hall_hidden2);%提取隱藏信息碼hidden = de2bi(hall_hidden2);word_len = bi2de(hidden(1:8,1)');%提取隱藏信息字節(jié)數(shù)len = double(word_len)*8 + 8;word_bin = hidden(9:len,1);%恢復(fù)

40、隱藏信息word_bin = reshape(word_bin,8,)'words = char(bi2de(word_bin)');disp('words = ' words);2、 依次實(shí)現(xiàn)本章介紹的三種變換域信息隱藏方法和提取方法，分析嵌密方法的隱蔽性以及嵌密后JPEG圖像的質(zhì)量變化和壓縮比變化。分析：只需考慮實(shí)現(xiàn)2、3方法，由于此時(shí)的信息是隱藏在變換域中，JPEG的熵編碼是對經(jīng)過量化并舍入為整數(shù)的DCT系數(shù)的完全無損可逆編碼，所以說這樣做是完全可以抵抗JPEG編碼的。首先使用信息位逐一替換掉若干量化后的DCT系數(shù)的最低位，得到的加密后的圖像如下右圖所示：

41、對比發(fā)現(xiàn)加密后圖像變得原圖像模糊了，最后嘗試恢復(fù)原來的信息，得到：圖像的PNSR值為：壓縮比 = 120*168*8/(16 + 23072 + 2173) 6.385進(jìn)一步，采用逐一將信息為追加在每個(gè)塊Zigzag順序的最后一個(gè)非零DCT系數(shù)之后，得到加密后的圖像如右下所示：加密前后差異主要體現(xiàn)在圖像的上端，房頂與云的部分變得比較模糊，嘗試恢復(fù)原來的信息，得到：圖像的PNSR值為：壓縮比 = 120*168*8/(25717 + 2187 + 16) 5.777源代碼如下：Tesk3_2_2.mclc;%產(chǎn)生待處理DCT系數(shù)矩陣load JpegCoeff.mat;load hall.mat

42、;%實(shí)現(xiàn)分塊a b = size(hall_gray);M = 8*ones(1,a/8);N = 8*ones(1,b/8);Block = mat2cell(hall_gray,M,N);%DCT變換、量化、zigzag掃描J = zeros(64,a*b/64); %a*b/64為分塊的數(shù)目for i = 1:a*b/64 Blocki = double(Blocki) - 128; %預(yù)處理 C = dct2(Blocki); %DCT變換 C_qual = round(C./QTAB); %量化 J(:,i) = zigzag(C_qual); %zigzag掃描endwords =

43、 'Matlab is an interesting and useful tool'w_bin = de2bi(uint8(words),8); %將隱藏信息轉(zhuǎn)換成二進(jìn)制形式w_bin = reshape(w_bin',1); %生成二進(jìn)制碼流w_len = length(words);hidden = de2bi(uint8(w_len),8)'w_bin;%前八位記錄隱藏信息字節(jié)數(shù)J_r = reshape(abs(J),1);J_bin = de2bi(J_r,8);J_bin(1:length(hidden),1) = hidden; %逐一替換每個(gè)量

44、化后的DCT系數(shù)的二進(jìn)制最低位J_hidden = bi2de(J_bin);J_hidden = reshape(J_hidden,size(J);J_hidden(J < 0) = - J_hidden(J < 0); %將負(fù)數(shù)部分恢復(fù)%恢復(fù)出加密后的圖像Block = cell(a/8,b/8);for i = 1:a*b/64 Q = Izigzag(J(:,i);%調(diào)用自定義反zigzag函數(shù)，恢復(fù)量化系數(shù)矩陣 C = Q.*QTAB; %反量化 Blocki = idct2(C);%反DCT變換 Blocki = Blocki + 128; %反預(yù)處理end%拼接恢復(fù)出

45、圖像hall_gray_hidden = cell2mat(Block);hall_gray_hidden = uint8(hall_gray_hidden);figure;subplot(1,2,1),imshow(hall_gray,0 255);subplot(1,2,2),imshow(hall_gray_hidden,0 255);%解碼恢復(fù)出加密的信息J_r = reshape(abs(J_hidden),1);J_bin = de2bi(J_r,8);w_len = bi2de(J_bin(1:8,1)')*8;%信息字節(jié)數(shù)乘以8，得到信息占用的總二進(jìn)制位數(shù)w_bin =

46、 J_bin(9:w_len + 8,1);w_bin = reshape(w_bin,8,)'words = char(bi2de(w_bin)');disp('words = ',words);tesk3_2_3.mclc;%產(chǎn)生待處理DCT系數(shù)矩陣load JpegCoeff.mat;load hall.mat;%實(shí)現(xiàn)分塊a b = size(hall_gray);M = 8*ones(1,a/8);N = 8*ones(1,b/8);Block = mat2cell(hall_gray,M,N);%DCT變換、量化、zigzag掃描J = zeros(6

47、4,a*b/64); %a*b/64為分塊的數(shù)目for i = 1:a*b/64 Blocki = double(Blocki) - 128; %預(yù)處理 C = dct2(Blocki); %DCT變換 C_qual = round(C./QTAB); %量化 J(:,i) = zigzag(C_qual); %zigzag掃描endwords = 'Matlab is an interesting and useful tool'w_bin = de2bi(uint8(words),8); %將隱藏信息轉(zhuǎn)換成二進(jìn)制形式w_bin = reshape(w_bin',1)

48、; %生成二進(jìn)制碼流w_len = length(words);hidden = de2bi(uint8(w_len),8)'w_bin;%前八位記錄隱藏信息字節(jié)數(shù)J_r = reshape(abs(J),1);J_bin = de2bi(J_r,8);%逐一替換每個(gè)量化后的最后一個(gè)非零DCT系數(shù)之后for n = 9:length(hidden) for x = 1:64 if(J(floor(n/8) + 1,x) = 0); break; end end J_bin(n,x) = hidden(n);endJ_hidden = bi2de(J_bin);J_hidden = reshape(J_hidden,size(J);J_hidden(J < 0) = - J_hidden(J < 0); %將負(fù)數(shù)部分恢復(fù)%恢復(fù)出加密后的圖像Block = cell(a/8,b/8);for i = 1