《物联网信息安全》实验指导书.doc

实验一 图像输入与输出基本操作（2学时）
实验二 基于DCT域的数字水印技术（2学时）
综合性实验 基于傅立叶域变换的数字水印技术（4学时）

实验前预备知识：
数字水印技术是利用数字产品普遍存在的冗余数据与随机性，将水印信息嵌入在数字产品本身中，从而起到保护数字产品版权或者完整性的一种技术。
现在学术界对数字水印算法的理解都是将一些不易察觉的具有随机特性的数据嵌入到图像频域或空域的系数上。从信号处理的角度看，嵌入水印可以看成是在强背景下迭加一个弱信号，由于人类视觉系统的分辨率受到一定的限制，只要迭加的信号幅度不超过HVS的对比门限，人眼就无法感觉到信号的存在，所以可以通过对原始图像进行一定调整，在不影响视觉效果的情况下嵌入一些水印信息。
数字水印系统的一般模型如图1所示：

图1 数字水印系统基本模型
水印嵌入器的输入量有三个：水印信号M，宿主信号S和密钥K。
水印信号M是指原始水印（图像或一个数字序列）通过一定的方法经过调制将嵌入到宿主信号中的数字信号。
宿主信号S是指被嵌入水印的信号（原始信号）。
密钥K则指用于提高水印系统安全性的密码信息，它独立于宿主信号。密钥有私有密钥和公共密钥之分，前者指攻击者在明确了水印嵌入方法但又不知道密钥的情况下，水印不会被破坏或盗取；后者是指攻击者对宿主信号（如内容标识、语言字幕等）不感兴趣的情况下，密钥也就不存在保密性，可以作为公共密钥。

实验一 图像输入与输出基本操作
一、实验题目：
图像输入与输出操作
二、实验目的
学习在MATLAB环境下对图像文件的I/O操作，为读取各种格式的图像文件和后续进行图像处理打下基础。
三、实验内容
利用MATLAB为用户提供的专门函数从图像格式的文件中读/写图像数据、显示图像，以及查询图像文件的信息。
四、预备知识
熟悉MATLAB开发环境。
五、实验原理
（1）图像文件的读取
利用imread函数可以完成图像文件的读取操作。常用语法格式为：
I=imread(‘filename’,‘fmt’)或I=imread(‘filename.fmt’)；
其作用是将文件名用字符串filename表示的、扩展名用字符串fmt（表示图像文件格式）表示的图像文件中的数据读到矩阵I中。当filename中不包含任何路径信息时，imread会从当前工作目录中寻找并读取文件。要想读取指定路径中的图像，最简单的方法就是在filename中输入完整的或相对的地址。 MATLAB支持多种图像文件格式的读、写和显示。因此参数fmt常用的可能值有：
‘bmp’ Windows位图格式
‘jpg’or‘jpeg’ 联合图像专家组格式
‘tif’or‘tiff’ 标志图像文件格式
‘gif’ 图形交换格式
‘pcx’ Windows画刷格式
‘png’ 可移动网络图形格式
‘xwd’ X Window Dump格式
例如，命令行
>>I=imread(‘lena.jpg’);
将JPEG图像lena读入图像矩阵I中。
（2） 图像文件的写入（保存）
利用imwrite完成图像的输出和保存操作，也完全支持也完全支持上述各种
图像文件的格式。其语法格式为：
imwrite(I,‘filename’,‘fmt’)或imwrite(I,‘filename.fmt’)；
其中的I、filename和fmt的意义同上所述。
注意事项：当利用imwrite函数保存图像时，MATLAB默认的保存方式是将其简化为uint8的数据类型。与读取文件类型类似，MATLAB在文件保存时还支持16位的PNG和TIFF图像。所以，当用户保存这类文件时，MATLAB就将其存储在uint16中。
（3）图像文件的显示
图像的现实过程是将数字图像从一组离散数据还原为一幅可见图像的过程。
MATLAB的的图像处理工具箱提供了多种图像显示技术。例如imshow可以直接从文件显示多种图像；image函数可以将矩阵作为图像 ；colorbar函数可以用来显示颜色条；montage函数可以动态显示图像序列。这里仅对常用的显示函数进行介绍。
①图像的显示
imshow函数是最常用的显示各种图像的函数，其调用格式如下：
imshow(I,N)；
imshow(I,N)用于显示灰度图像，其中I为灰度图像的数据矩阵，N为灰度级数目，默认值为256。
例如下面的语句用于显示一幅灰度图像：
>> I=imread(‘lena.jpg’)；
>> imshow(I);
如果不希望在显示图像之前装载图像，那么可以使用以下格式直接进行图像文件的显示：
imshow filename
其中，filename为要显示的图像文件的文件名。
实例1-1 显示一幅在当前目录下的.bmp格式的图像：
>>imshow lena.jpg
显示结果如图1.1所示。

