【语音隐写】基于matlab小波变换算法求解水印嵌入提取【含Matlab源码 513期】

一、简介

随着计算机和网络的飞速发展，人们的许多创作和成果都以数字形式进行存储和发布。然而，数字作品极易被非法拷贝、伪造和窜改，使得很多版权所有者不愿意利用网络公开其作品，从而阻碍其自身发展。目前，数字作品的版权保护不仅仅是立法问题，也是一个很重要的技术难题。从技术上看，数字媒体版权信息的嵌入和检测问题，是数字作品版权保护的两个关键问题，它综合了传统密码学的认证和鉴别问题的特点，又加入了稳健性要求。版权保护信息必须与被保护的数据密切结合，版权保护信息的鉴别过程必须具有了抗干扰能力，在这种情况下，数字水印技术应运而生。

1. 数字水印技术概述

提到水印，人们都会想到钞票中的水印。钞票水印具有两条特性，首先，水印在通常情况下不可见，只有在特殊的观察条件下才可见（钞票中水印在光下可见）。其次，水印信息必须与载体对象相关（在这里表示纸币的真实性）。

1.1 数字水印技术的特性

可证明性：能够为受到版权保护的数字产品提供完全可靠的证据。

不可见性：即被嵌入水印信息的数字产品不会出现明显的质量下降，隐藏的数据不易被察觉；另外，不能用统计的方法恢复出水印。

鲁棒性：添加的数字水印必须对施加于宿主图象的攻击具有一定的免疫能力，不能因为对宿主图象的某种操作而导致水印信息丢失。

1.2 数字水印技术的应用

水印技术的应用极为广泛，主要有以下7种应用领域：广播监控、所有者识别、所有权验证、交易跟踪、内容真伪鉴别、拷贝控制以及设备控制。

2. 理论基础

2.1 小波变换理论

自1986年以来，小波分析的理论、方法与应用的研究一直方兴未艾。作为一种数学工具，小波变换是对人们熟知的傅立叶变换和窗口傅立叶变换的一个重大突破，为信号分析、图像处理及其它非线性科学研究领域带来了革命性的影响。

人类视觉系统（HVS）的文理特性和照亮掩蔽特性表明，纹理越复杂，背景的亮度越亮，人类视觉对其轻微变化就越不敏感。大量的研究表明，人眼在处理图像信号时，将图像滤波成若干的子带信号，他们占据不同的频率范围，即图像在HVS中被认为是由不同频率范围的信息组成。其特征为：人眼对反映局部结构的边缘和轮廓不敏感；对低频信号，表现出较高的灵敏度。HVS在同一品大范围对不同方向纹理细节信号等表现出不同的灵敏度，这一特点与小波变换的多分辨率分析具有一定的相似性。小波变换是傅立叶变换的发展，是空间和频率的局部变换，它在频域和时域同时具有良好的局部化特征。小波变换在图像处理中的基本思想是把图像进行多分辨率分解成不同的空间和独立的频率带的子图像，然后对子图像的系数进行处理。

根据S. Mallat的塔式分解算法，图像经过小波变换后分解成四个子图：水平方向LH、垂直方向HL和对角线方向HH的中高频细节子图和低频逼近子图LL。低频部分还可以继续分解，产生三个高频带系列LHn、HLn、HHn（n=1，2，3）和一个低频带LL3（见图1）。

图1中的LL3表示小波变换分解级数决定的最大尺度、最小分辨率下对原始图像的最佳逼近，它的同级特征和原理图相似，图像大部分能量集中于此。高频带系列代表图像的边缘和纹理。

2.2 小波变换的应用

小波分析的应用领域十分广泛，它包括：数学领域的许多学科；信号分析、图象处理；量子力学、理论物理；军事电子对抗与武器的智能化；计算机分类与识别；音乐与语言的人工合成；医学成像与诊断；地震勘探数据处理；大型机械的故障诊断等方面；例如，在数学方面，它已用于数值分析、构造快速数值方法、曲线曲面构造、微分方程求解、控制论等。在信号分析方面的滤波、去噪声、压缩、传递等。在图象处理方面的图象压缩、分类、识别与诊断，去污等。在医学成像方面的减少B超、CT、核磁共振成像的时间，提高分辨率等。
⑴小波分析用于信号与图象压缩是小波分析应用的一个重要方面。它的特点是压缩比高，压缩速度快，压缩后能保持信号与图象的特征不变，且在传递中可以抗干扰。基于小波分析的压缩方法很多，比较成功的有小波包最好基方法，小波域纹理模型方法，小波变换零树压缩，小波变换向量压缩等。
⑵小波在信号分析中的应用也十分广泛。它可以用于边界的处理与滤波、时频分析、信噪分离与提取弱信号、求分形指数、信号的识别与诊断以及多尺度边缘检测等。
⑶在工程技术等方面的应用。包括计算机视觉、计算机图形学、曲线设计、湍流、远程宇宙的研究与生物医学方面。
从图像处理的角度看，小波变换存在以下几个优点：
⑴小波分解可以覆盖整个频域（提供了一个数学上完备的描述）
⑵小波变换通过选取合适的滤波器，可以极大的减小或去除所提取得不同特征之间的相关性
⑶小波变换具有“变焦”特性，在低频段可用高频率分辨率和低时间分辨率（宽分析窗口），在高频段，可用低频率分辨率和高时间分辨率（窄分析窗口）
⑷小波变换实现上有快速算法（Mallat小波分解算法）

