MATLAB高级编程与工程应用人脸识别实验报告+源代码

1、MATLAB高级编程与工程应用 实验四 图像处理第一章 基础知识1、MATLAB 提供了图像处理工具箱，在命令窗口输入help images 可查看该工具箱内的所有函数。请阅读并大致了解这些函数的基本功能。大致了解。2、利用MATLAB 提供的Image file I/O 函数分别完成以下处理：(a)以测试图像的中心为圆心，图像的长和宽中较小值的一半为半径画一个红颜色的圆；分析：直接利用半径条件，满足条件的点将红色元素置为255，绿色和蓝色元素置为0，于是得到如下图像：源代码：load(hall_color.mat);%首先获得三原数组R = hall_color(:,:,1);G = hal

2、l_color(:,:,2);B = hall_color(:,:,3);%将圆上的点改为红色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(hall_color1);(

3、b)将测试图像涂成国际象棋状的“黑白格”的样子，其中“黑”即黑色，“白”则意味着保留原图。用一种看图软件浏览上述两个图，看是否达到了目标。分析：首先设置标记flag在进行循环，对不同方格实行颜色更改就行。效果：源代码：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 = 1) for m = 15*(i - 1) +

4、 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、图像的预处理是将每个像素灰度值减去128 ，这个步骤是否可以在变换域进行？请在测试图像中截取一块验证你的结论。分析：可以在变换域进行，这个操作对应

5、于在变换域将直流分量减去128*8*8/8 = 1024，于是可以得到如下图像：原图为：直接将灰度值减去128得到图像为：通过改变变换域数据的方式得到图像为：对比发现，两种变换方式得到的预处理后图像都是一样的。另外需要说明一下，为了让预处理后的图像显示出来不全为黑色，我将原图像的灰度像素均乘以了一个2，这样可以更加方便的对比变换后的图像。2、请编程实现二维DCT ，并和MATLAB 自带的库函数dct2 比较是否一致。分析：由PDF中的补充知识我们可以知道，任意二维矩阵做DCT变换的时候都可以由它和与之大小相同的余弦序列矩阵加权表示，也就是说C = DPDT，那么我们先找出相应的D矩阵，于是我

6、写出程序循环得到了相应维度的D矩阵，然后可以直接使用矩阵乘法得到相应的二维DCT变换，这样与MATLAB自带的库函数dct2虽然在代码实现上有一定的差异（dct2中调用了一维DCT函数实现二维DCT），但是结果实际上是一致的。得到DCT系数矩阵的代码如下所示：%生成8*8 DCT系数矩阵clc;%参数定义及初始化N = 8;DCT = zeros(N,N);%计算DCT系数矩阵DCT(1,:) = sqrt(1/N);for i = 2:N; for j = 1:N DCT(i,j) = sqrt(2/N)*cos(i - 1)*(2*j - 1)*pi/(2*N); endend3、如果将D

7、CT 系数矩阵中右侧四列的系数全部置零，逆变换后的图像会发生什么变化？选取一块图验证你的结论。如果左侧的四列置零呢？分析：检测图像如下图所示：如果DCT系数矩阵中右侧四列的系数全部置零，逆变后的图像高频成分将会被抹去，得到图像如下图所示：比较两个图，可以发现在黑灰、白灰以及黑白交界处有些许变化，右下角的变化较为明显，之后我将两个图像矩阵作差，这样可以更加直观的看到除去高频成分后的差异：很明显的我们可以看到在上面所谓的交界处差值绝对值交大，而越往相邻颜色差异不大的块走，那个差值的绝对值越小。进一步调整。如果DCT系数矩阵中左侧四列的系数全部置零，逆变后的图像直流分量以及低频分量都被抹去，得到的图

8、像如下图所示：可以看到，完全就是一片黑的。查看逆变后的矩阵我发现，两次逆变后的矩阵相加之和就是原来的矩阵了，也就是说他们是互补的。4、若对DCT 系数分别做转置、旋转90 度和旋转180 度操作(rot90) ，逆变换后恢复的图像有何变化？选取一块图验证你的结论。分析：我选择如下图像作验证：首先，将DCT系数作转置，逆变换后恢复得到的图像如下图所示：明显我们可以看到逆变换后恢复的图像是原图像逆时针旋转了90度后的结果，原因是对DCT系数转置，原来各横行频率信息变为各竖列频率信息，对应于原图像的横竖调换。进一步，将DCT系数做旋转90度操作，逆变换后恢复得到的图像如下图所示：可以看到图像发生了较