图1.1.1 显示一幅图像文件中的图像
注意事项：该文件名必须带有合法的扩展名（指明文件格式），且该图像文件必须保存在当前目录下，或在MATLAB默认的目录下。
②添加色带
colorbar函数可以给一个坐标轴对象添加一条色带。如果该坐标轴对象包含一个图像对象，则添加的色带将指示出该图像中不同颜色的数据值。这对于了解被显示图像的灰度级特别有用。其调用格式为：
colorbar
实例1-2
>> I=imread(‘lena.jpg’)；
>> imshow(I);
>> colorbar;

图1.1.2 显示图像并加入颜色条
从上图可知，该图像是数据类型为uint8的灰度图像，其灰度级范围从0
－255。
③显示多幅图像
显示多幅图像最简单的方法就是在不同的图形窗口中显示它们。imshow总是在当前窗口中显示一幅图像，如果用户想连续显示两幅图像，那么第二幅图像就会替代第一幅图像。为了避免图像在当前窗口中的覆盖现象，在调用imshow函数显示下一幅图像之前可以使用figure命令来创建一个新的窗口。例如：
imshow(I1);
figure, imshow(I2);
figure, imshow(I3);
有时为了便于在多幅图像之间进行比较，需要将这些图像显示在一个图形窗口中。达到这一目的有两种方法：一种方法是联合使用imshow和subplot函数，但此方法在一个图形窗口只能有一个调色板；另一种方法是联合使用subimage和subplot函数，此方法可在一个图形窗口内使用多个调色板。
subplot函数将一个图形窗口划分为多个显示区域，其调用格式如下：
subplot(m,n,p)
subplot函数将图形窗口划分为m（行）×n（列）个显示区域，并选择第p个区域作为当前绘图区。
例1-3 用两排显示四幅图像，可以使用以下语句：
>> I1=imread(‘lena.bmp’)；
>> I2=imread(‘gs256.bmp’)；
>> I3=imread(‘lena.bmp’)；
>> I4=imread(‘gs256.bmp’)；
>> subplot（2,2,1）, imshow(I1);
>> subplot（2,2,2）, imshow(I2);
>> subplot（2,2,3）, imshow(I3);
>> subplot（2,2,4）, imshow(I4);

图1.1.3 在一个图形窗口中显示多幅图像
（4） 图像文件信息的查询
imfinfo函数用于查询图像文件的有关信息，详细地显示出图像文件的各种属性。其语法格式为：
info=imfinfo(‘filename’,‘fmt’)或info=imfinfo(‘filename.fmt’)
或imfinfo filename.fmt
imfinfo函数获取的图像文件信息依赖于文件类型的不同而不同，但至少应包
含以下内容：
文件名。如果该文件不在当前目录下，还包含该文件的完整路径。
文件格式。
文件格式的版本号。
文件最后一次修改的时间。
文件的大小。以字节为单位。
图像的宽度。
图像的高度。
每个像素所用的比特数。也叫像素深度。
图像类型。即该图像是真彩色图像、索引图像还是灰度图像。
例如，命令行
>>imfinfo bubbles25.jpg
会输出如下信息（注意，在这种情况下，有些域不包含信息）：
Filename: ‘bubbles25.jpg’
FileModDate: ‘04-Jan-2003 12:31:26’
FileSize: 13849
Format: ‘jpg’
FormatVersion: ‘’
Width: 714
Height: 682
BitDepth: 8
ColorType: ‘grayscale’
FormatSignature: ‘’
Comment: {}

六、实验步骤：
（1）利用imread函数完成对图像文件的读取操作。
（2）利用imwrite函数完成图像的写入（保存）操作。
（3）利用imshow函数显示图像。
（4）利用imfinfo函数查询图像文件的有关信息。

