1、实验一 图像输入与输出基本操作（2 学时） 实验二 基于 DCT 域的数字水印技术（2 学时） 综合性实验 基于傅立叶域变换的数字水印技术（4 学时）实验前预备知识：数字水印技术是利用数字产品普遍存在的冗余数据与随机性，将水印信息嵌入在数字产品本身中，从而起到保护数字产品版权或者完整性的一种技术。现在学术界对数字水印算法的理解都是将一些不易察觉的具有随机特性的数据嵌入到图像频域或空域的系数上。从信号处理的角度看，嵌入水印可以看成是在强背景下迭加一个弱信号，由于人类视觉系统的分辨率受到一定的限制，只要迭加的信号幅度不超过 HVS 的对比门限，人眼就无法感觉到信号的存在，所以可以通过对原始图像进行
2、一定调整，在不影响视觉效果的情况下嵌入一些水印信息。数字水印系统的一般模型如图 1 所示：图 1 数 字 水 印 系 统 基 本 模 型水印嵌入器的输入量有三个：水印信号 M，宿主信号 S 和密钥 K。水印信号 M 是指原始水印（图像或一个数字序列）通过一定的方法经过调制将嵌入到宿主信号中的数字信号。宿主信号 S 是指被嵌入水印的信号（原始信号） 。密钥 K 则指用于提高水印系统安全性的密码信息，它独立于宿主信号。密钥有私有密钥和公共密钥之分，前者指攻击者在明确了水印嵌入方法但又不知道密钥的情况下，水印不会被破坏或盗取；后者是指攻击者对宿主信号（如内容标识、语言字幕等）不感兴趣的情况下，密钥也
3、就不存在保密性，可以作为公共密钥。实验一 图像输入与输出基本操作一、实验题目：图像输入与输出操作二、实验目的学习在 MATLAB 环境下对图像文件的 I/O 操作，为读取各种格式的图像文件和后续进行图像处理打下基础。三、实验内容利用 MATLAB 为用户提供的专门函数从图像格式的文件中读 /写图像数据、显示图像，以及查询图像文件的信息。四、预备知识熟悉 MATLAB 开发环境。五、实验原理（1）图像文件的读取利用 imread 函数可以完成图像文件的读取操作 。常用语法格式为：I=imread(filename,fmt)或 I=imread(filename.fmt)；其作用是将文件名用字符串
4、 filename 表示的、扩展名用字符串 fmt（表示图像文件格式）表示的图像文件中的数据读到矩阵 I 中。当 filename 中不包含任何路径信息时，imread 会从当前工作目录中寻找并读取文件。要想读取指定路径中的图像，最简单的方法就是在 filename 中输入完整的或相对的地址。 MATLAB 支持多种图像文件格式的读、写和显示。因此参数 fmt 常用的可能值有：bmp Windows 位图格式jpgorjpeg 联合图像专家组格式tifortiff 标志图像文件格式gif 图形交换格式pcx Windows 画刷格式png 可移动网络图形格式xwd X Window Dump
5、格式例如，命令行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 和 TIF
6、F 图像。所以，当用户保存这类文件时，MATLAB 就将其存储在 uint16 中。（3）图像文件的显示图像的现实过程是将数字图像从一组离散数据还原为一幅可见图像的过程。MATLAB 的的图像处理工具箱提供了多种图像显示技术。例如 imshow 可以直接从文件显示多种图像；image 函数可以将矩阵作为图像 ；colorbar 函数可以用来显示颜色条；montage 函数可以动态显示图像序列。这里仅对常用的显示函数进行介绍。图像的显示imshow 函数是最常用的显示各种图像的函数，其调用格式如下： imshow(I,N)； imshow(I,N)用于显示灰度图像，其中 I 为灰度图像的数据矩阵
7、，N 为灰度级数目，默认值为 256。例如下面的语句用于显示一幅灰度图像： I=imread(lena.jpg)； imshow(I);如果不希望在显示图像之前装载图像，那么可以使用以下格式直接进行图像文件的显示：imshow filename其中，filename 为要显示的图像文件的文件名。实例 1-1 显示一幅在当前目录下的.bmp 格式的图像：imshow lena.jpg显示结果如图 1.1 所示。图 1.1.1 显示一幅图像文件中的图像注意事项：该文件名必须带有合法的扩展名（指明文件格式） ，且该图像文件必须保存在当前目录下，或在 MATLAB 默认的目录下。添加色带colorba
8、r 函数可以给一个坐标轴对象添加一条色带。如果该坐标轴对象包含一个图像对象，则添加的色带将指示出该图像中不同颜色的数据值。这对于了解被显示图像的灰度级特别有用。其调用格式为： colorbar 实例 1-2 I=imread(lena.jpg)； imshow(I); colorbar;图 1.1.2 显示图像并加入颜色条从上图可知，该图像是数据类型为 uint8 的灰度图像，其灰度级范围从 0255。显示多幅图像显示多幅图像最简单的方法就是在不同的图形窗口中显示它们。imshow 总是在当前窗口中显示一幅图像，如果用户想连续显示两幅图像，那么第二幅图像就会替代第一幅图像。为了避免图像在当前窗
9、口中的覆盖现象，在调用 imshow函数显示下一幅图像之前可以使用 figure 命令来创建一个新的窗口。例如：imshow(I1);figure, imshow(I2);figure, imshow(I3);有时为了便于在多幅图像之间进行比较，需要将这些图像显示在一个图形窗口中。达到这一目的有两种方法：一种方法是联合使用 imshow 和 subplot 函数，但此方法在一个图形窗口只能有一个调色板；另一种方法是联合使用 subimage和 subplot 函数，此方法可在一个图形窗口内使用多个调色板。subplot 函数将一个图形窗口划分为多个显示区域，其调用格式如下： subplot(m
10、,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 在一个图形窗口中显示多幅图
11、像（4） 图像文件信息的查询imfinfo 函数用于查询图像文件的有关信息，详细地显示出图像文件的各种属性。其语法格式为：info=imfinfo(filename,fmt)或 info=imfinfo(filename.fmt)或 imfinfo filename.fmtimfinfo 函数获取的图像文件信息依赖于文件类型的不同而不同，但至少应包含以下内容：文件名。如果该文件不在当前目录下，还包含该文件的完整路径。文件格式。文件格式的版本号。文件最后一次修改的时间。文件的大小。以字节为单位。图像的宽度。图像的高度。每个像素所用的比特数。也叫像素深度。图像类型。即该图像是真彩色图像、索引图像还
12、是灰度图像。例如，命令行imfinfo bubbles25.jpg会输出如下信息（注意，在这种情况下，有些域不包含信息）：Filename: bubbles25.jpgFileModDate: 04-Jan-2003 12:31:26FileSize: 13849Format: jpgFormatVersion: Width: 714Height: 682BitDepth: 8ColorType: grayscaleFormatSignature: Comment: 六、实验步骤：（1）利用 imread 函数完成对图像文件的读取操作。（2）利用 imwrite 函数完成图像的写入（保存）操作
13、。（3）利用 imshow 函数显示图像。（4）利用 imfinfo 函数查询图像文件的有关信息。实 验 二 基 于 DCT 域 的 数 字 水 印 技 术离散余弦变换，简称 DCT，是一种实数域变换，其变换核为实数的余弦函数，计算速度较快，是一种近似最佳变换，很适合于做图象压缩和随机信号处理，它对常用的图象压缩有较强的鲁棒性。基于 DCT 域的数字水印算法，一般而言有两种算法。一种是直接对整幅图像进行 DCT 计算，选取合适频段的系数，嵌入水印。另一种是，首先将整幅图像分成块，对每一块分别进行 DCT 计算，在每一块中选取合适频段的系数，将水印信息分散嵌入到每一块所选取的 DCT 系数中。这
14、种方法称为分块 DCT，这种方法非常适合于 88 的图像块 DCT 变换。目前，DCT 域的水印方法大多数是采用分块 DCT 方法。二维 DCT 概述：一般而言，数字图象 ,其二维NnMmnS0,),(DCT 变换用矩阵 表示，其定义式如下式所示：S（4-1）)21(cos)21(cos),()(2),( 210 NljMkijiSlkNMlkSij 二维 IDCT 定义 如下： 8（4-2）)21(cos)21(cos),()(2),( 210 ljkilkSjiSkl其中 并且,.1,., Nlj（4-3）,.21,)(1Mkc1,.2,10)(2Nl我们最常用到的是 88 分块 DCT，
15、也就是首先将大小为 MN 的图像分成的非重叠的大小为 88 的块，然后对每一块作64/)()/(8/NMDCT。根据式（4-1）和式（4-2） ，88DCT 的定义如下 ：（4-4）)162cos()162(cos),(2)(),(701 ljkijiSlklSijc88IDCT 的定义如下：（4-5）)162(cos)162(cos),(2)(),(701 ljkilkSjiSkl其中， 并且,.,（4-6）7,.21,)(1kc7,.21,0)(2l实验目的：学习使用基于 DCT 域数字水印嵌入算法和基于 DCT 域数字水印检测提取算法。实验内容和原理：（1）基于 DCT 域数字水印嵌入算
16、法：DCT 域的水印嵌入过程，就是首先对宿主图像 S 进行分块 DCT 运算，选取出要嵌入的 DCT 系数 C，并且对宿主图像用 HVS 模型进行分析，得到感知系数掩蔽模板，用其来控制水印嵌入强度，再将编码后的水印信号 W，用一定方式嵌入到选好的系数中。用嵌入水印的系数替代原来位置的系数，再进行逆DCT（Inverse Discrete Cosine Transform,IDCT）变换，就得到了含有水印的图像X。DCT 域水印系统的水印嵌入过程框图如图 1 所示：图 1 DCT 域水印嵌入过程方框图（2）基于 DCT 域数字水印检测提取算法：DCT 域水印的验证过程（即包括水印的检测和提取过程
17、） ，就是将待测图像Y（注意此时的图像也许已经遭受了图象处理或恶意攻击）进行 DCT 运算；然后按着嵌入水印时选取的 DCT 系数，选取其含有水印的系数，进行水印提取，然后利用相关检测法判断水印是否存在。如果水印检测器输出结果显示水印存在，则根据需要可以用水印解码器对提取的水印进行解码，恢复水印。水印检测提取过程如图 2图 2 DCT 域数字水印检测提取框图实验操作方法与步骤：（1）水印嵌入步骤如下：1. 将宿主图像 分成 88 的块，每一块都进),.21;,.21(), NjMijS行 DCT 变换，得到与宿主图象相同尺寸的 DCT 域系数矩阵。),.;,.21(),NjMijS2. 产生两
18、个不相关的伪随机序列，并设定一个密钥 key。3. 将每个 88 的 DCT 系数矩阵，从每一块的中频段取出个系数 。)/()864(NM),.21(),numKBkjiC4. 然后嵌入水印，嵌入规则为当水印图像元素为 0时， 按照式子（其中 alpha 为尺度因子，mark 为随机序列， 为原图)*1( marklphW W像的分块 DCT 系数） ，将一个随机序列与原始图像块的幅度谱对应元素进行乘性叠加；为1时，用另一个伪随机序列与幅度谱对应元素进行乘性叠加。5. 用得到的新的 DCT 系数对原来位置的 DCT 系数进行置换。6. 对新的 DCT 系数矩阵进行 DCT 逆变换，得到了嵌入水
19、印信号的图像。（ 2） 水印检测提取步骤如下：1. 对含有水印图像进行 88 的块 DCT 变换。2. 对每一块的 DCT 系数进行“之”字型排列，从我们已知的嵌入位置取出可能含有水印的系数 ， 将所有提取出的系数，按顺序组成新的系数序列),(jiCk。 ),.21(),wNlC3. 用密钥生成水印信号 )1,0(NW4. 将含水印的系数序列与水印信号进行相关计算，得到归一化相关系数 。NC实验结果及分析：以下是尺度因子 alpha=30 时的实验结果,以及受到各种攻击后所提取到的水印图像及其与原始水印图像的相关系数 :NC图 3 原始图像和嵌入水印后的图像此时，alpha=30图 4 没 有
20、 受 到 攻 击 的 图像及 提 取 到 的 水 印 图像此时相关系数 =1NC图 5 加入白噪声后的图像、原始水印和提取的水印图像此时相关系数 =0.99969NC图 6 高 斯 低 通 滤 波 后 的 图像及 提 取 的 水 印 图像此时相关系数 =0.99058NC图 7 剪 切 后 的 图像及 提 取 的 水 印 图像此时相关系数 =0.93909NC图 8 旋 转 10 后 的 图像及 提 取 的 水 印 图像此时相关系数 =0.48383NC结论:当图像 (宿主信息)受到加噪、高斯低通滤波、剪切等攻击后，其可视度、提取的水印图像与原始水印图像的相似度（都在 90%以上）都是非常令人
21、满意的；当受到旋转攻击后，其提取的水印图像与原始水印图像的相似度（不足 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),imsho
22、w(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)=1k=rand1;elsek=rand2; endcda1(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
23、,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);endend% 嵌入水印后图象 %wimage=blkproc(cda1,8,8,idct2); wimage_1=uint8(wimage);imwrite(wimage_1,withmark.bmp,bmp);subplo
24、t(2,3,2),imshow(wimage,),title(嵌入水印后图象 );% 进行攻击测试 %disp(1加入白噪声 );disp(2高斯低通滤波 );disp(3剪切图像 );disp(4旋转攻击 );disp(5直接检测 )begin=input(请选择攻击（ 1-5） ： )switch begin% 加入白噪声 %case 1Aimage1=wimage;Wnoise=20*randn(size(Aimage1);Aimage1=Aimage1+Wnoise;subplot(2,3,4),imshow(Aimage1,),title(加入白噪声后的图象 );att=Aimage
25、1;imwrite(att,whitenoiseimage.bmp);% 高斯低通滤波 %case 2Aimage2=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 3Aimage3=wimage;Aimage3(1:128,1:128)=256;subplot(2,3,4),imshow(Aimage3,
26、),title(剪切后的图象 );att=Aimage3;imwrite(att,cutpartimage.bmp);% 旋转攻击 %case 4Aimage4=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 5subplot(2,3,4),imshow(wima
27、ge,),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:mofor j=1:nox=(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);p
28、ass(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;elsetmark_1(i,j)=0;endendend% 计算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:mfor j=1:nnc_0=nc_0+g_mark(i,j)*o_mar
29、k(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);endendNC=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);综合设计性实验 基于傅立叶域变换的数字水印技术实验目的：学习使用基于傅立叶域的数字水印嵌入算法和
30、基于 DCT 域数字水印检测提取算法。实验内容和原理：傅立叶变换（Fourier Transform）是线性系统分析的一个有力工具，是研究信号的频谱方法，它架起了时域和频域之间的桥梁。使我们能够定量分析诸如数字化系统、采样点、电子放大器、卷积滤波器、噪声等的作用，把傅立叶变换的理论同其物理解释相结合，将大大有助于解决大多数图象处理问题。Fourier 分析理论十分完善，既可以处理连续信号也可以处理离散信号。 计算机只能处理离散信号， 且在数字图象处理中， 输 入 图 像 和 输 出 图 像 通 常 都是二维的， 一般表示成二维数字矩阵， 因此这里直接讨论二维离散傅立叶变换（DFT ）和二维快速
31、傅立叶变换（FFT ） 。（ 1） 二 维 离 散 傅 立 叶 变 换 DFT（ Discrete Fourier Transform）二维离散傅立叶变换对定义式如下：（ 1-1）,210;,210,),(),(10)(2 NvMueyxfvuFNyMNvyBxj （ 1-2）10)(2 ,;,),),(vuvxMiFxf 二维离散傅立叶变换的傅立叶谱、相位、功率谱分别如下：傅立叶谱： （ 1-2122),(),(),(vuIRvu3）相位： （ 1-),(,),(Iarctg4）功率谱： （ 1-),(),(),(),( 22vuIRvuFE5）式子（1-1）可以分离为：（ 1-6）1010
32、22),(),(NyMxuxMjvyj eyfevuF式子（1-2）可以分离为：（ 1-7）101022),(),(NvMuuxjvyj evFeyxf 由上两式可知，二维傅立叶变换可由连续两次一维傅立叶变换得到 ， 从而将二维 DFT 分解为水平和垂直两部分运算，上式中方括号中的项表示 在 图像的行上计算的 DFT，方括号外边的求和则实现结果数组在列上的 DFT，这种分解使我们可以用一维 FFT 来快速实现二维 DFT。（2）二维快速傅立叶变换 FFT（Fast Fourier Transform）FFT 的基本思想是：令序列 的长度 为 ，如果不满足， 在尾部补零，没有任何影响。)(nfN
33、M2按 n 奇偶把 分解为两个 N/2 点的子序列：（1-8）);2/(,10),()1 Nmfg（1-9）2(h那么（1-10）为 奇 数为 偶 数 nknNjnkNj efefkF22)()()(将（1-8）和（1-9）式代入上式得：（1-11）2211200()()()NNjkmjkmNmFkgehe上式右边的两部分恰好是长度（周期）为 N/2 的 的傅立叶变换)(,1hg，所以：)(,1kHG（1-12）)2(,10),()(121 kHekGFkNj （1-13）)1(,),()(2(121Nkj 这样一个长度为 N 的 DFT 就分解为两个长度为 N/2 的 DFT，然后进行 N/
34、2次复数的蝶式运算，再运用分解-递归思想，分解 次 ，每一次均有 N/2Mlog个蝶式运算，所以 FFT 的时间复杂度为 。)l(OFFT 在图像运用中的思想是，先将原图 像 进行转置 ， 按行对转置后的图像矩阵做一维 FFT， 将此变换所得的中间矩阵再转置，再按对转置后的中间矩阵做一维 FFT，最后得到的就是二维 FFT。实验操作方法与步骤：（1）嵌入算法的步骤如下：1.将宿主图像 分成 88 的块，每一块都进行),.21;,.21(), NjMijfDFT 变换，得到与宿主图象相同尺寸的 DFT 域图象系数矩阵 。(,)kFij2. 产生两个不相关的伪随机序列，并设定一个密钥 key。3.
35、 将每个 88 的 DFT 系数矩阵，从每一块的中频段取出个系数 。)/()864(NMnumKBkjiC,.21),(4. 由于 DFT 域的幅度谱具有对称性，为了水印嵌入后保持这种对称性 ，也为了确保恢复图像像素值为实数，嵌入水印时采用对称嵌入。嵌入规则为当水印图像元素为 0时，按照式子 （其中 alfa 为尺度因子，)*1( arklfWmark 为随机序列， 为原图像的分块 DFT 系数） ， 将一个随机序列与原始图像块的幅度谱对应元素进行乘性叠加；为1时，用另一个伪随机序列与幅度谱对应元素进行乘性叠加。5. 用得到的新的 DFT 系数对原来位置的 DFT 系数进行置换。6. 对每一图
36、像块进行 DFT 逆变换，得到含水印图像 。(2)提取算法步骤如下：1. 将宿主图像 分成 88 的块，每一块都进行NjMijf ,.21;,.21),(DFT 变换，得到与宿主图像相同尺寸的 DFT 域图象系数矩阵。NjMijFk,.;,.1),(2. 对每一块的 DFT 系数进行“之”字型排列，从我们已知的嵌入位置取出可能含有水印的系数 ，将所有提取出的系数，按顺序组成新的系数序列),(jiCk。 wlC,.21),(3. 计算嵌入水印图像块的幅度谱与两个伪随机序列的相关性，假设嵌入的水印图像矩阵为0时叠加的是随机序列 A，为 1时叠加的是随机序列 B，则提取时如果嵌入水印幅度谱与随机序列
37、 A 的相关性大于其 B 的相关性，那么这个位置嵌入的是水印图像矩阵0元素。按这个规则即可提取出水印图像。4. 将含水印的系数序列与水印信号进行相关计算，得到归一化相关系数 。NC实验结果与分析:以下是尺度因子 alpha=260 时的实验结果,以及受到各种攻击后所提取到的水印图像及其与原始水印图像的相关系数 :NC图 1 原始图像与嵌入水印后的图像图 2 没有受到攻击的图像和直接提取的水印图像此时相关系数 =1NC图 3 加入白噪声后的图像、原始水印图像及提取的水印图像此时相关系数 =0.96566NC图 4 高斯低通滤波后的图像及提取的水印图像此时相关系数 =0.97311NC图 5 剪切
38、后的图像及提取的水印图像此时相关系数 =0.96907NC图 6 旋转后的图像及提取的水印图像此时相关系数 =0.91301NC结论：由于 DFT 域对全局性的旋转、平移和缩放变换具有不变性的特点，基于 DFT 域的数字水印技术要比基于 DCT 域的数字水印技术的鲁棒性要强的多，特别是对于旋转攻击，前者的优势更加突出！程 序 参 考 ：% 嵌入水印算法 %clear all; clc;% 读取水印图象 % omark=double(imread(muxiao.bmp)/255;mo=size(omark,1);no=size(omark,2);% 以下生成水印信息 %mark=omark;al
39、pha=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:
40、mo for jj=1:no x=(ii-1)*8;y=(jj-1)*8;if mark(ii,jj)=1k=rand1;elsek=rand2; endfft1(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)
41、=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);endendfft_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
42、);subplot(2,3,2),imshow(wimage,),title(嵌入水印后的图象) ;% 进行各种攻击 %disp(1加入白噪声) ;disp(2高斯低通滤波 );disp(3剪切图像 );disp(4旋转图象 );disp(5直接检测 );begin=input(请选择攻击（ 1-5） );switch begin%加入白噪声% %case 1Aimage1=wimage;Wnoise=12*randn(size(Aimage1);Aimage1=Aimage1+Wnoise;subplot(2,3,4),imshow(Aimage1,),title(加入白噪声后的图象 );a
43、tt=Aimage1;imwrite(att,whitenoiseimage.bmp);% 高斯低通滤波 %case 2Aimage2=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 3Aimage3=wimage;Aimage3(1:128,1:128)=256;subplot(2,3,4),imshow
44、(Aimage3,),title(剪切后的图象 );att=Aimage3;imwrite(att,cutpartimage.bmp);% 旋转攻击 %case 4Aimage4=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 5subplot(2,3,4),imshow(w
45、image,),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:mofor jj=1:nox=(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
46、);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(ii,jj)=1;elsetmark_1(ii,jj)=0;endendend% 计算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 ii=1:mfor jj=1:nnc_0=nc_0+g_mark(ii,jj)*o_mark(ii,jj);nc_1=nc_1+o_mark(ii,jj)*o_mark(ii,jj);nc_2=nc_2+g_mark(ii,jj)*g_mark(ii,jj);endendNC=nc_0/sqrt(nc_1*nc_2);% 显示提取水印图象
