https://blog.csdn.net/weixin_39707121/article/details/79407070
https://www.jiamisoft.com/blog/2914-shuzishuiyinjiamisuanfa.html
%
% 音频水印技术 教学版演示程序
%
% 车生兵 著
%
% 采用一级小波变换低频系数量化嵌入方式
%
% 变量说明：
%
% delt： 量化步长； w: 水印图像； rem: 量化后余数；
% p: 攻击概率； w1: 提取水印； LL1： 调整后的低频量化系数；
% I0： 原始音频； I1： 载体音频； I11： 嵌入水印后的局部块；
% [m,n]: 水印尺寸； n0: 量化系数； p: 攻击概率；
% [LL LH HL HH]: 小波变换后得到的低频、中频、高频系数；
%
% Step 1: 数据准备
delt=2; %选取大规模实验获得的最佳值
p=0.5; %选取大规模实验获得的最佳值
[Y,fs,nbits]=wavread('oriaud.wav');
Y1=Y;
tmp1=(Y(1:256*256,1)+1)*128;
I0=reshape(tmp1,256,256);
I1=I0; %载体音频初始值
w=imread('w.bmp'); %获取水印图像
[m,n]=size(w); %获取水印图像的尺寸
% Step 2: 水印嵌入
for i=1:1:m %按照从上到下，从左至右顺序嵌入水印信息
for j=1:1:n
I=I0(2*i-1:2*i,2*j-1:2*j); %获取嵌入w(i,j)水印值的原始图像块I
[LL,LH,HL,HH]=dwt2(I,'haar'); %获得haar小波变化系数
n0=floor(LL/delt); %获得小波变换低频系数的量化系数
rem=LL-n0*delt; %获得小波变换低频系数量化后的余数
if (mod(n0,2)==w(i,j)) %判断：量化系数与待嵌入水印系数的奇偶性是否吻合？
LL1=n0*delt+p*delt; %吻合：量化系数不变，调整低频系数值
else %不吻合：量化系数修改
if(rem>p*delt) %判断：余数大于攻击最佳调整值p*delt？
LL1=(n0+1)*delt+p*delt; %是：量化系数增加1，调整低频系数值
else
LL1=(n0-1)*delt+p*delt; %否：量化系数减少1，调整低频系数值
end
end
I11=idwt2(LL1,LH,HL,HH,'haar'); %逆小波变换，获得初始载体数据
I1(2*i-1:2*i,2*j-1:2*j)=I11; %获得整形后的载体数据，这里是用round()进行整形
end
end
% Step 3：提取水印
[M,N]=size(I1); %获得载体音频尺寸
m=M/2; %计算提取水印图像行数
n=N/2; %计算提取水印图像列数
w1=w; %提取水印图像赋初始值
for i=1:1:m %按照从上到下，从左至右顺序提取水印信息
for j=1:1:n
I=I1(2*i-1:2*i,2*j-1:2*j); %获取载体数据块
[LL,LH,HL,HH]=dwt2(I,'haar'); %做haar小波变换
n0=floor(LL/delt); %提取量化系数
if (mod(n0,2)==0) %判断：量化系数奇偶性是否为0？
w1(i,j)=0; %是：提取出的水印值为0
else
w1(i,j)=1; %否：提取出的水印值为1
end
end
end
% Step 4：显示结果
tmp1=(I1./128-1);
Y1(1:256*256,1)=reshape(tmp1,256*256,1);
imwrite(w1,'w1.bmp');
wavwrite(Y1,fs,nbits,'embeded.wav');
figure,subplot(2,2,1);plot(Y);title('原始音频');
subplot(2,2,2);plot(Y1);title('载体音频');
subplot(2,2,3);imshow(w);title('水印图像');
subplot(2,2,4);imshow(w1);title('提取水印');
% 程序结束
%
% 视频水印 演示程序
%
% 车生兵 著
%
% 截取视频中的图像
clc;
clear;
close all;
% 将电影转成图片序列
mov=aviread('avi\search.avi');
fi = AVIINFO('avi\search.avi');
% [m,n]=size(mov(1).cdata);
fnum=size(mov,2); % 读取电影的帧数，“2”表示返回的为列数
for i=1:fnum
strtemp_rgb=strcat('oringalframe\',int2str(i),'.bmp'); % 将每帧转成bmp的图片
imwrite(mov(i).cdata,mov(i).colormap,strtemp_rgb); % 写入真彩图像
end
%
% 把Index类型图像转换成RGB图像
%
for i=1:fnum
strtemp_rgb=strcat('oringalframe\',int2str(i),'.bmp'); % 将每帧转成bmp的图片
strtemp_rgb1=strcat('rgbframe\',int2str(i),'.bmp'); % 将每帧转成bmp的图片
[X,MAP]=imread(strtemp_rgb);
RGB = ind2rgb(X,MAP);
imwrite(RGB,strtemp_rgb1); % 写入真彩图像
end
% 创建文件夹
% 文件夹"FlowerGarden_rgb_embed"存放写入了水印的视频帧图像
mkdir('rgb_embed');
copyfile('rgbframe','rgb_embed');
% 设置用来选择待嵌入视频帧keyframe的随机数key,注意，关键帧的选择要随机，而且1<=key<=fi.NumFrames
key=5;
strtemp_rgb2=strcat('rgb_embed\',int2str(key),'.bmp');
keyframe=imread(strtemp_rgb2);
I0=rgb2ycbcr(keyframe);
[m,n]=size(I0(:,:,1));
%
I11=I0(:,:,1);
%
% Step 1: 数据准备
delt=12; %选取大规模实验获得的最佳值
p=0.5; %选取大规模实验获得的最佳值
I1=I11; %载体图像初始值
w=imread('w.bmp'); %获取水印图像
[m,n]=size(w); %获取水印图像的尺寸
% Step 2: 水印嵌入
for i=1:1:m %按照从上到下，从左至右顺序嵌入水印信息
for j=1:1:n
I=I0(2*i-1:2*i,2*j-1:2*j); %获取嵌入w(i,j)水印值的原始图像块I
[LL,LH,HL,HH]=dwt2(I,'haar'); %获得haar小波变化系数
n0=floor(LL/delt); %获得小波变换低频系数的量化系数
rem=LL-n0*delt; %获得小波变换低频系数量化后的余数
if (mod(n0,2)==w(i,j)) %判断：量化系数与待嵌入水印系数的奇偶性是否吻合？
LL1=n0*delt+p*delt; %吻合：量化系数不变，调整低频系数值
else %不吻合：量化系数修改
if(rem>p*delt) %判断：余数大于攻击最佳调整值p*delt？
LL1=(n0+1)*delt+p*delt; %是：量化系数增加1，调整低频系数值
else
LL1=(n0-1)*delt+p*delt; %否：量化系数减少1，调整低频系数值
end
end
I11=idwt2(LL1,LH,HL,HH,'haar'); %逆小波变换，获得初始载体数据
I1(2*i-1:2*i,2*j-1:2*j)=round(I11); %获得整形后的载体数据，这里是用round()进行整形
end
end
% Step 3：提取水印
% 应该从载体视频帧中获取嵌入帧数据，在这里直接从嵌入程序中直接获取，由于写入质量fi.Quality为100，所以从载体视频中读取出来也是一样的
[M,N]=size(I1); %获得载体图像尺寸
m=M/2; %计算提取水印图像行数
n=N/2; %计算提取水印图像列数
w1=w; %提取水印图像赋初始值
for i=1:1:m %按照从上到下，从左至右顺序提取水印信息
for j=1:1:n
I=I1(2*i-1:2*i,2*j-1:2*j); %获取载体数据块
[LL,LH,HL,HH]=dwt2(I,'haar'); %做haar小波变换
n0=floor(LL/delt); %提取量化系数
if (mod(n0,2)==0) %判断：量化系数奇偶性是否为0？
w1(i,j)=0; %是：提取出的水印值为0
else
w1(i,j)=1; %否：提取出的水印值为1
end
end
end
% Step 4：显示结果
imwrite(w1,'w1.bmp');
figure,subplot(1,2,1);imshow(w);title('水印图像');
subplot(1,2,2);imshow(w1);title('提取水印');
%
%I12(:,:,1)=I11;
I12(:,:,1)=I1;
I12(:,:,2)=I0(:,:,2);
I12(:,:,3)=I0(:,:,3);
I1=ycbcr2rgb(I12);
% 写入嵌入水印的帧
imwrite(I1,strtemp_rgb2);
% 将嵌入水印的视频帧序列转换成avi视频
% 定义一个avi文件
aviobj=avifile('avi\rgb_embed.avi');
aviobj.Quality =fi.Quality;
aviobj.compression='None';
aviobj.colormap=gray(256);
aviobj.Fps=fi.FramesPerSecond;
% 使用addframe把图片写入视频
for i=1:fnum
picdata=imread(strcat('rgb_embed','\',num2str(i),'.bmp')); % 获得一帧图像
aviobj=addframe(aviobj,uint8(picdata)); % 加到电影剪辑文件中
end
aviobj=close(aviobj); % 关闭文件，结束数值仿真模拟过程
%
% 数字图像水印技术 教学版演示程序
%
% 车生兵 著
%
% 采用一级小波变换低频系数量化嵌入方式
%
% 变量说明：
%
% delt： 量化步长； w: 水印图像； rem: 量化后余数；
% p: 攻击概率； w1: 提取水印； LL1： 调整后的低频量化系数；
% I0： 原始图像； I1： 载体图像； I11： 嵌入水印后的局部块；
% [m,n]: 水印尺寸； n0: 量化系数； p: 攻击概率；
% [LL LH HL HH]: 小波变换后得到的低频、中频、高频系数；
%
% Step 1: 数据准备
delt=12; %选取大规模实验获得的最佳值
p=0.5; %选取大规模实验获得的最佳值
II=imread('lena.jpg'); %获取原始图像
I0=rgb2gray(II); %只对灰度图像做水印嵌入与提取
I1=I0; %载体图像初始值
w=imread('w.bmp'); %获取水印图像
[m,n]=size(w); %获取水印图像的尺寸
% Step 2: 水印嵌入
for i=1:1:m %按照从上到下，从左至右顺序嵌入水印信息
for j=1:1:n
I=I0(2*i-1:2*i,2*j-1:2*j); %获取嵌入w(i,j)水印值的原始图像块I
[LL,LH,HL,HH]=dwt2(I,'haar'); %获得haar小波变化系数
n0=floor(LL/delt); %获得小波变换低频系数的量化系数
rem=LL-n0*delt; %获得小波变换低频系数量化后的余数
if (mod(n0,2)==w(i,j)) %判断：量化系数与待嵌入水印系数的奇偶性是否吻合？
LL1=n0*delt+p*delt; %吻合：量化系数不变，调整低频系数值
else %不吻合：量化系数修改
if(rem>p*delt) %判断：余数大于攻击最佳调整值p*delt？
LL1=(n0+1)*delt+p*delt; %是：量化系数增加1，调整低频系数值
else
LL1=(n0-1)*delt+p*delt; %否：量化系数减少1，调整低频系数值
end
end
I11=idwt2(LL1,LH,HL,HH,'haar'); %逆小波变换，获得初始载体数据
I1(2*i-1:2*i,2*j-1:2*j)=round(I11); %获得整形后的载体数据，这里是用round()进行整形
end
end
% Step 3：提取水印
[M,N]=size(I1); %获得载体图像尺寸
m=M/2; %计算提取水印图像行数
n=N/2; %计算提取水印图像列数
w1=w; %提取水印图像赋初始值
for i=1:1:m %按照从上到下，从左至右顺序提取水印信息
for j=1:1:n
I=I1(2*i-1:2*i,2*j-1:2*j); %获取载体数据块
[LL,LH,HL,HH]=dwt2(I,'haar'); %做haar小波变换
n0=floor(LL/delt); %提取量化系数
if (mod(n0,2)==0) %判断：量化系数奇偶性是否为0？
w1(i,j)=0; %是：提取出的水印值为0
else
w1(i,j)=1; %否：提取出的水印值为1
end
end
end
% Step 4：显示结果
imwrite(w1,'w1.bmp');
imwrite(I1,'I1.bmp');
figure,subplot(2,2,1);imshow(I0);title('原始图像');
subplot(2,2,2);imshow(I1);title('载体图像');
subplot(2,2,3);imshow(w);title('水印图像');
subplot(2,2,4);imshow(w1);title('提取水印');
% 程序结束
% 水印嵌入过程
I=rgb2gray(imread('lena.jpg'));
W=imread('water.bmp');
[m,n]=size(I);
[h,l]=size(W);
delta=12.0;
n=0.0;
rem=0.0;
p=0.5;
I1=I;
for (i=1:1:h)
for (j=1:1:l)
I0=double(I(((i-1)*2+1):(i*2),((j-1)*2+1):(j*2)));
[LL LH HL HH]=dwt2(I0,'haar');
n=floor(LL/delta);
rem=mod(LL,delta);
if (W(i,j)==mod(n,2))
LL1=n*delta+p*delta;
else
if (rem>p*delta)
LL1=(n+1)*delta+p*delta;
if (LL1>(510-delta))
LL1=(n-1)*delta+p*delta;
end
else
LL1=(n-1)*delta+p*delta;
if (LL1<delta)
LL1=(n+1)*delta+p*delta;
end
end
end
I01=idwt2(LL1,LH,HL,HH,'haar');
I11=floor(I01);
I12=mod(I01,1);
xsh=round(sum(sum(I12)));
xx=reshape(I11,1,4);
[yy,ii]=sort(xx);
max=4;
for (r=xsh:-1:1)
for (s=1:1:4)
if (ii(1,s)==max)
break;
end
end
xx(1,s)=xx(1,s)+1;
max=max-1;
end
xxx=reshape(xx,2,2);
I1(((i-1)*2+1):(i*2),((j-1)*2+1):(j*2))=xxx;
end
end
I1=uint8(I1);
imwrite(I,'lena.bmp');
imwrite(I1,'embded.bmp');
figure,subplot(2,2,1),imshow(I);title('Original');
subplot(2,2,2),imshow(I1);title('Embeded');
subplot(2,2,3),imshow(W);title('Water');
% 水印提取过程
x=imread('lena.jpg');
I=imread('embded.bmp');
W=imread('water.bmp');
[m,n]=size(I);
[h,l]=size(W);
delta=12.0;
n=0.0;
p=0.5;
W1=W;
for (i=1:1:h)
for (j=1:1:l)
I0=double(I(((i-1)*2+1):(i*2),((j-1)*2+1):(j*2)));
[LL LH HL HH]=dwt2(I0,'haar');
n=floor(LL/delta);
W1(i,j)=mod(n,2);
end
end
figure,subplot(2,2,1),imshow(x);title('Original');
subplot(2,2,2),imshow(I);title('Embeded');
subplot(2,2,3),imshow(W);title('Water');
subplot(2,2,4),imshow(W1);title('Extract');