实验二 基于DCT域的数字水印技术
离散余弦变换，简称DCT，是一种实数域变换，其变换核为实数的余弦函数，计算速度较快，是一种近似最佳变换，很适合于做图象压缩和随机信号处理，它对常用的图象压缩有较强的鲁棒性。
基于DCT域的数字水印算法，一般而言有两种算法。一种是直接对整幅图像进行DCT计算，选取合适频段的系数，嵌入水印。另一种是，首先将整幅图像分成块，对每一块分别进行DCT计算，在每一块中选取合适频段的系数，将水印信息分散嵌入到每一块所选取的DCT系数中。这种方法称为分块DCT，这种方法非常适合于8×8的图像块DCT变换。目前，DCT域的水印方法大多数是采用分块DCT方法。
二维DCT概述：一般而言，数字图象,其二维DCT变换用矩阵表示，其定义式如下式所示：
（4-1）
二维IDCT定义如下：
（4-2）
其中 并且
（4-3）
我们最常用到的是8×8分块DCT，也就是首先将大小为M×N的图像分成的非重叠的大小为8×8的块，然后对每一块作DCT。根据式（4-1）和式（4-2），8×8DCT的定义如下：
（4-4）
8×8IDCT的定义如下：
（4-5）
其中， 并且
（4-6）
实验目的：学习使用基于DCT域数字水印嵌入算法和基于DCT域数字水印检测提取算法。
实验内容和原理：
（1）基于DCT域数字水印嵌入算法：
DCT域的水印嵌入过程，就是首先对宿主图像S进行分块DCT运算，选取出要嵌入的DCT系数C，并且对宿主图像用HVS模型进行分析，得到感知系数掩蔽模板，用其来控制水印嵌入强度，再将编码后的水印信号W，用一定方式嵌入到选好的系数中。用嵌入水印的系数替代原来位置的系数，再进行逆DCT（Inverse Discrete Cosine Transform,IDCT）变换，就得到了含有水印的图像X。DCT域水印系统的水印嵌入过程框图如图1所示：

图1 DCT域水印嵌入过程方框图
（2）基于DCT域数字水印检测提取算法：
DCT域水印的验证过程（即包括水印的检测和提取过程），就是将待测图像Y（注意此时的图像也许已经遭受了图象处理或恶意攻击）进行DCT运算；然后按着嵌入水印时选取的DCT系数，选取其含有水印的系数，进行水印提取，然后利用相关检测法判断水印是否存在。如果水印检测器输出结果显示水印存在，则根据需要可以用水印解码器对提取的水印进行解码，恢复水印。水印检测提取过程如图2

图2 DCT域数字水印检测提取框图
实验操作方法与步骤：
（1）水印嵌入步骤如下：
1. 将宿主图像分成8×8的块，每一块都进行DCT变换，得到与宿主图象相同尺寸的DCT域系数矩阵。
2. 产生两个不相关的伪随机序列，并设定一个密钥key。
3. 将每个8×8的DCT系数矩阵，从每一块的中频段取出个系数。
4. 然后嵌入水印，嵌入规则为当水印图像元素为‘0’时，按照式子（其中alpha为尺度因子，mark为随机序列，为原图像的分块DCT系数），将一个随机序列与原始图像块的幅度谱对应元素进行乘性叠加；为‘1’时，用另一个伪随机序列与幅度谱对应元素进行乘性叠加。
5. 用得到的新的DCT系数对原来位置的DCT系数进行置换。
6. 对新的DCT系数矩阵进行DCT逆变换，得到了嵌入水印信号的图像。
（2）水印检测提取步骤如下：
1. 对含有水印图像进行8×8的块DCT变换。
2. 对每一块的DCT系数进行“之”字型排列，从我们已知的嵌入位置取出可能含有水印的系数，将所有提取出的系数，按顺序组成新的系数序列。
3. 用密钥生成水印信号
4. 将含水印的系数序列与水印信号进行相关计算，得到归一化相关系数。
实验结果及分析：
以下是尺度因子alpha=30时的实验结果,以及受到各种攻击后所提取到的水印图像及其与原始水印图像的相关系数:

图3 原始图像和嵌入水印后的图像
此时，alpha=30

