%由高斯正态分布序列g1产生36×4的水印信%号w0,w0由(0,1)组成。clearrandn('state',1106);g1=randn(36,4);fori=1:36forj=1:4ifg1(i,j)>=0w0(i,j)=1;elsew0(i,j)=0;end;end;end;figure;imshow(w0);title('水印');%对水印信号w0进行(7,4)汉明编码,得到一%36×7的分组码x0。x0=w0;fori=1:36s=8*x0(i,1)+4*x0(i,2)+2*x0(i,3)+x0(i,4);switchscase0x0(i,5)=0;x0(i,6)=0;x0(i,7)=0;case1x0(i,5)=0;x0(i,6)=1;x0(i,7)=1;case2x0(i,5)=1;x0(i,6)=1;x0(i,7)=0;case3x0(i,5)=1;x0(i,6)=0;x0(i,7)=1;case4x0(i,5)=1;x0(i,6)=1;x0(i,7)=1;case5x0(i,5)=1;x0(i,6)=0;x0(i,7)=0;case6x0(i,5)=0;x0(i,6)=0;x0(i,7)=1;case7x0(i,5)=0;x0(i,6)=1;x0(i,7)=0;case8x0(i,5)=1;x0(i,6)=0;x0(i,7)=1;case9x0(i,5)=1;x0(i,6)=1;x0(i,7)=0;case10x0(i,5)=0;x0(i,6)=1;x0(i,7)=1;case11x0(i,5)=0;x0(i,6)=0;x0(i,7)=0;case12x0(i,5)=0;x0(i,6)=1;x0(i,7)=0;case13x0(i,5)=0;x0(i,6)=0;x0(i,7)=1;case14x0(i,5)=1;x0(i,6)=0;x0(i,7)=0;case15x0(i,5)=1;x0(i,6)=1;x0(i,7)=1;end;end;%对x0进行行向位扩展,得到一个由(-1,1)组成%的扩展序列y。cr为扩展因子。cr=256;fori=1:252ifx0(i)==1y(i,1:cr)=1;elsey(i,1:cr)=-1;end;end;y(253:256,:)=0;%以下产生伪随机序列p。为此先设定密钥(1114)%并产生高斯正态分布序列g2,再由g2产生由(-1,1)%组成的伪随机序列p。randn('state',1114);g2=randn(256,256);fori=1:256forj=1:256ifg2(i,j)>0p(i,j)=1;elsep(i,j)=-1;end;end;end;yp=y.*p;%下面设定的每类块基准噪声阈值jnd1是通过%反复实验确定出来的。t1=1.7;t2=2.1;t3=2.5;t4=2.9;t5=3.3;t6=3.7;t7=4.1;t8=4.5;%读入原图象并转换成双精度,k是8×8图像%块数。f0=imread('cameraman.tif');f0=double(f0);[c,s]=size(f0);k=c*s/64;%计算每块的方差std2和能量e,并对能量e按%5升序排列。cf0=im2col(f0,[8,8],'distinct');std1=std(cf0);std2=std1.^2;fori=1:ke(i)=sum(cf0(:,i).^2);end;[e1,ind1]=sort(e);%按能量e的索引ind1顺序重排方差std2,将结%果存于std3。j=1;fori=1:kz=ind1(i);std3(j)=std2(z);j=j+1;end;%设定分类界限。m1=median(std3(1:k/4));m2=median(std3(k/4+1:k/2));m3=median(std3(k/2+1:3*k/4));m4=median(std3(3*k/4+1:k));n1=e1(k/4);n2=e1(2*k/4);n3=e1(3*k/4);n4=e1(k);%按能量e和方差std2将原图像块分成8类,并%给每类块赋基准噪声阈值。fori=1:kife(i) < =n1 ifstd2(i) > =m1jnd1(i)=t1;elsejnd1(i)=t2;end;elseif(n1 < e (i)&e(i)< =n2) ifstd2(i) > =m2jnd1(i)=t3;elsejnd1(i)=t4;end;elseif(n2 < e (i)&e(i)< =n3) ifstd2(i) > =m3jnd1(i)=t5;elsejnd1(i)=t6;end;elseif(n3 < e (i)&e(i)< =n4) ifstd2(i) > =m4jnd1(i)=t7;elsejnd1(i)=t8;end;end;end;%计算每类块的噪声阈值jnd,它等于基准噪声%阈值jnd1和附加噪声阈值jnd2之和。deta=0.0035;jnd2=deta*cf0;fori=1:kjnd(:,i)=jnd1(i)+jnd2(:,i);end;%嵌入已调制的水印信号yp。并重构成含水印%的图像f1。recf0=reshape(cf0,256,256);rejnd=reshape(jnd,256,256);recf1=recf0+rejnd.*yp;%重构嵌入水印的图像f1。cf1=reshape(recf1,64,1024);f1=col2im(cf1,[8,8],[256,256],'distinct');%从含水印的图像f1中提取出36×7的分组码%x1。fori=1:252sk(i)=sum((recf1(i,:)-recf0(i,:)).*p(i,:));end;fori=1:252ifsign(sk(i))==-1rex1(i)=0;elserex1(i)=1;end;end;x1=reshape(rex1,36,7);%对提取出来的分组码x1进行纠错解码,最后%得到一个36×4的水印w1。fori=1:36s1(i)=x1(i,1)+x1(i,2)+x1(i,3)+x1(i,5);s2(i)=x1(i,2)+x1(i,3)+x1(i,4)+x1(i,6);s3(i)=x1(i,1)+x1(i,2)+x1(i,4)+x1(i,7);s1(i)=mod(s1(i),2);s2(i)=mod(s2(i),2);s3(i)=mod(s3(i),2);if(s1(i)==0 &s2 (i)==0)&(s3(i)==1)x1(i,7)=~x1(i,7);elseif(s1(i)==0 &s2 (i)==1)&(s3(i)==0)x1(i,6)=~x1(i,6);elseif(s1(i)==0 &s2 (i)==1)&(s3(i)==1)x1(i,4)=~x1(i,4);elseif(s1(i)==1 &s2 (i)==0)&(s3(i)==0)x1(i,5)=~x1(i,5);elseif(s1(i)==1 &s2 (i)==0)&(s3(i)==1)x1(i,1)=~x1(i,1);elseif(s1(i)==0 &s2 (i)==1)&(s3(i)==0)x1(i,3)=~x1(i,3);elseif(s1(i)==1 &s2 (i)==1)&(s3(i)==1)x1(i,2)=~x1(i,2);end;end;w1=x1(:,1:4);%计算原图像f0和嵌入水印的图像f1的信噪比%snr。v1=sum(sum(f0.^2));v2=sum(sum((f0-f1).^2));snr=10*log10(v1/v2);%计算原水印w0和提取出来的水印w1的相关%系数r。r=corr2(w0,w1);%显示原始图像f0和嵌入水印的图像f1,结果%示于图1、图2。figureimshow(f0,[]);title('原图像');figureimshow(f1,[]);title('嵌入水印的图像');%将含水印图像f1归一化,以便于攻击处理。m=max(max(f1));f=f1/m;%1.JPEG压缩。imwrite(f,'attackf.jpg','jpg','quality',30);attackf=imread('attackf.jpg');attackf=double(attackf)/255;%2.高斯低通滤波。h=fspecial('gaussian',3,1);attackf=filter2(h,f);%3.直方图均衡化。attackf=histeq(f);%4.图像增亮。attackf=imadjust(f,[],[0.4,1]);%5.图像变暗。attackf=imadjust(f,[],[0,0.85]);%6.增加对比度。attackf=imadjust(f,[0.3,0.6],[]);%7.降低对比度。attackf=imadjust(f,[],[0.2,0.8]);%8.添加高斯噪声。attackf=imnoise(f,'gaussian',0,0.01);%9.增加黑白像素点。attackf=imnoise(f,'salt&pepper',0.06);%10添加乘积性噪声。attackf=imnoise(f,'speckle',0.08);%显示攻击后的水印图像attackf,为了简单起见%这里只给出最后一种攻击即添加乘性噪声后的%水印图像,其它攻击情况相同。figureimshow(attackf,[]);title('攻击后的水印图像');