3.DWT变换域数字水印技术

3.1数字水印嵌入技术

图像的水印技术根据水印嵌入的方式可以分为两类：时/空域技术（水印被直接嵌入在图像的亮度值中）和变换域技术（将图像做某种数学变换，然后水印被嵌入于变换系数中）。早期人们对水印的研究基本上是基于时空域的，算法相对简单，实时性较强，但在鲁棒性上不如变换域算法。目前变换域方法正日益普遍，有DCT、DWT、DFT变换域算法。

变换域算法的优点是：1、水印信息分布到空间域的所有像素上，有利于提高水印的不可见性。2、能方便的与HVS(人类视觉系统)的某些特性结合。3、很好的鲁棒性，对图像压缩、常用的图像滤波以及噪声均有一定的抵抗力。离散余弦变换是从图像空间到频率空间的全局变换，由于离散余弦变换的全局本质，在变换空间中任何一个数据的误差都会影响到图像中的每一个像素。利用小波变换把原始图像分解成多频段的图像，能适应人眼的视觉特性且使得水印的嵌入和检测可分多个层次进行，小波变换域数字水印方法见具有时空域方法和DCT变换域方法的优点。

在一系列信号处理后，如果观察者的主观感觉图像的变化不大，那么图像处理前后低、中频的小波系数的改变幅度同样有限，另外系数幅值改变的方向(变大或变小)在多数情况下也不同，因此，低、中频系数的平均改变幅度十分有限。

本算法选取部分低、中频系数并分成一定大小的系数块，通过量化系数块的平均值来嵌入水印序列。

3.1.1 选择系数

设X(m，n)是一幅大小为M*N灰度级为2“的灰度图像(1≤m≤M ，1≤n≤N)，其中M，N，a为正整数。对X(m，n)进行l层(l为正整数)小波分解，得到3×l个细节图像和一个低频近似图像，用Xk，l(mi，nj)k=h，v ，d；l=1，2，…，l；mi=1，2，…，M／2l ；nj=1，2，…，N／2l 表示选择的小波系数，其中l表示分解的层次，k=h，v，d分别表示第l层水平、垂直和对角方向的子图像。考虑到量化低频子图可能产生较大失真，因此不在其中嵌入水印，而选择除低频外的中频系数[3]。

3.1.2 分块并计算每块的平均值

根据嵌入的信息量和对算法鲁棒性的要求，块越大，水印的鲁棒性越好，但嵌入的水印比特少。把XK，l (mi，nj)分成一定大小的块，用Block (s，t)表示XK，l (mi，nj)中大小为s×t的系数块，其中s=1，2，…，mi ，t=1，2，…，nj，b为正整数，代表该块的编号。其平均值为Ave =∑Block（s，t）/（s*t）：其中 ∑Block 为块内系数幅值的累计和。

3.1.3 量化

水印序列w的嵌入是通过对Ave 的量化完成的，例如：量化成奇数代表嵌入“1”，量化成偶数相当于嵌入“0”。根据对鲁棒性和隐藏性的折中考虑，设量化间隔△l，l=1，2，…，l表示分解层数，对于低频的第l层，由于系数幅值极大，可以作较大间隔的量化，对第l-1，…，1层次作间隔逐渐减小的量化。

根据wi ={0，1}将Ave 量化到与之最近的奇、偶点。用Dat(i，j)表示Block 中的一个小波系数，量化后的该系数用Dat ′(i，j)表示，其中i=1，2，…，s；j=1，2，…，t。

设T=Ave／△l，Turdat=rem([T]，2)其中[]表示四舍五入取整，rem 表示求[T] 除以2的余数。

若Turdat与wi相同，则量化的小波系数为

Dat′(i，j)=Dat(i，j)+[T]×△l-Ave

若Turdat与wi不同，小波系数按下列量化

Dat′(i，j)=Dat(i，j)+（[T]+1）×△l–Ave，T≥[T]

Dat′(i，j)=Dat(i，j)+（[T]-1）×△l –Ave，T<[T]

3.1.4 重构

使用相同的小波基，通过小波逆变换生成含有水印的图像，并将小波基、小波分解层数、选择的系数区域、分块方法、量化间隔、奇偶对应关系记录形成密钥[4]。

3.2 数字水印提取技术

水印的提取由嵌入方式来决定，它是嵌入方式的逆过程。首先对要检测的图像进行小波变换，根据密钥确定嵌入水印的位置，并对水印进行置乱处理的逆运算。可以通过计算归一化互相关函数Nc和峰值信噪比Rpsnr 来度量该水印算法的鲁棒性和透明性。

计算参数:（公式自己查找）