图4 没有受到攻击的图像及提取到的水印图像
此时相关系数=1

图5 加入白噪声后的图像、原始水印和提取的水印图像
此时相关系数=0.99969

图6 高斯低通滤波后的图像及提取的水印图像
此时相关系数=0.99058

图7 剪切后的图像及提取的水印图像
此时相关系数=093909

图8 旋转10后的图像及提取的水印图像
此时相关系数=0.48383
结论:当图像(宿主信息)受到加噪、高斯低通滤波、剪切等攻击后，其可视度、提取的水印图像与原始水印图像的相似度（都在90%以上）
都是非常令人满意的；当受到旋转攻击后，其提取的水印图像与原始水印图像的相似度（不足50%）却令人感到失望，不过仍可证明水印信息的存在。
程序参考：
%%%%%%%%%%%%% 嵌入算法 %%%%%%%%%%%%%%%
clear all;
clc;
%%%%%% 读取水印图像 %%%%%%%%
omark=double(imread('muxiao.bmp'))/255;
mo=size(omark,1);
no=size(omark,2);
%%%%% 以下生成水印信息 %%%%%%
mark=omark;
alpha=30;
rand1=randn(1,8);
rand2=randn(1,8);
%%%%%% 读取宿主图像 %%%%%%%%
cimage=imread('lena.bmp');
figure(1);
subplot(2,3,1),imshow(cimage,[]),title(‘原始的宿主图像’);
[mc,nc]=size(cimage);
cda0=blkproc(cimage,[8,8],'dct2');
%%%%%%% 嵌入水印 %%%%%%%%%
cda1=cda0;
for i=1:mo
for j=1:no
x=(i-1)*8;y=(j-1)*8;
if mark(i,j)==1
k=rand1;
else
k=rand2;
end
cda1(x+1,y+8)=cda0(x+1,y+8)+alpha*k(1);
cda1(x+2,y+7)=cda0(x+2,y+7)+alpha*k(2);
cda1(x+3,y+6)=cda0(x+3,y+6)+alpha*k(3);
cda1(x+4,y+5)=cda0(x+4,y+5)+alpha*k(4);
cda1(x+5,y+4)=cda0(x+5,y+4)+alpha*k(5);
cda1(x+6,y+3)=cda0(x+6,y+3)+alpha*k(6);
cda1(x+7,y+2)=cda0(x+7,y+2)+alpha*k(7);
cda1(x+8,y+1)=cda0(x+8,y+1)+alpha*k(8);
end
end
%%%%%% 嵌入水印后图象 %%%%%%
wimage=blkproc(cda1,[8,8],'idct2');
wimage_1=uint8(wimage);
imwrite(wimage_1,'withmark.bmp','bmp');
subplot(2,3,2),imshow(wimage,[]),title(‘嵌入水印后图象’);
%%%%%%%%%%%%%% 进行攻击测试 %%%%%%%%%%%
disp(‘1-->加入白噪声’);
disp(‘2-->高斯低通滤波’);
disp(‘3-->剪切图像’);
disp(‘4-->旋转攻击’);
disp(‘5-->直接检测’)
begin=input(‘请选择攻击（1-5）：’)
switch begin
%%%%%%% 加入白噪声 %%%%%%%%
case 1
Aimage1=wimage;
Wnoise=20*randn(size(Aimage1));
Aimage1=Aimage1+Wnoise;
subplot(2,3,4),imshow(Aimage1,[]),title(‘加入白噪声后的图象’);
att=Aimage1;
imwrite(att,'whitenoiseimage.bmp');
%%%%%%% 高斯低通滤波 %%%%%%%
case 2
Aimage2=wimage;
H=fspecial('gaussian',[4,4],0.5);
Aimage2=imfilter(Aimage2,H);
subplot(2,3,4),imshow(Aimage2,[]),title(‘高斯低通滤波后的图象’);
att=Aimage2;
imwrite(att,'gaussianimage.bmp')
%%%%%%%% 剪切攻击 %%%%%%%%
case 3
Aimage3=wimage;
Aimage3(1:128,1:128)=256;
subplot(2,3,4),imshow(Aimage3,[]),title(‘剪切后的图象’);
att=Aimage3;
imwrite(att,'cutpartimage.bmp');
%%%%%%%% 旋转攻击 %%%%%%%%
case 4
Aimage4=wimage;
Aimage4=imrotate(Aimage4,10,'bilinear','crop');
Aimage_4=mat2gray(Aimage4);
subplot(2,3,4),imshow(Aimage_4,[]),title(‘旋转10 度后的图象’);
att=Aimage_4;
imwrite(att,'rotatedimage.bmp');
%%%%%%% 没有受到攻击 %%%%%%%

