一、引言
该篇包括三部分，1）引言、2）图像变化技术简介和代码实现 、3）基于图像变换技术的数字水印技术及代码实现。
数字水印是一种有效的数字产品版权保护和数据安全维护技术, 是信息隐藏领域的一个重要分支, 也是密码学的一种有益的补充技术。近年来它引起了人们的广泛关注。图像隐形水印是其中最主要的研究方向, 按照嵌入位置可分为空间域方法和变换域( DCT, DFT 和DWT 等) 方法。为了使水印信息有更好的鲁棒性（抵抗攻击的能力）众多学者采用变换域中的嵌入方法，前面我有两篇博客涉及到了有关DCT的数字水印方法和算法，接下来将分别介绍基于DCT、DFT、DWT三种基础图像变换技术的数字水印方法及代码实现。
二、DFT、DCT、DWT图像变换技术
图像变换技术是为了能让图像用正交函数或正交矩阵表示而对原图像所作的变换，该变换是二维线性可逆的。一般称原始图像为空间域图像，称变换后的图像为转换域图像（也称为频率域），转换域图像可反变换为空间域图像。经过图像变换后，一方面能够更有效地反映图像自身的特征，另一方面也可使能量集中在少量数据上，更有利于图像的存储、传输及处理。从数学角度来说，图像变换是把图像中的像素表达成“另外一种形式”满足实际需求，大多数情况下需要对变换处理后的图.像进行逆变换，从而获取处理以后的图像。这种变换同样可以应用于其他有关物理或数学的各种问题中，并可以采用其他形式的变量。在数字水印领域中，常用的基础图像变换技术有傅里叶变换（DFT）、离散余弦变换（DCT）、小波变换（DWT）。
2.1 图像傅里叶变换（DFT）
傅里叶变换是以时间为自变量的“信号”与频域为自变量的“频普”函数之间的某种变换关系。从纯粹的数学意义上看，傅立叶变换是将一个图像函数转换为一系列周期函数来处理的;从物理效果看，傅立叶变换是将图像从空间域转换到频率域，其逆变换是将图像从频率域转换到空间域。实际上对图像进行二维傅立叶变换得到频谱图，就是图像梯度的分布图，傅立叶频谱图上看到的明暗不一的亮点，实际上图像上某--点与邻域点差异的强弱，即梯度的大小，也即该点的频率大小。如果频谱图中暗的点数更多，那么实际图像是比较柔和的:反之，如果频谱图中亮的点数多，那么实际图像一定是尖锐的，边界分明且边界两边像素差异较大。
傅里叶变换包括一维连续傅里叶变换、二维连续傅里叶变换、一维离散傅里叶变换、二维离散傅里叶变换。在图像处理领域中常用二维离散傅里叶变换。其定义如下式(1)：
(1)
其中u = 0, 1, ... , M-1; v = 0, 1, ... , N-1 。 F（u, v）的二维离散傅里叶反变换定义为下式（2）：
(2)
图像的二维傅里叶结果通常是复数形式，故在进行图像处理的过程将结果进行平方处理，其具体的数学计算公式以及其他三种傅里叶变化的定义式可查阅《数字图像处理》——冈萨雷斯。
定义式是抽象的表达了傅里叶的变换过程，但是落实图像是怎样的呢，接下来我们用MATLAB工具进行图像的傅里叶变换。其代码及结果如下：
close all;clear all;clc;I = imread('pepper.bmp');J = rgb2gray(I);K_1 = fft2(J);L_1 = abs(K_1/256);K_2 = fftshift(K_1);L_2 = abs(K_2/256);figure;subplot(2,2,1);imshow(I); title('原始图像');subplot(2,2,2);imshow(J); title('灰度图像');subplot(2,2,3);imshow(uint8(L_1)); title('灰度图像傅里叶频谱');subplot(2,2,4);imshow(uint8(L_2)); title('平移后的频谱');
其运行结果如下图：
2.2 离散余弦变换（DCT）
离散余弦变换(DCT)是一组不同频率和幅值的余弦函数和来近似一副图像，实际上是傅里叶变换的实数部分。由于离散余弦变量对于一副图像，其大部分可视化信息都集中在少数的变换系数上。因此，离散余弦变量是数据压缩常用的一个变换编码方法，它能将高相关数据能量集中，使得它非常适用于图像压缩，例如国际压缩标准的JPEG格式中就采用了离散余弦变换。
在傅立叶变换过程中，如果被展开的函数是实偶函数，那么其傅立叶变换中只包含余弦项，基于傅立叶变换的这一特点，人们提出了离散余弦变换。DCT变换先将图像函数变换成偶函数形式，再对其进行二维离散傅立叶变换，因此DCT变换可以看成是一种简化的傅立叶变换。离散余弦变换叶分为一维离散余弦变换和二维离散余弦变换，在图像处理中用到二维变化，所以介绍下其二维定义，其他定义可查阅《数字图像处理》——冈萨雷斯。其二维离散余弦变换定义如下式：
(3)
其中，
二维离散余弦反变换定义式：
（4）
用MATLAB工具进行图像的傅里叶变换。其代码及结果如下：
close all;clear all;clc;I = imread ('boats.bmp');J = dct2(I);figure('Name','离散余弦变换');subplot(1,2,1);imshow(I); title('原始图像');subplot(1,2,2);imshow(log(abs(J))); title('离散余弦变换系数图像');
其运行结果如下图2：
2.3 小波变换（DWT）
小波变化是对傅里叶变换和短时傅里叶变换的一个突破，其改变就在于，将无限长的三角函数基换成了有限长的会衰减的小波。小波变化的理论较多，更多内容可查阅《数字图像处理》——冈萨雷斯相关章节内容，以下从一个叫简单的角度来解释小波变换。
由上图也可以看出为什么叫小波，其波形较小。
其一维连续小波的定义式为：
(5)
连续小波变换具有如下的重要性质： 线性、平移不变性、 冗余性
在 MATLAB 小波变换工具箱中提供了多个母小波族函数如下表所示，不同的母小波具有不同的特点，其对图像的处理结果也不同。
接下来以“haar”母小波为例，进行处理例，其代码如下：
close all;clear all;clc;I = imread('lena.jpg');I = rgb2gray(I);I = im2double(I);[ca1,ch1,cv1,cd1]=dwt2(I,'haar');[ca2,ch2,cv2,cd2]=dwt2(ca1,'haar');figure(1);imshow(I);figure('Name','载体小波分解')% subplot(121)% imshow([ca1,ch1;cv1,cd1])% title('一级小波分解')% subplot(122)% imshow([ca2,ch2;cv2,cd2])% title('二级小波分解')%% 或者subplot(121)imagesc([wcodemat(ca1),wcodemat(ch1);wcodemat(cv1),wcodemat(cd1)])title('一级小波分解')subplot(122)imagesc([wcodemat(ca2),wcodemat(ch2);wcodemat(cv2),wcodemat(cd2)])
其运行结果如下图：
三、 基于DCT、DFT、DWT图像变换的数字水印技术代码实现
引言部分提到数字水印是密码学的一个重要分支之一。其实现主要分为三部分，水印标识的嵌入、图像攻击、水印提取。其过程可由下图5-7表示：
鉴于对文章篇幅的控制，本文对水印的嵌入和提取两部分进行实现，图像的攻击在后续进行更新总结。
数字水印的嵌入与提取算法有多种多样，那怎样评价数字水印嵌入效果和提取效果的好坏呢？业界是这样评价的：针对嵌入效果，我们用峰值信噪比（PSNR）进行评价，数值越大嵌入效果越好。针对水印的提取效果，我们用归一化相关系数（NC）进行评价，数字越大效果越好。那么PSNR和NC的MATLAB代码怎样实现呢，见下代码
PSNR函数MATLAB代码：
%nc(归一化相关系数)function dNC = nc(ImageA,ImageB)if (size(ImageA,1) ~= size(ImageB,1)) or (size(ImageA,2) ~= size(ImageB,2))error('ImageA <> ImageB');dNC = 0;return ;endImageA=double(ImageA);ImageB=double(ImageB);M = size(ImageA,1);N = size(ImageA,2);d1=0 ;d2=0;d3=0;for i = 1:Mfor j = 1:Nd1=d1+ImageA(i,j)*ImageB(i,j) ;d2=d2+ImageA(i,j)*ImageA(i,j) ;d3=d3+ImageB(i,j)*ImageB(i,j) ;endenddNC=d1/(sqrt(d2)*sqrt(d3));return
NC函数MATLAB 实现代码如下：
% PSNR (峰值信噪比)function dPSNR = psnr(ImageA,ImageB)if (size(ImageA,1) ~= size(ImageB,1)) or (size(ImageA,2) ~= size(ImageB,2))error('ImageA <> ImageB');dPSNR = 0;return ;endImageA=double(ImageA);ImageB=double(ImageB);M = size(ImageA,1);N = size(ImageA,2);d = 0 ;for i = 1:Mfor j = 1:Nd = d + (ImageA(i,j) - ImageB(i,j)).^2 ;endenddPSNR = 10*log10((M*N*max(max(ImageA.^2)))/d) ;return
在后面的代码中将直接调用这两个函数，不再重复编写。
3.1 基于DCT变换的数字水印嵌入与提取流程的代码实现
基于DCT数字水印代码叫简单，直接上代码，其具体需要注意的内容在代码注释中有具体讲解。
% The watermark image is embedded into the color image by DCT transformation % 将水印图像做DCT变换嵌入到彩色图像中 clc;clear all;wtype = 'dct2';iwtype = 'idct2';originalImage=imread('lena256.bmp');alpha=0.1; % Embedded factor 嵌入因子 dim_i=size(originalImage);rm=dim_i(1);cm=dim_i(2);% Embedding formula 设嵌入公式v‘=v（1+aXk）中的a=0.1% Display carrier image显示原始图片subplot(2,2,1);imshow(originalImage);title('Original Image');% Read and display watermark images 水印图片watermark=imread('mark32.bmp');watermark = rgb2gray(watermark);subplot(2,2,2);imshow(watermark);title('Watermark Image');% DCT transform of watermark image对水印图片做DCT变换watermark_dct=blkproc(watermark,[8,8],wtype);waterseria=watermark_dct(:); % 列向量% Extract green channel提取绿色分量image=imread('lena256.bmp'); % 重新读取一次原始图像，为了分离和复原颜色通道更方便。image_g=image(:,:,2);% DCT transformed for green channel 对绿色分量进行DCT变换dct_image_g=blkproc(image_g,[8,8],wtype); dct_image1=dct_image_g;k=1;for i=1:rm/8for j=1:cm/8x=(i-1)*8;y=(j-1)*8;ave=(dct_image_g(x+3,y+5)+dct_image_g(x+5,y+5)+dct_image_g(x+4,y+4)+dct_image_g(x+4,y+6)+dct_image_g(x+2,y+7))/5;dct_image1(x+4,y+5)=ave+alpha*waterseria(k);k=k+1;endendimage2=blkproc(dct_image1,[8,8],iwtype); % 反DCT变换%将嵌入水印的绿色分量加到原图中image(:,:,2)=image2;% 计算PSNRPSNR=psnr(image_g,image2);image = uint8(image);%加水印之后的图片subplot(2,2,3);imshow(image)name='嵌入水印图像';title(strcat(num2str(name),' k=',num2str(alpha),' PSNR=',num2str(PSNR)));imwrite(image,'withmark.bmp','bmp');%% 水印提取%提取绿色分量outpicture = imread('withmark.bmp');out_image=outpicture(:,:,2); %提取绿色分量%水印提取过程outdct_image=blkproc(out_image,[8,8],wtype); %dct变换 %提取频谱k=1;for i=1:rm/8for j=1:cm/8x=(i-1)*8;y=(j-1)*8;ave=(outdct_image(x+3,y+5)+outdct_image(x+5,y+5)+outdct_image(x+4,y+4)+outdct_image(x+4,y+6)+outdct_image(x+2,y+7))/5;outwaterseria(k)=(outdct_image(x+4,y+5)-ave)/alpha;k=k+1;endendk=1;for i=1:rm/8for j=1:cm/8outwatermark_dct(j,i)=outwaterseria(k);k=k+1;end endoutwatermark=blkproc(outwatermark_dct,[8,8],iwtype);outwatermark = uint8(outwatermark);%% 计算NC（归一化相关系数）%%NC=nc(outwatermark,watermark)subplot(2,2,4);imshow(outwatermark);name = 'Extract image';title(strcat(num2str(name), ' NC=',num2str(NC)));
其运行结果如下图：
3.2基于DFT变换的数字水印嵌入与提取流程的代码实现
该代码中用到了图像置乱，考虑到置乱周期问题 ，所设计的代码要求有两：1）水印图像的尺寸为40*40dpi；2）原始载体图像的 （长*宽）/ （分块的平方）要大于等于40*40. 以下为代码。由于代码较长，遵循代码模块化设计原则，将代码分为嵌入embed和提取recover两部分，其具体需要注意的内容在代码注释中有具体讲解，如下：
embed 嵌入：
%基于傅立叶域的数字水印%注意:水印必须为40*40的二值图像 %因为40阶的二维arnold置乱周期为30,所以嵌入时置乱8次,提取时置乱22.可以根据自己的需要更改.%嵌入源码clcclear all;% 保存开始时间start_time=cputime;iTimes=8;%置乱次数k=1.5; % 设置嵌入强度系数blocksize=8;% 块的大小filter_m=[ 1,1,1,1,1,1,1,1; % 滤波矩阵1,1,1,1,1,1,1,1;1,1,0,0,0,0,1,1;1,1,0,0,0,0,1,1;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0;];% 读入载体图像original_image=imread('lena.jpg');file_name_g=original_image(:,:,2);%提取绿色分量cover_object=double(file_name_g)/255;% 载体图像矩阵的行数与列数Mc=size(cover_object,1); Nc=size(cover_object,2); % 最大嵌入信息量max_message=Mc*Nc/(blocksize^2);% 读入水印图像mark=imread('mark40.bmp');mark=rgb2gray(mark);mark=im2bw(mark);message=double(mark);%水印图像矩阵的行数与列数Mm=size(message,1);Nm=size(message,2);% 检查水印信息是否过大if Mm*Nm>max_message % 载体图像的行*列要大于102400 error('水印信息过大')end%对水印图像进行Arnold置乱if Mm~=Nmerror('水印矩阵必须为方阵');endif Mm~=40error('必须为40*40大小,或者修改置乱次数');endtempImg=message;for n=1:iTimes % 次数for u=1:Mmfor v=1:Nmtemp=tempImg(u,v);ax=mod(u+v,Mm)+1;ay=mod(u+2*v,Nm)+1;outImg(ax,ay)=temp;endendtempImg=outImg;endmessage_vector=reshape(outImg,1,Mm*Nm);% 将cover_object（原图绿色通道矩阵）写入watermarked_imagewithmark_image=cover_object;%置随机数发生器的状态为1100key=1100;rand('state',key);% 产生伪随机序列pn_sequence_zero=round(2*(rand(1,sum(sum(filter_m)))-0.5));pn_sequence_one=round(2*(rand(1,sum(sum(filter_m)))-0.5)); % 将图像分块x=1;y=1;h=waitbar(0,'嵌入水印，请等待');for (kk = 1:length(message_vector))% 做傅立叶变换fft_block=fft2(cover_object(y:y+blocksize-1,x:x+blocksize-1));%计算幅值abs_block=fftshift(abs(fft_block));%计算相位angle_block=angle(fft_block);% 当message_vector=0且filter_m=1时用伪随机序列pn_sequence_zero叠加abs_block% 当message_vector=1且filter_m=1时用伪随机序列pn_sequence_one叠加abs_blockll=1;if (message_vector(kk)==0)for ii=1:blocksizefor jj=1:blocksizeif (filter_m(ii,jj)==1)abs_block_o=abs_block(ii,jj);abs_block(ii,jj)=abs_block(ii,jj)*(1+k*pn_sequence_zero(ll));abs_block(blocksize-ii+1,blocksize-jj+1)=abs_block(blocksize-ii+1,blocksize-jj+1)+abs_block(ii,jj)-abs_block_o;ll=ll+1;endendendelse for ii=1:blocksizefor jj=1:blocksizeif (filter_m(ii,jj)==1)abs_block_o=abs_block(ii,jj);abs_block(ii,jj)=abs_block(ii,jj)*(1+k*pn_sequence_one(ll));abs_block(blocksize-ii+1,blocksize-jj+1)=abs_block(blocksize-ii+1,blocksize-jj+1)+abs_block(ii,jj)-abs_block_o;ll=ll+1;endendendend% 进行傅立叶逆变换abs_block=fftshift(abs_block);withmark_image(y:y+blocksize-1,x:x+blocksize-1)=abs(ifft2(abs_block.*exp(i*angle_block))); % 移动到下一块if (x+blocksize) >= Ncx=1;y=y+blocksize;elsex=x+blocksize;endwaitbar(kk/length(message_vector),h);endclose(h);% 转换为uint8，并写入watermarked_image_g=(withmark_image*255);original_image(:,:,2)=watermarked_image_g;imwrite(original_image,'withwatermarked_image.bmp','bmp');%计算运行时间elapsed_time=cputime-start_time%计算psnrPSNR=psnr(cover_object,withmark_image);% 显示水印，嵌入水印图像与原始图像figure(1)subplot(2,2,2)imshow(message);title('原始水印');subplot(2,2,3);imshow(tempImg);title('置乱水印');subplot(2,2,4)imshow(original_image);name='嵌入水印图像';title(strcat(num2str(name),' k=',num2str(k),' PSNR=',num2str(PSNR)));subplot(2,2,1)imshow(original_image)title('原始图像');
其运行结果如下入：
recover 提取代码：
%提取clcclear all;% 保存开始时间start_time=cputime;iTimes=22; %置乱次数blocksize=8; % 设置块的大小filter_m=[ 1,1,1,1,1,1,1,1; % 滤波矩阵1,1,1,1,1,1,1,1;1,1,0,0,0,0,1,1;1,1,0,0,0,0,1,1;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0;];% 读入嵌入水印图像withmark_1=imread('withwatermarked_image.bmp');withmark=withmark_1(:,:,2);imshow(withmark_1)watermarked_image=double(withmark)/255;% 嵌入水印图像矩阵的行数与列数Mw=size(watermarked_image,1);Nw=size(watermarked_image,2); % 最大可嵌入信息量max_message=Mw*Nw/(blocksize^2);% 读入原始水印mark=imread('mark40.bmp');mark=rgb2gray(mark);orig_watermark=double(mark);% 原始水印矩阵的行数与列数Mo=size(orig_watermark,1);No=size(orig_watermark,2);%置随机数发生器的状态为1100key=1100;rand('state',key);% 产生伪随机序列pn_sequence_zero=round(2*(rand(1,sum(sum(filter_m)))-0.5));pn_sequence_one=round(2*(rand(1,sum(sum(filter_m)))-0.5)); % 将图像分块x=1;y=1;h=waitbar(0,'提取水印，请等待');for (kk = 1:max_message)% 傅立叶变换fft_block_w=fft2(watermarked_image(y:y+blocksize-1,x:x+blocksize-1));abs_block_w=abs(fftshift(fft_block_w));ll=1;for ii=1:blocksizefor jj=1:blocksizeif (filter_m(ii,jj)==1)sequence(ll)=abs_block_w(ii,jj); ll=ll+1;endendend% 计算sequence与pn_sequence_zero和pn_sequence_one的相关系数correlation_zero(kk)=corr2(pn_sequence_zero,sequence);correlation_one(kk)=corr2(pn_sequence_one,sequence);% 移动到下一块if (x+blocksize) >= Nwx=1;y=y+blocksize;elsex=x+blocksize;endwaitbar(kk/max_message,h);endclose(h);% 如果correlation_zero>correlation_one，那么message_vector=0，反之为1.for (kk=1:Mo*No)if correlation_zero(kk)>correlation_one(kk)message_vector(kk)=0;elsemessage_vector(kk)=1;endend% Arnold置乱tempImg=reshape(message_vector(1:Mo*No),Mo,No);message_arnold=tempImg;for n=1:iTimes % 次数for u=1:Mofor v=1:Notemp=tempImg(u,v);ax=mod(u+v,Mo)+1;ay=mod(u+2*v,No)+1;outImg(ax,ay)=temp;endendtempImg=outImg;endmessage=outImg;message=uint8(message);imwrite(outImg,'marked.jpg','jpg')% 计算运行时间elapsed_time=cputime-start_time,%计算NC（归一化相关系数）NC=nc(message,orig_watermark)% 显示提取水印与原始水印figure(3)subplot(1,2,1);imshow(message,[]);name='提取水印';title(strcat(num2str(name),'NC=',num2str(NC)));subplot(1,2,2)imshow(orig_watermark,[])title('原始水印');
其运行效果如图：
3.2基于DWT变换的数字水印嵌入与提取流程的代码实现
该部分代码实现中用到了随机数种子，上面的基于DFT数字水印方法也用到了随机数。所以在使用grng(seed)需要注意之前是否有使用过随机数，具体注意事项在代码注释有解释。该部分代码分为三部分，主函数、嵌入函数、提取函数。其代码如下：
clc;clear;close all%% 主函数I = imread('xuxian.jpg'); %读取载体图像I = rgb2gray(I); %转换为灰度图W = imread('logo.tif'); %读取水印图像W=W(12:91,17:96); %剪裁为长宽相等figure('Name','载体图像')imshow(I);title('载体图像')figure('Name','水印图像')imshow(W);title('水印图像')ntimes=23; %密钥1， Arnold置乱次数rngseed=59433;%密钥2，随机数种子flag=1;%是否显示中间图像[Iw,psnr]=setdwtwatermark(I,W,ntimes,rngseed,flag); %水印嵌入[Wg,nc]=getdwtwatermark(Iw,W,ntimes,rngseed,flag);%水印提取%% 嵌入函数function [Iw,psnr]=setdwtwatermark(I,W,ntimes,rngseed,flag) %小波水印嵌入type=class(I); %数据类型I=double(I); %强制类型转换为doubleW=logical(W); %强制类型转换为logical[mI,nI]=size(I);[mW,nW]=size(W);if mW~=nW%由于Arnold置乱只能ui方正图像进行处理error('SETDWTWATERMARK:ARNOLD','ARNOLD置乱要求水印图像长宽必须相等！')end%对载体图像进行小波分解%一级Harr小波分解%低频，水平，垂直，对角线[ca1,ch1,cv1,cd1]=dwt2(I,'haar');%二级小波分解[ca2,ch2,cv2,cd2]=dwt2(ca1,'haar');if flagfigure('Name','载体小波分解')subplot(121)imagesc([wcodemat(ca1),wcodemat(ch1);wcodemat(cv1),wcodemat(cd1)])title('一级小波分解')subplot(122)imagesc([wcodemat(ca2),wcodemat(ch2);wcodemat(cv2),wcodemat(cd2)])title('二级小波分解')end%对水印图像进行预处理%初始化置乱数组Wa=W;%对水印进行Arnold变换H=[1,1;1,2]^ntimes;for i=1:nWfor j=1:nWidx=mod(H*[i-1;j-1],nW)+1;Wa(idx(1),idx(2))=W(i,j);endendif flagfigure('Name','水印置乱效果')subplot(121)imshow(W)title('原始水印')subplot(122)imshow(Wa)title(['置乱水印，变换次数=',num2str(ntimes)]);end%小波数字水印的嵌入%初始化嵌入水印的ca2系数ca2w=ca2;%从ca2中随机选择mW*nW个系数rng(rngseed); % 确保您在执行此程序之前没有生成过其他随机数，如果有可使用 rng('default')或者重新启动 MATLABidx=randperm(numel(ca2),numel(Wa));% 将水印信息嵌入到ca2中for i=1:numel(Wa)%二级小波系数c=ca2(idx(i));z=mod(c,nW);%添加水印信息if Wa(i) %水印对应二进制位1if z<nW/4f=c-nW/4-z;elsef=c+nW*3/4-z;endelse %水印对应二进制位0if z<nW*3/4f=c+nW/4-z;elsef=c+nW*5/4-z;endend%嵌入水印后的小波系数ca2w(idx(i))=f;end%根据小波系数重构图像ca1w=idwt2(ca2w,ch2,cv2,cd2,'haar');Iw=idwt2(ca1w,ch1,cv1,cd1,'haar');Iw=Iw(1:mI,1:nI);%计算水印图像峰值信噪比mn=numel(I);Imax=max(I(:));psnr=10*log10(mn*Imax^2/sum((I(:)-Iw(:)).^2));%输出嵌入水印图像最后结果I=cast(I,type);Iw=cast(Iw,type);if flagfigure('Name','嵌入水印的图像')subplot(121)imshow(I);title('原始图像')subplot(122);imshow(Iw);title(['添加水印，PSNR=',num2str(psnr)]);end%% 提取函数function [Wg,nc]=getdwtwatermark(Iw,W,ntimes,rngseed,flag) %小波水印提取[mW,nW]=size(W);if mW~=nWerror('GETDWTWATERMARK:ARNOLD','ARNOLD置乱要求水印图像长宽必须相等！')endIw=double(Iw);W=logical(W);ca1w=dwt2(Iw,'haar');ca2w=dwt2(ca1w,'haar');Wa=W;rng(rngseed);idx=randperm(numel(ca2w),numel(Wa));for i=1:numel(Wa)c=ca2w(idx(i));z=mod(c,nW);if z<nW/2Wa(i)=0;elseWa(i)=1;endendWg=Wa;H=[2 -1;-1,1]^ntimes;for i=1:nWfor j=1:nWidx=mod(H*[i-1;j-1],nW)+1;Wg(idx(1),idx(2))=Wa(i,j);endend%提取和原始水印相关系数的计算nc=sum(Wg(:).*W(:))/sqrt(sum(Wg(:).^2))/sqrt(sum(W(:).^2));if flagfigure('Name','数字水印提取结果')subplot(121)imshow(W)title('原始水印')subplot(122)imshow(Wg)title(['提取水印，NC=',num2str(nc)]);end
其运行结果如下图：
关于数字水印的攻击实验我会在后期内容更新，也有更多数字水印算法。另，本文所述内容和代码如有错误，也请同志们多多指教，当然也会有更优化的方法，请多多交流 !也希望该篇文章能为像我一样同为数字水印初学者提供些许帮助。
参考文献：
[1]Soumitra Roy, Arup Kumar Pal. . An indirect watermark hiding in discrete cosine transform–singular value decomposition domain for copyright protection. Royal Society Open Science. Roy S, Pal AK. : 1-22.
[2]Ibrahim, S., Afrakhteh, M. and Salleh, M., “Adaptive watermarking for printed document authentication,” Proc. ICCIT, Seoul, Korea, .
[3] 杨丹，赵海滨，龙哲. MATLAB图像处理实例详解[M].清华大学出版社.
[4] 冈萨雷斯等. 数字图像处理[M].[美].电子工业出版社.
[5] 夏侯玮明. 基于多变换域的鲁棒数字水印算法[D].南昌大学..6
[6]N. M. Makbol and B.E. Khoo,A hybrid robust image watermarking scheme using integer wavelet transform, singular value decomposition and arnold transform[j]Advances in Visual Informatics. Springer, , pp. 36–47.
[7] 谢琨. 基于半色调技术的数字水印技术研究[D].西安电子科技大学..9
基于离散余弦变换（DCT）傅里叶变换（DFT）小波变换（DWT）的彩色图像数字水印的嵌入 提取简介及MATLAB实现