9、大变化，分析其原因主要是DCT系数旋转90度后，频域内原来的竖列直流信息的行列最高频分量变成了图像的直流分量，这个分量较小，而原图像直流分量变成了各列最高频分量的直流分量，这个分量比较大，还原回去的新图像各竖列的波动较大，而整体灰度较暗。下面， 将DCT系数做旋转90度操作，逆变换后恢复得到的图像如下图所示：观察发现图像的变化更大了，黑白交错得更厉害，分析其原因在于将DCT系数旋转180度后，原图像直流分量编程各横行竖列的最高频率分量，这个量是比较大的，而原图像横行竖列的最高频率分量变为了图像的直流分量，这个量是相当小的，所以说最终逆变后的图像沿横竖两个方向的灰度波动都交大，并且图像整体灰度较

10、暗。源代码如下：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、如果认为差分编码是一个系统，请绘出这个系统的频率响应，说明它是一个_高通_（低通、高通、带通、带阻）滤波器。DC 系数先进行差分编码再进行熵编码，说明DC 系数的_低_频率分量更多。分析：首先写出系统的差分方程如下：y(n) = x(n-1) x(n)对方程进行Z

11、变换后得到如下传递函数：然后画出频率相应得到：由此可知系统呈现高通特性，说明DC系数的低频分量更多，通过差分编码后可以尽可能多的减短码长。源代码如下：image_2_5.mclc,close all;a = 1;b = -1 1;tf(b,a);figure;freqz(b,a);6、DC预测误差的取值和Category 值有何关系？如何利用预测误差计算出其Category ？分析：Category值就是DC预测误差的绝对值的二进制位长，可以得到一个转换公式：7、你知道哪些实现Zig-Zag扫描的方法？请利用MATLAB的强大功能设计一种最佳方法。分析：Zigzag扫描实现将每一个8*8矩阵的

12、DCT系数按照特定的顺序存储成一列数据，但是在实现扫描的Z字型移动时，会遇到较多麻烦，在编程实现时，可以建立一个映射关系，采取顺序扫描，离散映射的方法快速获得列矢量各元素。下表表示矩阵各系数对应编号：12671516282935814172730434913182631424410121925324145541120243340465355212334394752566122353848515760623637495058596364 按照二维顺序扫描此8*8图像，然后按照对应编号将每个分量映射到矢量的对应位置。 写出函数zigzag()源代码如下所示：zigzag.m:function Z