1. 归一化互相关函数Nc

2. 峰值信噪比Rpsnr

二、源代码

clear all;
clc;
key=35;
%Arnold置换次数，作为密钥
Orignalmark=double(imread('suda64.bmp'));  %读入64*64的水印图片
[wrow,wcol]=size(Orignalmark);  
if wrow~=wcol 
    error('wrow~=wcol error');
end
%--- 测试密钥key是否超出范围---------
n=check_arnold(wrow);
if (key+1)>n
    error('arnold key error');
end
Arnoldw=arnold(Orignalmark,wrow,key); %对水印图像进行Arnold转化
[X,fs,bits]=wavread('laile.wav'); %读入音频文件
%X=imread('lena.bmp');
figure;    
subplot(2,1,1); 
%imshow(X),title('原始图像');
%X=double(X);
plot(X);     %显示音频文件波形
title('原始音频信号'); 
%sound(X,fs,bits);
%pause;
%水印嵌入--------------------------------------------------
[c,l]=wavedec(X,2,'db4'); %用db4小波对读入的声音文件进行2级小波分解 
ca2=appcoef(c,l,'db4',2); %提取2级小波分解的低频系数和高频系数 
cd2=detcoef(c,l,2); 
cd1=detcoef(c,l,1); 
lca=length(ca2);  %低频长度
blocksize=fix(lca/(wrow*wcol)); %每块的大小
water_vector=reshape(Arnoldw,1,wrow*wcol);  %将置乱后的水印转化为一维的
wlength=wrow*wcol;  %水印的长度
a=0.25;  %量化步长
j=1;
for i=1:wlength
    Block=ca2(j:j+blocksize-1);
    [U,S,V]=svd(double(Block)); 
    cc=floor(S(1,1)/a);  
    if(Arnoldw(i)==1)           %嵌入奇数倍
        if(mod(cc,2)==0)
            cc=cc+1;
        end
        S(1,1)=a*cc;
    end
    if(Arnoldw(i)==0)            %嵌入偶数倍
        if(mod(cc,2)==1)  
            cc=cc+1;
        end
        S(1,1)=a*cc;
    end
    Blockw=U*S*V';          %SVD 逆变换还原 
    ca2(j:j+blocksize-1)=Blockw; 
    j=j+blocksize;
end
c1=[ca2',cd2',cd1']';
MarkedX=waverec(c1,l,'db4');%b为量化嵌入水印后的音频数据 
subplot(2,1,2); 
plot(MarkedX);
%imshow(MarkedX),title('嵌入水印后图像');
title('嵌入水印后音频'); 
%sound(MarkedX,fs,bits);
wavwrite(MarkedX,fs,bits,'Markedwave.wav'); 
%imwrite(MarkedX,'Markedwave.bmp')
extractedmark=extract(MarkedX,wrow,wcol);
watermark=arnold(extractedmark,wrow,n-key);
figure;
subplot(121);
imshow(Orignalmark),title('原始水印');
subplot(122);
imshow(watermark),title('提取水印');
psnrvalue=psnr(Orignalmark,watermark)
Nc=nc(Orignalmark,watermark)
function wa=arnold(w,wrow,key)
wa=w;
for k=1:key
    for x=1:wrow
        for y=1:wrow
             x1=x+y;
             y1=x+2*y;
             if x1>wrow
                 x1=mod(x1,wrow);
             end
             if y1>wrow
                 y1=mod(y1,wrow);
             end
             if x1==0
                 x1=wrow;
             end
             if y1==0
                 y1=wrow;
             end
             wa(x1,y1)=w(x,y);
        end
    end
     w=wa;
end


  

 

  
1
  
2
  
3
  
4
  
5
  
6
  
7
  
8
  
9
  
10
  
11
  
12
  
13
  
14
  
15
  
16
  
17
  
18
  
19
  
20
  
21
  
22
  
23
  
24
  
25
  
26
  
27
  
28
  
29
  
30
  
31
  
32
  
33
  
34
  
35
  
36
  
37
  
38
  
39
  
40
  
41
  
42
  
43
  
44
  
45
  
46
  
47
  
48
  
49
  
50
  
51
  
52
  
53
  
54
  
55
  
56
  
57
  
58
  
59
  
60
  
61
  
62
  
63
  
64
  
65
  
66
  
67
  
68
  
69
  
70
  
71
  
72
  
73
  
74
  
75
  
76
  
77
  
78
  
79
  
80
  
81
  
82
  
83
  
84
  
85
  
86
  
87
  
88
  
89
  
90
  
91
  
92
  
93
  
94
  
95
  
96
  
97
  
98
  
99
 

三、运行结果

四、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1]韩纪庆,张磊,郑铁然.语音信号处理（第3版）[M].清华大学出版社，2019.
[2]柳若边.深度学习:语音识别技术实践[M].清华大学出版社，2019.

文章来源: qq912100926.blog.csdn.net，作者：海神之光，版权归原作者所有，如需转载，请联系作者。

原文链接：qq912100926.blog.csdn.net/article/details/114743860