2）水印分段操作：
因为原音频信号进行小波变换后，利用量化的方法将其分段，所以讲水印音频也进行相应的分段，最后对应相加。
N=L/LS
这里水印分段不一定取得正数段，所以对其取整，将余数归为最后一段。
3）原始音频信号的一维离散小波变换：
选择合适的小波基进行一维小波三级分解。
DL=DWT（AL）=cA3ÅcD3ÅcD2ÅcD1
其中：
cA3和cD3是三级分解的近似分量和细节分量；cD2和cD1是二级和一级小波分解的细节分量。
由于小波分解的近似分量是信号的低频部分，往往是最重要的，水印嵌入在这部分可以增强水印的稳健性。
因此，提取这部分小波系数来进行下一步的变换。
4）特征区的检测：
因为离散余弦变换的中低频系数集中了信号的大部分特征，同时也是数值较大的部分，所以将水印信号嵌入在此，一是水印的嵌入对其影响不大；二是水印的稳健性也会加强。
因此，通过排序将满足水印长
度的最大的离散余弦系数作为嵌入水印的特征点。
设这些点所组成的序列为
Ck={MAX（CL’），0£k£K}
5）水印信号的嵌入：
这里通过修改系数来进行水印的嵌入,设C*为嵌入水印后的音频信号
* ìCk（1+a）if v（k）=1
Ck=í
（+a）if
v（k）=0
（0≤k≤K）
îCk1
这里的a是大于0的比例因子，通过调节它的大小，在具有听觉不可见性的同时，保证所嵌入的水印信息强度足够大，便于水印信息的正确提取。
在此期间其他的离散余弦系数值不变。
6）离散小波逆变换：
以C*’代替5式的cA3作离散小波变换，得到嵌入水印后的小波变换域的表示式为
D’L=C*’ÅcD3ÅcD2ÅcD1
将其倒置
A’L=D
将A’L代替AL最终得到含有水印的音频信号
Aw=AH+A’L
最后将其输出。
2.3.2水印提取
1）读取嵌入水印的音频信号。
对信号进行三级小波分解，再次提取其低频分量ca31和其他主要分量cd31、cd21、cd11。
2）读取原音频信号，对进行一级小波分解，读取低频分量ca3。
3）对低频分量ca31和ca3进行分段，使用嵌入算法中分段的方式，找到水印嵌入的数据点。
使用嵌入水印音频信号的低频分量减去元信号低频分量，即可得到原水印信号c3。
c3=ca31–ca3
2.4算法的流程图
置乱
数字水印信号
数据嵌入
特征点检测
源数字音频信号
小波D
级分量
其他D级小分量
分段
嵌入水印点的音频
分段
图2-1算法的流程图
3.课程设计的步骤及结果
3.1读入原始信号
clareset;
Fs=44100;
FILE1='Rihanna-TakeABow.wav';
[y,Fs,bits]=wavread（FILE1）;%音频信号y，采样率fs,采样精度bitsy_fft=fft（y,Fs）;
y_fft_f=2*sqrt（y_fft.*conj（y_fft））;axes（handles.axes1）;
plot（y）;gridon;axistight;
title（'原始音频信号的时域波形'）;
xlabel（'time（s）'）;ylabel（'幅度'）;
axes（handles.axes2）;plot（y_fft_f）;gridon;axistight;
title（'原始音频信号的频域波形'）;
xlabel（'f（Hz）'）;ylabel（'幅度'）;
sound（y,Fs）;%播放原始语音
3.2读取原始音频
clareset;
Fs=44100;
FILE1='Rihanna-TakeABow.wav';
[y,Fs,bits]=wavread（FILE1）;%音频信号y，采样率fs,采样精度bits
%用db4小波对原始音频信号进行3级小波分解
[c,l]=wavedec（y,3,'db4'）;%3级小波分解，低频部分为相似，高频部分为细节
%提取3级小波分解的低频系数和高频系数
ca3=appcoef（c,l,'db4',3）;%提取三级小波分解的最低频分
cd3=detcoef（c,l,3）;%提取三级小波分解的次低频分cd2=detcoef（c,l,2）;
cd1=detcoef（c,l,1）;
x=ca3;%提取三级小波分解的最低频部分
%找到插入位置，检测特征点s=max（abs（x））*0.2;i=find（abs（x）>s）;lx=length（x（i））;
%读取水印音频FILE2='test_new.wav';[mark,Fs,bits]=wavread（FILE2）;mark=mark（1:
lx）;mark_fft=fft（mark,Fs）;
mark_fft_f=2*sqrt（mark_fft.*conj（mark_fft））;axes（handles.axes1）;
plot（mark）;
gridon;axistight;
title（'水印音频信号的时域波形'）;
xlabel（'time（s）'）;ylabel（'幅度'）;axes（handles.axes2）
plot（mark_fft_f）;gridon;axistight;
title（'水印音频信号的频域波形'）;
xlabel（'f（Hz）'）;ylabel（'幅度'）;sound（mark,Fs）;
3.3水印嵌入
clareset;
Fs=44100;
FILE1='Rihanna-TakeABow.wav';
[y,Fs,bits]=wavread（FILE1）;%音频信号y，采样率fs,采样精度bits
%用db4小波对原始音频信号进行3级小波分解
[c,l]=wavedec（y,3,'db4'）;%3级小波分解，低频部分为相似，高频部分为细节
%提取3级小波分解的低频系数和高频系数
ca3=appcoef（c,l,'db4',3）;%提取三级小波分解的最低频分
cd3=detcoef（c,l,3）;%提取三级小波分解的次低频分cd2=detcoef（c,l,2）;
cd1=detcoef（c,l,1）;
x=ca3;%提取三级小波分解的最低频部分
%找到插入位置，检测特征点s=max（abs（x））*0.2;i=find（abs（x）>s）;lx=length（x（i））;
%读取水印音频FILE2='test_new.wav';[mark,Fs,bits]=wavread（FILE2）;mark=mark（1:
lx）;
%水印信号嵌入ss=mark（1:
lx）;rr=ss*0.02;x（i）=x（i）.*（1+rr'）;
%小波重构，生成加入了水印信号的音频信号c1=[x',cd3',cd2',cd1'];
s1=waverec（c1,l,'db4'）;
%把加入了水印的原始音频信号作为final1.wav保存FILE3='final1.wav'
y1_fft_f=2*sqrt（y1_fft.*conj（y1_fft））;axes（handles.axes1）;
plot（y1）;gridon;axistight;
title（'嵌入水印后的原始音频信号的时域波形'）;
xlabel（'time（s）'）;ylabel（'幅度'）;axes（handles.axes2）;plot（y1_fft_f）;gridon;axistight;
title（'嵌入水印后的原始音频信号的频域波形'）;
xlabel（'f（Hz）'）;ylabel（'幅度'）;ly1=length（y1）;y1=y1（1:
0.5*ly1）;
sound（y1,Fs）;
3.4水印提取
clareset;
Fs=44100;
FILE1='Rihanna-TakeABow.wav';
[y,Fs,bits]=wavread（FILE1）;%音频信号y，采样率fs,采样精度bits
%用db4小波对原始音频信号进行3级小波分解
[c,l]=wavedec（y,3,'db4'）;%3级小波分解，低频部分为相似，高频部分为细节
%提取3级小波分解的低频系数和高频系数
ca3=appcoef（c,l,'db4',3）;%提取三级小波分解的最低频分
cd3=detcoef（c,l,3）;%提取三级小波分解的次低频分cd2=detcoef（c,l,2）;
cd1=detcoef（c,l,1）;
x=ca3;%提取三级小波分解的最低频部分
%找到插入位置，检测特征点s=max（abs（x））*0.2;i=find（abs（x）>s）;lx=length（x（i））;
%读取含有水印的原始音频FILE3='final1.wav';[y1,Fs1,bits1]=wavread（FILE3）;
%用db4小波对含有水印的原始音频信号进行3级小波分解[c1,l1]=wavedec（y1,3,'db4'）;
%提取3级小波分解的低频系数和高频系数ca31=appcoef（c1,l1,'db4',3）;cd31=detcoef（c1,l1,3）;cd21=detcoef（c1,l1,2）;cd11=detcoef（c1,l1,1）;
x1=ca31;
%水印信号提取z（i）=x1（i）-x（i）;
s2=z（i）./x（i）';s2=s2/0.02;
%把提取的水印信号作为final2.wav保存FILE4='final2.wav';wavwrite（s2,Fs,bits,FILE4）;[mark1,Fs,bits]=wavread（FILE4）;mark1_fft=fft（mark1）;
mark1_fft_f=2*sqrt（mark1_fft.*conj（mark1_fft））;axes（handles.axes1）
plot（mark1）;gridon;axistight;
title（'提取出的水印音频信号的时域波形'）;
xlabel（'time（s）'）;ylabel（'幅度'）;axes（handles.axes2）;plot（mark1_fft_f）;gridon;axistight;
title（'提取出的水印音频信号的频域波形'）;
xlabel（'f（Hz）'）;ylabel（'幅度'）;sound（mark1,Fs）;
3.5低通滤波
%读取含有水印的原始音频FILE3='final1.wav';[y1,Fs,bits]=wavread（FILE3）;[b,a]=butter（1,0.1）;s3=filtfilt（b,a,y1）;
%把经过滤波后的含有水印的原始音频信号作为final3.wav保存FILE5='final3.wav';
wavwrite（s3,Fs,bits,FILE5）;
[y_filter,Fs,bits1]=wavread（FILE5）;y_filter_fft=fft（y_filter）;axes（handles.axes1）;plot（y_filter）;gridon;axistight;
title（'经过滤波后的含有水印的原始音频信号的时域波形'）;
xlabel（'time（s）'）;ylabel（'幅度'）;axes（handles.axes2）;plot（y_filter_fft）;gridon;axistight;
title（'经过滤波后的含有水印的原始音频信号的频域波形'）;
xlabel（'f（Hz）'）;ylabel（'幅度'）;ly2=length（y_filter）;y_filter=y_filter（1:
0.5*ly2）;sound（y_filter,Fs）;
3.6低通滤波后提取水印
clareset;
Fs=44100;
FILE1='Rihanna-TakeABow.wav';
[y,Fs,bits]=wavread（FILE1）;%音频信号y，采样率fs,采样精度bits
%用db4小波对原始音频信号进行3级小波分解
[c,l]=wavedec（y,3,'db4'）;%3级小波分解，低频部分为相似，高频部分为细节
%提取3级小波分解的低频系数和高频系数
ca3=appcoef（c,l,'db4',3）;%提取三级小波分解的最低频分
cd3=detcoef（c,l,3）;%提取三级小波分解的次低频分cd2=detcoef（c,l,2）;
cd1=detcoef（c,l,1）;
x=ca3;%提取三级小波分解的最低频部分
%找到插入位置，检测特征点s=max（abs（x））*0.2;i=find（abs（x）>s）;lx=length（x（i））;
%读取经过滤波后的含有水印的原始音频FILE5='final3.wav';[y_filter,Fs,bits1]=wavread（FILE5）;
%用db4小波对经过滤波后的含有水印的原始音频信号进行3级小波分解
[c2,l1]=wavedec（y_filter,3,'db4'）;
%提取3级小波分解的低频系数和高频系数ca32=appcoef（c2,l1,'db4',3）;cd32=detcoef（c2,l1,3）;
cd22=detcoef（c2,l1,2）;cd12=detcoef（c2,l1,1）;x2=ca32;
%水印信号提取z（i）=x2（i）-x（i）;
s3=z（i）./x（i）';s3=s3/0.02;
%把提取的水印信号作为final4.wav保存FILE6='final4.wav';wavwrite（s3,Fs,bits,FILE6）;[mark_filter,Fs,bits1]=wavread（FILE6）;mark_filter_fft=fft（mark_filter）;axes（handles.axes1）;
plot（mark_filter）;gridon;axistight;
title（'经过滤波器处理后提取的水印音频信号的时域波形'）;
xlabel（'time（s）'）;ylabel（'幅度'）;axes（handles.axes2）;plot（mark_filter_fft）;gridon;axistight;
title（'经过滤波器处理后提取的水印音频信号的频域波形'）;
xlabel（'f（Hz）'）;ylabel（'幅度'）;sound（mark_filter,Fs）;
3.7界面的设计
1）在Commandwindow窗口下输入guide命令，并回车；
2）在CreatNewGUI界面下选择BlankGUI（Default），点击OK；
3）在新的界面下，点击Axe按钮，生成一个区域，再点击PushB
升级会员 