13、= 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 = zeros(64,1);for i = 1:8 for j = 1:8 Z(order(i,j) = C(i,j); endendend8、对测试图像分块、DCT和量化，将量

14、化后的系数写成矩阵的形式，其中每一列为一个块的DCT系数Zig-Zag扫描后形成的列矢量，第一行为各个块的DC系数。分析：首先将测试图像分成互补重叠的8*8小块，然后依次对每一个8*8小块进行DCT变换，得到8*8的DCT系数矩阵，系数矩阵除以量化矩阵得到根据权重量化的DCT系数，接着调用7问中的zigzag()函数将量化后的DCT系数按顺序存储成为64*1的列向量，最后将所有列向量按行并起来即可得到用来表征整个图像的DCT系数矩阵，这个整体图像DCT系数矩阵的第一行即为DC系数，其余为AC系数。原代码如下：image_2_8.mclc,close all,clear all;load Jpe

15、gCoeff.mat;load hall.mat;%实现分块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为分块数目for i = 1:a*b/64 Blocki = double(Blocki) - 128; %预处理 C = dct2(Blocki); %DCT变换 C_qual = round(C./QTAB); %量化 J(:,i) = zigzag(C_qua

16、l); %zigzag扫描end9、请实现本章介绍的JPEG编码（不包括写JFIF文件），输出为DC系数的码流、AC系数的码流、图像高度和图像宽度，将这四个变量写入Jpegcodes.mat文件。分析：按照pdf中所给的编码方式分别写出DC编码和AC编码的函数，然后将使用8问方法处理过后的矢量调用这两个函数分别得到了DC系数的码流和AC系数的码流。源代码如下：clc,close all,clear all;load JpegCoeff.mat;load hall.mat;%实现分块a b = size(hall_gray);M = 8*ones(1,a/8);N = 8*ones(1,b/8)

17、;Block = mat2cell(hall_gray,M,N);%DCT变换、量化、zigzag扫描 J = zeros(64,a*b/64); %a*b/64为分块的数目for i = 1:a*b/64 Blocki = double(Blocki) - 128; %预处理 C = dct2(Blocki); %DCT变换 C_qual = round(C./QTAB); %量化 J(:,i) = zigzag(C_qual); %zigzag扫描end%熵编码DC = DCencode(J(1,:);%调用函数DCencode实现直流编码AC = ACencode(J); %调用函数Ce

18、ncode实现交流编码%保存save Jpegcodes.mat DC AC a b;DCencode.mfunction DC = DCencode(dc)%读入DCTABload JpegCoeff.mat;%进行差分运算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 = zeros(0);for i = 1:b %计算各category，用于定位DCTAB中对应行数 if dc_dif(i) = 0; category = floor(log

19、2(abs(dc_dif(i) + 1; else category = 0; end %提取category的huffman码 len = DCTAB(category + 1,1); dc_huff = DCTAB(category + 1,2:len + 1); %生成相应二进制码 dc_bin = de2bi(abs(dc_dif(i),left-msb); if(dc_dif(i) 0) dc_bin = dc_bin; end %生成完整编码 if(category = 0) DC = DC,dc_huff; else DC = DC,dc_huff,dc_bin; endende

20、ndACencode.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系数 place = find(J(2:end,i); %定位此序列中向量非0系数所在下标 last = place(end) + 1; %校正下标值，应加1 Run = 0; for j = 2:last if(J(j,i) = 0) R

21、un = 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);%得到对应huffman编码 ac_bin = de2bi(abs(J(j,i),left-msb); %转化为二进制数 if(J(j,i) 0) %对负数取反 ac_bin = ac_bin; end AC = AC,ac_huff,ac_bin; Run = 0; en

22、d end end AC = AC,EOB; %插入块结束符endend10、计算压缩比（输入文件长度/输出码流长度），注意转换为相同进制。分析：输出码流长度 = DC码流长度 + AC码流长度 + 高的二进制位长 + 宽的二进制位长，由此可以得到输出码流长度 = 2173 + 23072 + 8 + 8 = 25261，而输入文件长度 = 120*168*8 = 161280，最终我们可以得到压缩比 = 161280/25261 6.385。11、请实现本章介绍的JPEG解码，输入是你生成的Jpegcodes.mat文件。分别用客观(PSNR)和主观方式评价编解码效果如何。分析：按照编码的逆

23、顺序处理了自己生成的Jpegcodes.mat文件，生成如下图像：按照课件中的方法编程计算出了PSNR系数的大小：由此不管是从主管观察还是客观的用PSNR系数来判断，都可以看到解码出来的图像与原图像相比的效果是不太好的。源代码如下：image_2_11.mclc,close all,clear all;%载入相关数据load JPEGCodes.mat;load JpegCoeff.mat;%解码J_dc = DCdecode(DC,a,b);J_ac = ACdecode(AC,a,b,ACTAB);J = cell2mat(J_dc;J_ac);%第一行是直流分量，之后为交流分量%反zig

24、zag扫描，反量化，反DCT变换Block = cell(a/8,b/8);for i = 1:a*b/64 Q = Izigzag(J(:,i);%调用自定义反zigzag函数，恢复量化系数矩阵 C = Q.*QTAB; %反量化 Blocki = idct2(C);%反DCT变换 Blocki = Blocki + 128; %反预处理end%拼接恢复图像hall_gray2 = cell2mat(Block);hall_gray2 = uint8(hall_gray2);save hall_gray2.mat hall_gray2;%画出图像进行对比load hall.matfigure

25、;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 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

26、 64;C = zeros(8,8);for i = 1:64 C(i) = J(order(i);endendDCdecode.mfunction J = DCdecode(DC,a,b)%解码恢复预测误差序列J_dif = ones(1,a*b/64);i = 1;j = 1;while i = length(DC); %计算category if(DC(i) = 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

27、= 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) category = 6; i = i + 4; elseif(DC(i + 4) = 0) category = 7; i = i + 5; elseif(DC(i + 5) = 0) category = 8; i = i + 6;

28、 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 %二进制转化为十进制 if category = 0 J_dif(j) = 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,l

29、eft-msb); j = j + 1; else J_dif(j) = bi2de(DC_bin,left-msb); j = j + 1; end i = i + category; endend%恢复DC系数J = zeros(1,a*b/64);J(1) = J_dif(1);for i = 2:a*b/64 J(i) = 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

30、/64);blocknum = 1;%块位置编号i = 1; %列位置编号j = 1; %熵编码位置编号while(j length(AC) if(all(AC(j:j + 3) = EOB) blocknum = blocknum + 1; %进入下一个块 i = 1; %重置列位置编号 j = j + 4; %熵编码位置往后推进 elseif(all(AC(j:j + 10) = ZRL) i = i + 16;%相当于填16个零 j = j + 11; else for x = 1:160 len = ACTAB(x,3); %将符合条件的列找出来，从而确定原数据信息 if(j + le

31、n = length(AC)&all(AC(j:j + len - 1) = ACTAB(x,4:len + 3) Run = ACTAB(x,1); Size = ACTAB(x,2); i = i + Run;%填充相应个数的零 j = j + len; AC_bin = AC(j:j + Size - 1); if(AC_bin(1) = 0)%对负数的情况单独处理 AC_bin = AC_bin; J(i,blocknum) = - bi2de(AC_bin,left-msb); i = i + 1; else J(j,blocknum) = bi2de(AC_bin,left-msb

32、); i = i + 1; end j = j + Size; break;%找到对应列后可以直接跳出循环，提高转换效率 end end endendend12、将量化步长减小为原来的一半，重做编解码。同标准量化步长的情况比较压缩比和图像质量。分析：只需要将代码中的QTAB整体除以2即可。代码在文件image_2_12.m中。最终恢复出来的图像为：压缩比 = 120*168*8/(2576 + 34164 + 16) 4.388直观上看，两种步长的画质都差不多，但是PSNR有细微差别：13、看电视时偶尔能看到美丽的雪花图像（见snow.mat），请对其编解码。和测试图像的压缩比和图像质量进行比

33、较，并解释比较结果。分析：只需要将12问中的hall_gray变量赋值为snow，并且将步长改回QTAB即可，得到文件image_2_13.m，运行文件可得到如下图像:压缩比 = 120*168*8/(43546 + 1426 + 16) 3.585观察数据发现AC系数的编码数较多，这是因为图像块内黑白交替过快，高频分量较多，导致无法有效进行压缩，这是导致压缩比偏小的主要原因。第三章 信息隐藏1、 实现本章介绍的空域隐藏方法和提取方法。验证其抗JPEG编码能力。分析：假设要隐藏的数据为：“Matlab is an interesting and useful tool”。隐藏后图像效果如下所示

34、：可以看到图像是几乎没有发生变化的，所以说空域隐藏信息的隐蔽性比较好。但是经过JPEG编码之后再解码，就无法从解码后的图像中还原得到信息，最后恢复得到的信息如下所示：words = ?i?4?t?%j4?%Ual?j? K?c?$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); %将隐藏信息依次转换成8位二进制码w_bin = r

35、eshape(w_bin.,1); %生成二进制码流word_len = length(words); %记录隐藏信息字节数hidden = de2bi(uint8(word_len),8);w_bin;%信息前八位用于记录信息字节数%将隐藏信息添加到像素值最低位hall_bin = de2bi(hall_gray);%由于函数默认将二进制方向调换，所以最低位被放在了首列hall_bin(1:length(hidden),1) = hidden;hall_hidden = bi2de(hall_bin);hall_hidden = reshape(hall_hidden,size(hall_g

36、ray);figure;subplot(1,2,1),imshow(hall_gray);subplot(1,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为分块的数目for i = 1:a*b/64 Blocki = double(Blocki) - 128; %预处理 C = dct2(Blocki);

37、 %DCT变换 C_qual = round(C./QTAB); %量化 J(:,i) = zigzag(C_qual); %zigzag扫描end%熵编码DC = DCencode(J(1,:);%调用函数DCencode实现直流编码AC = ACencode(J); %调用函数ACencode实现交流编码 %解码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

38、 i = 1:a*b/64 Q = Izigzag(J(:,i);%调用自定义反zigzag函数，恢复量化系数矩阵 C = Q.*(QTAB/2); %反量化 Blocki = idct2(C);%反DCT变换 Blocki = Blocki + 128; %反预处理end%拼接恢复出图像hall_hidden2 = cell2mat(Block);hall_hidden2 = uint8(hall_hidden2);%提取隐藏信息码hidden = de2bi(hall_hidden2);word_len = bi2de(hidden(1:8,1);%提取隐藏信息字节数len = doubl

39、e(word_len)*8 + 8;word_bin = hidden(9:len,1);%恢复隐藏信息word_bin = reshape(word_bin,8,);words = char(bi2de(word_bin);disp(words = words);2、 依次实现本章介绍的三种变换域信息隐藏方法和提取方法，分析嵌密方法的隐蔽性以及嵌密后JPEG图像的质量变化和压缩比变化。分析：只需考虑实现2、3方法，由于此时的信息是隐藏在变换域中，JPEG的熵编码是对经过量化并舍入为整数的DCT系数的完全无损可逆编码，所以说这样做是完全可以抵抗JPEG编码的。首先使用信息位逐一替换掉若干量化后

40、的DCT系数的最低位，得到的加密后的图像如下右图所示：对比发现加密后图像变得原图像模糊了，最后尝试恢复原来的信息，得到：图像的PNSR值为：压缩比 = 120*168*8/(16 + 23072 + 2173) 6.385进一步，采用逐一将信息为追加在每个块Zigzag顺序的最后一个非零DCT系数之后，得到加密后的图像如右下所示：加密前后差异主要体现在图像的上端，房顶与云的部分变得比较模糊，尝试恢复原来的信息，得到：图像的PNSR值为：压缩比 = 120*168*8/(25717 + 2187 + 16) 5.777源代码如下：Tesk3_2_2.mclc;%产生待处理DCT系数矩阵load

41、JpegCoeff.mat;load hall.mat;%实现分块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为分块的数目for i = 1:a*b/64 Blocki = double(Blocki) - 128; %预处理 C = dct2(Blocki); %DCT变换 C_qual = round(C./QTAB); %量化 J(:,i) = zigzag(C

42、_qual); %zigzag扫描endwords = Matlab is an interesting and useful tool;w_bin = de2bi(uint8(words),8); %将隐藏信息转换成二进制形式w_bin = reshape(w_bin,1); %生成二进制码流w_len = length(words);hidden = de2bi(uint8(w_len),8);w_bin;%前八位记录隐藏信息字节数J_r = reshape(abs(J),1);J_bin = de2bi(J_r,8);J_bin(1:length(hidden),1) = hidden;

43、 %逐一替换每个量化后的DCT系数的二进制最低位J_hidden = bi2de(J_bin);J_hidden = reshape(J_hidden,size(J);J_hidden(J 0) = - J_hidden(J 0); %将负数部分恢复%恢复出加密后的图像Block = cell(a/8,b/8);for i = 1:a*b/64 Q = Izigzag(J(:,i);%调用自定义反zigzag函数，恢复量化系数矩阵 C = Q.*QTAB; %反量化 Blocki = idct2(C);%反DCT变换 Blocki = Blocki + 128; %反预处理end%拼接恢复出图

44、像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);%解码恢复出加密的信息J_r = reshape(abs(J_hidden),1);J_bin = de2bi(J_r,8);w_len = bi2de(J_bin(1:8,1)*8;%信息字节数乘以8，得到信息占用的总二进制位数w_bin = J_bin(

45、9:w_len + 8,1);w_bin = reshape(w_bin,8,);words = char(bi2de(w_bin);disp(words = ,words);tesk3_2_3.mclc;%产生待处理DCT系数矩阵load JpegCoeff.mat;load hall.mat;%实现分块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为分块的数目for

46、 i = 1:a*b/64 Blocki = double(Blocki) - 128; %预处理 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); %将隐藏信息转换成二进制形式w_bin = reshape(w_bin,1); %生成二进制码流w_len = length(words);hidden =

47、de2bi(uint8(w_len),8);w_bin;%前八位记录隐藏信息字节数J_r = reshape(abs(J),1);J_bin = de2bi(J_r,8);%逐一替换每个量化后的最后一个非零DCT系数之后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); %将负数部分恢复%恢复出加密后的图像Block = cell(a/8,b/8);for i = 1:a*b/64 Q =