case 5
subplot(2,3,4),imshow(wimage,[]),title(‘直接提取的图像’);
att=wimage;
imwrite(att,'directimage.bmp');
end
%%%%%%%%%%%%% 提取水印 %%%%%%%%%%%%%%
tmark_0=att;
tmark_0=blkproc(tmark_0,[8,8],'dct2');
pass=zeros(1,8);
for i=1:mo
for j=1:no
x=(i-1)*8;y=(j-1)*8;
pass(1)=tmark_0(x+1,y+8);
pass(2)=tmark_0(x+2,y+7);
pass(3)=tmark_0(x+3,y+6);
pass(4)=tmark_0(x+4,y+5);
pass(5)=tmark_0(x+5,y+4);
pass(6)=tmark_0(x+6,y+3);
pass(7)=tmark_0(x+7,y+2);
pass(8)=tmark_0(x+8,y+1);
if (corr2(pass,rand1)>corr2(pass,rand2))
tmark_1(i,j)=1;
else
tmark_1(i,j)=0;
end
end
end
%%%%%%%%%% 计算NC（归一化相关系数）%%%%%%%%
g_mark=double(tmark_1);
o_mark=double(omark);
[m,n]=size(g_mark);
nc_0=0;
nc_1=0;
nc_2=0;
for i=1:m
for j=1:n
nc_0=nc_0+g_mark(i,j)*o_mark(i,j);
nc_1=nc_1+o_mark(i,j)*o_mark(i,j);
nc_2=nc_2+g_mark(i,j)*g_mark(i,j);
end
end
NC=nc_0/sqrt(nc_1*nc_2);
%%%%%% 显示提取水印 %%%%%%%%
subplot(2,3,5),imshow('muxiao.bmp'),title(‘原始水印图像’);
subplot(2,3,6),imshow(tmark_1,[]);
name=‘提取得水印图像’;
title(strcat(num2str(name),'NC=',num2str(NC)));

综合设计性实验 基于傅立叶域变换的数字水印技术
实验目的：学习使用基于傅立叶域的数字水印嵌入算法和基于DCT域数字水印检测提取算法。
实验内容和原理：
傅立叶变换（Fourier Transform）是线性系统分析的一个有力工具，是研究信号的频谱方法，它架起了时域和频域之间的桥梁。使我们能够定量分析诸如数字化系统、采样点、电子放大器、卷积滤波器、噪声等的作用，把傅立叶变换的理论同其物理解释相结合，将大大有助于解决大多数图象处理问题。
Fourier分析理论十分完善，既可以处理连续信号也可以处理离散信号。计算机只能处理离散信号，且在数字图象处理中，输入图像和输出图像通常都是二维的，一般表示成二维数字矩阵，因此这里直接讨论二维离散傅立叶变换（DFT）和二维快速傅立叶变换（FFT）。
（1）二维离散傅立叶变换DFT（Discrete Fourier Transform）
二维离散傅立叶变换对定义式如下：
（1-1）
（1-2）
二维离散傅立叶变换的傅立叶谱、相位、功率谱分别如下：
傅立叶谱： （1-3）
相位： （1-4）
功率谱： （1-5）
式子（1-1）可以分离为：
（1-6）
式子（1-2）可以分离为：
（1-7）
由上两式可知，二维傅立叶变换可由连续两次一维傅立叶变换得到，从而将二维DFT分解为水平和垂直两部分运算，上式中方括号中的项表示在图像的行上计算的DFT，方括号外边的求和则实现结果数组在列上的DFT，这种分解使我们可以用一维FFT来快速实现二维DFT。
（2）二维快速傅立叶变换FFT（Fast Fourier Transform）
FFT的基本思想是：
令序列的长度为，如果不满足，在尾部补零，没有任何影响。按n奇偶把分解为两个N/2点的子序列：
（1-8）
（1-9）
那么
（1-10）
将（1-8）和（1-9）式代入上式得：
（1-11）
上式右边的两部分恰好是长度（周期）为N/2的的傅立叶变换，所以：
（1-12）
（1-13）
这样一个长度为N的DFT就分解为两个长度为N/2的DFT，然后进行N/2次复数的蝶式运算，再运用分解-递归思想，分解次，每一次均有N/2个蝶式运算，所以FFT的时间复杂度为。
FFT在图像运用中的思想是，先将原图像进行转置，按行对转置后的图像矩阵做一维FFT，将此变换所得的中间矩阵再转置，再按对转置后的中间矩阵做一维FFT，最后得到的就是二维FFT。

