watermark1.m

用MATLAB实现数字音频水印的嵌入和提取

% 加入白噪声的音频水印程序
clear;
kk.wave = wavread('3.wav');   %读入原始音频文件
y=kk.wave;
[c,l]=wavedec(y,3,'db4');               %三级小波分解
ca3=appcoef(c,l,'db4',3);
cd3=detcoef(c,l,3);
cd2=detcoef(c,l,2);
cd1=detcoef(c,l,1);
x=ca3;                              %提取低频系数
len=length(y);
x1=x;
s=max(abs(x))*0.2;
i=find(abs(x)>s);lx=length(x(i));         %找出大于最大值0.2倍的序列
figure;
subplot(2,2,1);
plot(ca3);                            %画出低频系数图
title('低频系数图形');
subplot(2,2,2);
plot(cd3);
title('cd3');
subplot(2,2,3);
plot(cd2);
axis([0 10e4 -0.5 0.5]);
title('cd2');
subplot(2,2,4);
plot(cd1);
title('cd1');
randn('seed',10);                      %产生随机高斯序列
mark=randn(1,lx);
ss=mark;
rr=ss*0.1;                           %设置水印嵌入强度
x(i)=x(i).*(1+2*rr');                   %嵌入水印
c1=[x',cd3',cd2',cd1'];
s1=waverec(c1,l,'db4');
file1='已加水印.wav';
dd=length(s1);
ee=reshape(s1,dd/2,2);
wavwrite(ee,file1);
figure;
subplot(3,1,1);plot(y);                   %画出原信号图
axis([0 18e4 -2 2]);
title('原信号的图');                    
subplot(3,1,2);plot(ss);
title('水印图');
subplot(3,1,3);plot(s1);                  %画出嵌入了水印的信号图
title('加入了水印的声音信号')
kk.wave = wavread('已加水印'); 
yc=kk.wave;
fz=sum(y.*y);                           %计算嵌入了水印的信号的信噪比
fm=sum((y-yc).*(y-yc));
SNR=-10*log(fm/fz)
yyy=randn(1,dd);                            %加入白噪声                          
b=sqrt(0.01); 
yyy=b*yyy; 
s1=s1+yyy;
ee=reshape(s1,dd/2,2);
wavwrite(ee,file1);
kk.wave = wavread('已加水印');            %读入声音文件
yr=kk.wave;
[cr,lr]=wavedec(yr,3,'db4');
car3=appcoef(cr,lr,'db4',3);
cdr3=detcoef(cr,lr,3);
cdr2=detcoef(cr,lr,2);
cdr1=detcoef(cr,lr,1);
xr=car3;
figure;
rrr=((xr(i)./x1(i))-1)/2;
subplot(2,1,1);plot(rr);
title('原水印的图');
subplot(2,1,2);plot(rrr);
axis([0 16e3 -0.4 0.4]);
title('加入白噪声后提取水印的图');
whos('rrr');
figure;
rr1=reshape(rrr,1,lx);
syc=rr1-rr;
plot(syc);                                  %画出水印差别图
axis([0 16e3 -0.3 0.3]);
title('水印之差');
rrr=((xr(i)./x1(i))-1)/2;
d=length(rr);
x=0;
y=0;
z=0;
for i=1:d
        x=x+rr(i)*rrr(i);
        y=y+rr(i)^2;
        z=z+rrr(i)^2;
    end
p=x/((y^0.5)*(z^0.5))                         %计算相关度