代码仓库：代码、嵌入提取使用的图像、jpeg_tool库、实验报告_Gitee。
实验环境：MATLAB 2022a。
LSB空域隐写
原理
我们知道，一个像素点是由R(red)B(blue)G(green)三原色组成的，通过调配这三种颜色，我们可以得到所有的颜色。
比如白色(255,255,255)，二进制就是bin(11111111,11111111,11111111)，黑色(0,0,0)，二进制就是(00000000,00000000,00000000)。
既然如此，我们将每个二进制的最后一位给替换成别的，比如将(255,255,255)替换成(254,254,254)，bin(11111110,11111110,11111110)，肉眼根本分辨不出。
因此我们将需要的信息转成二进制，再将每一位替换掉三元组的最后一位，便完成了LSB隐写。
值对现象原理
但是，LSB隐写有明显的统计特征——“值对现象”。
因为密文具有随机性，
可认为密文的二进制比特流中的0和1是均匀的。
将二进制的最后一位替换成密文的比特后，
原图像的灰度值的奇数和偶数将会趋于相同。
比如，颜色值为100的本来有10个，为101的有100个，在嵌入随机密文后，它们的个数有可能趋于相同——也就是颜色值为100的和101的都变成55个。
实验内容
载体选取：BMP格式灰度图像。
嵌入信息：JPG图像。
嵌入大小：227680bits。
若BMP图像不方便寻找，可用MATLAB将jpg图像转换成BMP：
% 将jpg转换成bmp图像
A=imread('avatar.jpg');
imwrite(A,'cover.bmp','bmp');
嵌入信息程序：
clear % 清空变量
close all % 关闭打开的图像窗口
Picture = imread('cover.bmp');
Double_Picture = Picture;
Double_Picture = double(Double_Picture);
%读取秘密信息文件为二进制数字流，为嵌入图像做准备
wen.txt_id = fopen ('avatar2.jpg','r');
[msg, len] = fread(wen.txt_id, 'ubit1');
%根据LSB算法，将秘密信息的二进制数字流隐藏入连续的像素中
[m, n]=size(Double_Picture);
p=1;
for f2 = 1:n
for f1 = 1:m
Double_Picture(f1, f2) = Double_Picture(f1, f2)-mod(Double_Picture(f1,f2), 2)+msg(p,1);
if p == len
break;
end
p = p+1;
end
if p == len
break;
end
end
%将得到的隐密图像保存为stego.bmp，并利用Matlab将载体图与隐密图画在同一对话框中进行比较
Double_Picture=uint8(Double_Picture);
imwrite(Double_Picture, 'stego1.bmp');
subplot(121);imshow(Picture);title('未嵌入信息的图像');
subplot(122);imshow(Double_Picture);title('嵌入信息的图像');
fclose(wen.txt_id);
运行结果：
提取程序：
clear %清空变量
% 利用Matlab自带函数读取隐密图像stego. bmp，得到隐密图像的信息并将图像转化为二进制
Picture = imread( 'stego1.bmp' );
Picture = double(Picture);
[m, n] = size(Picture);
% 打开存放秘密信息的文件，若没有则新建一个文件。顺序提取图像相应像素LSB的秘密信息，存储在打开的文件中并保存
frr = fopen('message.jpg', 'w');
len = 227680;
%test = [];
p = 1;
for f2 = 1:n
for f1 = 1:m
%lowbit = bitand(Picture(f1, f2), 1);
%test = [test lowbit]; % 注意！如果做数组cat操作，会严重增加提取程序运行时长
fwrite(frr, bitand(Picture(f1, f2), 1), 'ubit1');
if p == len
break;
end
p=p+1;
end
if p == len
break;
end
end
%fwrite(frr, test);
fclose(frr);
提取结果存储为message.jpg。
注意，不要尝试将内容全部写入一个数组，然后再统一用
fwrite写入，因为matlab的数组cat操作很慢，数据量太大了耗时很长。我多加了一行cat操作，结果给老师检查的时候翻车了qwq。
值对现象展示程序：
I=imread('cover.bmp');
O=imread('stego1.bmp');
I=I(1:200000);
O=O(1:200000);
%subplot(121);
x=100;
histogram(I, 0:1:x);
set(gca, 'xtick', 0:2:x); % 横坐标每隔2显示刻度
grid on;
hold on;
%subplot(122);
histogram(O, 0:1:x);
set(gca, 'xtick', 0:2:x); % 横坐标每隔2显示刻度
grid on;
值对现象运行结果：
上图中，蓝色是嵌入前，橙色是嵌入后。现象极其明显。
DCT域隐写
实验要求，完成3种嵌入方案：JSTEG、F4和F5。
载体选取：JPEG格式灰度图像。
嵌入信息：随机生成的0-1比特流。
嵌入长度：JSteg：28684bits；F4：45000bits；F5:32700bits。(长度与图像本身能够用于嵌入的AC系数个数、算法对AC系数的利用率有关)
嵌入思路：
本实验使用jpeg_toolbox库对JPEG图像进行读写操作，因此，若jpeg_toolbox文件夹尚未添加到当前路径，将会出现如下报错：
解决办法：右键文件夹-添加到路径-选定的文件夹。
JSteg
JSteg信息隐藏算法是LSB替换思想在DCT域的实现。
嵌入过程的关键步骤：将原始图像的AC系数中最低的位平面“替换”为要隐藏的秘密信息。这里的“替换”遵循如下规则：
（1）忽略-1、0、1；
（2）若AC系数为2i，秘密比特为0，该系数不变；
（3）若AC系数为2i，秘密比特为1，该系数变为2i+1；
（4）若AC系数为2i+1，秘密比特为0，该系数变为2i；
（5）若AC系数为2i+1，秘密比特为1，该系数不变。
（6）AC系数为负数时，其二进制的实际含义是正数。(如表2-1所示)
在本课程中，DCT系数为负时，认为它是相应正数的补码。
因此，负数的奇偶性和其他正常的JSteg说法不一样。
如，AC系数为(-16)?? = (01111)?。当秘密比特为0时，载密系数是(-17)?? = (01110)?。
嵌入程序stego_JSteg.m：
% 更改嵌入算法时，需要将下文JSteg_simulation替换成其他函数
% 并修改算法名称name、嵌入信息长度messageLen等变量
COVER='cover.jpg';
STEGO='stego.jpg';
name='JSteg';
messageLen=28684;
message=randi([0 1],1,messageLen);%0000); %生成随机数，作为隐藏信息
save('message','message','-ascii'); %保存秘密信息
tic;
[nzAC]=JSteg_simulation(COVER,STEGO,message);
T=toc;
fprintf('-----------------------------------\n');
fprintf('%s simulation finished\n', name);
fprintf('cover image: %s\n',COVER);
fprintf('stego image: %s\n',STEGO);
fprintf('number of nzACs in cover: %i\n',nzAC);
fprintf('bits of message: %d\n',length(message));
fprintf('elapsed time: %.4f seconds\n',T);
fprintf('-----------------------------------\n');
嵌入程序调用的函数JSteg_simulation.m：
function [ncAC]=JSteg_simulation(COVER,STEGO,message)
try
jobj=jpeg_read(COVER); %读取cover图片
PrimeDCT=jobj.coef_arrays{1};%读取DCT系数
DCT=PrimeDCT;
catch
error('ERROR (problem with the cover image)');
end
AC_Location=DCT; % 复制DCT
AC_Location(1:8:end,1:8:end)=false; % 将DC系数置0
AC_Location(abs(AC_Location)<=1)=0; % 将绝对值小于等于1的位置置为0
AC_Location=find(AC_Location); %找出DCT中不为DC系数、绝对值大于1的位置
ncAC=numel(AC_Location); % 得到绝对值大于1的AC系数个数
messageLen=length(message);
%信息过长
if(messageLen>ncAC)
error('ERROR (too long message)');
end
i_DCT=1;
for i_MSG=1:messageLen
if(i_DCT>ncAC) % 没有更多可供注入密文的AC系数
fprintf('Max messageLength is %d.\n', i_MSG-1);
error('ERROR (too long message)');
end
DCTInfo=DCT(AC_Location(i_DCT));
if(DCTInfo>0)
DCT(AC_Location(i_DCT))=DCTInfo+message(i_MSG)-mod(DCTInfo,2);
else
DCT(AC_Location(i_DCT))=DCTInfo+message(i_MSG)-mod(DCTInfo+1,2);
end
i_DCT=i_DCT+1;
end
%%% save the resulting stego image
try
jobj.coef_arrays{1} = DCT;
jobj.optimize_coding = 1;
jpeg_write(jobj,STEGO);
catch
error('ERROR (problem with saving the stego image)');
end
%显示图像
subplot(2,2,1);imshow(COVER);
title('未嵌入信息的图像');
subplot(2,2,2);histogram(PrimeDCT);axis([-10,10,0,2*1e4]);
title('嵌入前的图像DCT系数直方图');
subplot(2,2,3);imshow(STEGO);
title('嵌入信息的图像');
subplot(2,2,4);histogram(DCT);axis([-10,10,0,2*1e4]);
title('嵌入后的图像DCT系数直方图');
运行结果：
提取程序extract_JSteg.m：
% 更改提取算法时，需要将下文JSteg_extract替换成其他函数
% 并修改算法名称name、嵌入信息长度messageLen等变量
STEGO='stego.jpg';
name='JSteg';
messageLen=28684;
tic;
messageHiden=JSteg_extract(STEGO,messageLen);
T=toc;
save('messageHiden','messageHiden','-ascii'); %保存提取出来的秘密信息
fprintf('-----------------------------------\n');
fprintf('%s extract finished\n', name);
fprintf('elapsed time: %.4f seconds\n',T);
fprintf('-----------------------------------\n');
提取程序调用的提取函数JSteg_extract.m：
function message=JSteg_extract(STEGO,messageLen)
try
jobj=jpeg_read(STEGO); %读取stego图片
DCT=jobj.coef_arrays{1};%读取DCT系数
catch
error('ERROR (problem with the cover image)');
end
AC_Location=DCT; % 复制DCT
AC_Location(1:8:end,1:8:end)=false; % 将DC系数置0
AC_Location(abs(AC_Location)<=1)=0; % 将绝对值小于等于1的位置置为0
AC_Location=find(AC_Location); %找出DCT中不为DC系数、绝对值大于1的位置
i_DCT=1;
for i_MSG=1:messageLen
DCTInfo=DCT(AC_Location(i_DCT));
if(DCTInfo>0)
message(1,i_MSG)=mod(DCTInfo,2);
else
message(1,i_MSG)=mod(DCTInfo+1,2);
end
i_DCT=i_DCT+1;
end
后面的两个算法不贴出完整代码，代码仓库里面都有。也懒得解释实验现象了，实验报告都给得很清楚。
简单介绍一下剩下两个实验的原理吧。
F4
F4信息隐藏算法是JSteg的改进。
在F4算法中，当嵌入系数的最低位和密文不符时，将系数的绝对值一律减1，而非对最低位直接取反，从而较好地避免了值对现象。
它的替换，和JSteg的区别是：
①它对1和-1也进行操作。当嵌入系数的最低位和密文不符时，1和-1都将变成0。此时，该位的信息嵌入被视为无效，需要继续取下一位AC系数进行嵌入。
②嵌入时，统一直接减1。如原值是6，对应二进制是110，在JSteg中，嵌入后是111；在F4中，嵌入后是101。
体现出来的变换规则如下表2-2：
F5
F5的特点是矩阵编码，能够减小数据的修改量，并置乱DCT系数（本实验无置乱要求）。
矩阵编码：将k比特秘密消息嵌入到(2^k-1)个AC系数中，只修改1个位置。
其实，它的原理就和海明码纠错原理一样。
海明码（也叫汉明码）具有一位纠错能力。和海明码有些不同的是，矩阵编码不要求将AC系数计算得到的校验码添加到AC系数中去纠错，而是以计算得到的校验码作为秘密信息。
假设一共有H1,H2,H3,H4,H5,H6,...,H15这十五个数，它们的位置对应着二进制0001,0010,0011,…,1111。
将所有二进制最低位为1位置的数选出来异或，也就是b1=H1⊕H3⊕H5⊕H7⊕H9⊕H11⊕H13⊕H15，能得到1比特信息。
同理，将所有二进制第二位为1位置的数选出来，也就是b2=H2⊕H3⊕H6⊕H7⊕H10⊕H9⊕H9，也能得到1比特信息。
…
综上，一共能得到4比特信息。
倘若这4比特和4比特密文相符合，那就不用改密文。
假如不同，b1、b2、b3、b4的某些与密文不符合的可能组合一共有15种，这15种组合分别对应着需要修改的AC系数的位置。如下图所示。
这就是F5的修改原则——每(2^k-1)位AC系数，只需要最多改1位，就能得到k比特的密文信息。其中AC系数的分组大小可自行选择，修改的方式与F4保持一致，也是绝对值直接减1。
在本实验中，我选择以3个AC系数为一组，一次嵌入密文比特量k为2，分别定义为a1,a2,a3,b1,b2。嵌入过程的关键步骤是矩阵编码，遵循如下规则：
（1）忽略0；
（2）当系数被修改成0时，换位置重新嵌入。
（3）b1=a1⊕a2, b2=a2⊕a3，则不修改数据；
（4）b1≠a1⊕a2, b2=a2⊕a3，则修改a1；
（5）b1＝a1⊕a2, b2≠a2⊕a3，则修改a3；
（6）b1≠a1⊕a2, b2≠a2⊕a3，则修改a2；
（7）正奇数和负偶数代表秘密消息1，正偶数和负奇数代表秘密消息0；
（8）AC系数为负数时，其二进制的实际含义是正数。
提取秘密消息时，只需令b1=a1⊕a2, b2=a2⊕a3。
简单贴一下F5的信息隐藏和提取关键代码：
信息隐藏的关键代码：
AC_Location=DCT; % 复制DCT
AC_Location(1:8:end,1:8:end)=false; % 将DC系数置0
AC_Location=find(AC_Location); %找出DCT中不为DC系数、不为0的位置
ncAC=numel(AC_Location); % 得到AC系数总数
messageLen=length(message); % 得到密文的数量
%信息过长
if(messageLen/2>ncAC/3)
error('ERROR (too long message)');
end
i_DCT=1;
i_MSG=1;
while i_MSG+1<=messageLen
if(i_DCT+2>ncAC) % 没有更多可供注入密文的AC系数
fprintf('Max messageLength is %d.\n', i_MSG-2);
error('ERROR (too long message)');
end
DCTInfo=DCT(AC_Location(i_DCT:i_DCT+2));
DCTInfo(DCTInfo<0)=DCTInfo(DCTInfo<0)+1; % 将所有负数+1，改变最低位
xor1=bitxor(mod(DCTInfo(1),2),mod(DCTInfo(2),2)); % a1按位异或a2
xor2=bitxor(mod(DCTInfo(2),2),mod(DCTInfo(3),2)); % a2按位异或a3
tobe_change=i_DCT;
if(message(i_MSG)==xor1) % 判断异或值，确定修改位
if(message(i_MSG+1)==xor2)
tobe_change=-1;
else
tobe_change=tobe_change+2;
end
else
if(message(i_MSG+1)~=xor2)
tobe_change=tobe_change+1;
end
end
if(tobe_change~=-1) % 需要修改
if(DCT(AC_Location(tobe_change))>0) % 正数减1
DCT(AC_Location(tobe_change))=DCT(AC_Location(tobe_change))-1;
else % 负数加1
DCT(AC_Location(tobe_change))=DCT(AC_Location(tobe_change))+1;
end
if(DCT(AC_Location(tobe_change))==0) % 如果嵌入后是0，这两位密文重新嵌入
i_MSG=i_MSG-2;
AC_Location(tobe_change)=[]; % 删除这一位置
ncAC=ncAC-1;
i_DCT=i_DCT-3; % 复用未被修改的AC系数
end
end
i_DCT=i_DCT+3;
i_MSG=i_MSG+2; % for循环内循环索引改变不生效
end
若AC系数分组较大，如选择分组大小为15，则不建议继续采用
if-else写“确定修改位”部分，而应该反向打表。
首先将异或值不同的位按序直接转换成二进制，比如b1,b2,b3不同，应直接计算得到1110，然后查1110对应的是ACT系数7，则直接修改第七个ACT系数。
提取的关键代码：
AC_Location=DCT; % 复制DCT
AC_Location(1:8:end,1:8:end)=false; % 将DC系数置0
AC_Location=find(AC_Location); %找出DCT中不为DC系数、不为0的位置
i_DCT=1;
for i_MSG=1:2:messageLen
DCTInfo=DCT(AC_Location(i_DCT:i_DCT+2));
DCTInfo(DCTInfo<0)=DCTInfo(DCTInfo<0)+1; % 将所有负数+1后，改变最低位
message(1,i_MSG)=bitxor(mod(DCTInfo(1),2),mod(DCTInfo(2),2)); % a1按位异或a2
message(1,i_MSG+1)=bitxor(mod(DCTInfo(2),2),mod(DCTInfo(3),2)); % a2按位异或a3
i_DCT=i_DCT+3;
end
可以看到，提取算法非常简单。其实嵌入算法也是很简单的，就是在F4的基础上，加了一个异或判断修改位。