实验操作方法与步骤：
（1）嵌入算法的步骤如下：
1.将宿主图像分成8×8的块，每一块都进行DFT变换，得到与宿主图象相同尺寸的DFT域图象系数矩阵。
2. 产生两个不相关的伪随机序列，并设定一个密钥key。
3. 将每个8×8的DFT系数矩阵，从每一块的中频段取出个系数。
4. 由于DFT域的幅度谱具有对称性，为了水印嵌入后保持这种对称性，也为了确保恢复图像像素值为实数，嵌入水印时采用对称嵌入。嵌入规则为当水印图像元素为‘0’时，按照式子（其中alfa为尺度因子，mark为随机序列，为原图像的分块DFT系数）将一个随机序列与原始图像块的幅度谱对应元素进行乘性叠加；为‘1’时，用另一个伪随机序列与幅度谱对应元素进行乘性叠加。
5. 用得到的新的DFT系数对原来位置的DFT系数进行置换。
6. 对每一图像块进行DFT逆变换，得到含水印图像。
(2)提取算法步骤如下：
1. 将宿主图像分成8×8的块，每一块都进行DFT变换，得到与宿主图像相同尺寸的DFT域图象系数矩阵。
2. 对每一块的DFT系数进行“之”字型排列，从我们已知的嵌入位置取出可能含有水印的系数，将所有提取出的系数，按顺序组成新的系数序列。
3. 计算嵌入水印图像块的幅度谱与两个伪随机序列的相关性，假设嵌入的水印图像矩阵为‘0’时叠加的是随机序列A，为‘1’时叠加的是随机序列B，则提取时如果嵌入水印幅度谱与随机序列A的相关性大于其B的相关性，那么这个位置嵌入的是水印图像矩阵‘0’元素。按这个规则即可提取出水印图像。
4. 将含水印的系数序列与水印信号进行相关计算，得到归一化相关系数。
实验结果与分析:
以下是尺度因子alpha=260时的实验结果,以及受到各种攻击后所提取到的水印图像及其与原始水印图像的相关系数:

图1 原始图像与嵌入水印后的图像

图2 没有受到攻击的图像和直接提取的水印图像
此时相关系数=1

图3 加入白噪声后的图像、原始水印图像及提取的水印图像
此时相关系数=096566

图4 高斯低通滤波后的图像及提取的水印图像
此时相关系数=0.97311

图5 剪切后的图像及提取的水印图像
此时相关系数=0.96907

图6 旋转后的图像及提取的水印图像
此时相关系数=0.91301
结论：由于DFT域对全局性的旋转、平移和缩放变换具有不变性的特点，基于DFT域的数字水印技术要比基于DCT域的数字水印技术的鲁棒性要强的多，特别是对于旋转攻击，前者的优势更加突出！
程序参考：
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% 嵌入水印算法 %%%%%%%%%%%%%%
clear all;
clc;
%%%%%% 读取水印图象 %%%%%%%
omark=double(imread('muxiao.bmp'))/255;
mo=size(omark,1);
no=size(omark,2);
%%%%% 以下生成水印信息 %%%%%%
mark=omark;
alpha=200;
R1=randn(1,8);
rand1=round(R1);
R2=randn(1,8);
rand2=round(R2);
%%%%%%% 读取宿主图像 %%%%%%%
cimage=imread('lena.bmp');
figure(1);
subplot(2,3,1),imshow(cimage,[]),title('原始的宿主图像');
[mc,nc]=size(cimage);
cimage=double(cimage);
fft0=blkproc(cimage,[8,8],'fft2');
%%%%%%% 嵌入水印 %%%%%%%%%
fft_abs=abs(fft0);
fft_angle=angle(fft0);
fft1=fft_abs;
for ii=1:mo
for jj=1:no
x=(ii-1)*8;y=(jj-1)*8;
if mark(ii,jj)==1
k=rand1;
else
k=rand2;
end
fft1(x+1,y+8)=fft_abs(x+1,y+8)+alpha*k(1);
fft1(x+2,y+7)=fft_abs(x+2,y+7)+alpha*k(2);
fft1(x+3,y+6)=fft_abs(x+3,y+6)+alpha*k(3);
fft1(x+4,y+5)=fft_abs(x+4,y+5)+alpha*k(4);
fft1(x+5,y+4)=fft_abs(x+5,y+4)+alpha*k(5);
fft1(x+6,y+3)=fft_abs(x+6,y+3)+alpha*k(6);
fft1(x+7,y+2)=fft_abs(x+7,y+2)+alpha*k(7);
fft1(x+8,y+1)=fft_abs(x+8,y+1)+alpha*k(8);
end
end
fft_abs_0=fft1;
fft_abs_1=(fft_abs_0.*exp(i*fft_angle));
wimage_0=blkproc(fft_abs_1,[8,8],'ifft2');
wimage=abs(wimage_0);
%%%% 显示嵌入水印图象 %%%%%
wimage_1=uint8(wimage);
imwrite(wimage_1,'withmark.bmp','bmp');
subplot(2,3,2),imshow(wimage,[]),title('嵌入水印后的图象');
%%%%%%%%%%%%% 进行各种攻击 %%%%%%%%%%%%
disp('1-->加入白噪声');
disp('2-->高斯低通滤波');
disp('3-->剪切图像');
disp('4-->旋转图象');
disp('5-->直接检测');
begin=input('请选择攻击（1-5）');
switch begin
%%%%%%%%加入白噪声%%%%%%%%
case 1
Aimage1=wimage;
Wnoise=12*randn(size(Aimage1));
Aimage1=Aimage1+Wnoise;
subplot(2,3,4),imshow(Aimage1,[]),title('加入白噪声后的图象');
att=Aimage1;
imwrite(att,'whitenoiseimage.bmp');
%%%%%% 高斯低通滤波 %%%%%%%%
case 2
Aimage2=wimage;
H=fspecial('gaussian',[4,4],0.6);
Aimage2=imfilter(Aimage2,H);
subplot(2,3,4),imshow(Aimage2,[]),title('高斯低通滤波后的图象');
att=Aimage2;
imwrite(att,'gaussianimage.bmp')
%%%%%%% 剪切攻击 %%%%%%%%%
case 3
Aimage3=wimage;
Aimage3(1:128,1:128)=256;
subplot(2,3,4),imshow(Aimage3,[]),title('剪切后的图象');
att=Aimage3;
imwrite(att,'cutpartimage.bmp');
%%%%%%% 旋转攻击 %%%%%%%%%
case 4
Aimage4=wimage;
Aimage4=imrotate(Aimage4,10,'bilinear','crop');
Aimage_4=mat2gray(Aimage4);
subplot(2,3,4),imshow(Aimage_4,[]),title('旋转10度后的图象');
att=Aimage_4;
imwrite(att,'rotatedimage.bmp');
%%%%%%% 没有受到攻击 %%%%%%%
case 5
subplot(2,3,4),imshow(wimage,[]),title('直接提取的图象');
att=wimage;
imwrite(att,'directimage.bmp');
end
%%%%%%%%%%%%%%% 提取水印 %%%%%%%%%%%%%
tmark_0=att;
tmark_0=blkproc(tmark_0,[8,8],'fft2');
tmark_0=abs(tmark_0);
pass=zeros(1,8);
for ii=1:mo
for jj=1:no
x=(ii-1)*8;y=(jj-1)*8;
pass(1)=tmark_0(x+1,y+8);
pass(2)=tmark_0(x+2,y+7);
pass(3)=tmark_0(x+3,y+6);
pass(4)=tmark_0(x+4,y+5);
pass(5)=tmark_0(x+5,y+4);
pass(6)=tmark_0(x+6,y+3);
pass(7)=tmark_0(x+7,y+2);
pass(8)=tmark_0(x+8,y+1);
if (corr2(pa