第11章 小波分析用于信号压缩

信号压缩是信号处理中一类重要的课题之一。MATLAB小波工具箱提供了一种通过对小波域系数进行处理实现压缩的有损压缩算法。
信号压缩的动机源于对模拟信号的数字化。在离散化的过程中，大部分应用是使用统一的采样率，这样就引入很大的信息冗余，以声音信号为例，其携带大部分的信息都集中在低频部分，而一般电话的采样率都是统一的，所以在很多时间使用的采样率高于需要的精度，造成了信息的冗余。
本章从实践角度出发，首先介绍压缩的基本原理与方法；其次，介绍有关MATLAB小波工具箱中的压缩函数，并且将不惜篇幅地列举压缩的算例供读者参考；最后，介绍语音信号的压缩原理与实例。
学习目标：
(1)了解小波压缩的原理和方法
(2)熟练掌握MATLAB中小波压缩函数
(3)熟练掌握一维小波工具箱对信号压缩
(4)掌握语音信号压缩的方法
11.1 信号压缩
与信号去噪相似，压缩领域中由于小波的特殊优点，应用其进行压缩也受到了许多学者的重视，并获得了非常好的效果。
11.1.1 小波压缩概述
随着数据压缩理论的不断发展和日益成熟，一些研究人员逐渐开始对声音、文字和信号等信号的压缩技术进行研究，先后经历了经典压缩方法和现代压缩方法两个阶段。
经典压缩方法遵循香农信息论的基本理论，可分为无损压缩和有损压缩两类压缩方法。
(1)无损压缩方法(如Huffman编码、算术编码、游程编码等方法)设法改变信源的概率分布，使信号的概率分布尽可能地非均匀，再用最佳编码方法重新对每个信号分配码字，使平均码长逼近信源熵。该类方法的压缩效率都以其信息熵为上限。
(2)有损压缩方法(如预测编码、变换域编码、混合编码、矢量量化等方法)设法去除信源之间的相关性，使之成为或差不多成为不相关的信源，该类压缩方法也受信息熵的约束。
随着小波变换(Wavelet Transform，WT)、分形几何理论、数学形态学等数学理论和相关学科(如模式识别、人工智能、神经网络、感知生理心理学等)的深入发展，相继出现了新颖高效的现代压缩方法。该方法包括子带编码、小波变换编码、神经网络编码、分形编码、模型基编码等。随着信号压缩技术在相关领域中的广泛应用，信号压缩的国际标准相继被制定，如在静止信号压缩方面已经制定了 CCIT T.81、ISO 10918(JPEG)及JPEG2000等标准。
从信号分析角度来讲，傅立叶变换是信号和数字信号处理的理论基础，小波分析将信号和数字信号处理带入到一个崭新的领域。
多尺度几何分析是继小波分析后的新一代信号分析工具，它具有多分辨、局部化和多方向性等优良特性，更适合于处理高维信号。
为了缓解人们在处理信息的巨量需求所造成的信号采样、传输和存储的巨大压力，近年来许多小波分析方面的权威专家又提出了压缩感知理论。与 Nyquist 采样定理不同，它指出，只要信号在某个变换域是稀疏的，那么就可以用一个与变换基不相关的观测矩阵将变换所得高维信号投影到一个低维空间上，然后通过求解一个优化问题就可以从这些少量的投影中以高概率重构出原信号，可以证明这样的投影包含了重构信号的足够信息。
在压缩感知理论中，信号与信号的采样和压缩同时以低速率进行，使传感器的采样和计算成本大大降低，而信号的恢复过程是一个优化计算的过程。因此，该理论指出了将模拟信号直接采样压缩为数字形式的有效途径，具有直接信息采样特性。
由于从理论上讲任何信号都具有可压缩性，只要能找到其相应的稀疏表示空间，就可以有效地进行压缩采样，这一理论必将给信号采样方法带来一次新的革命。压缩感知理论的引人之处还在于它对应用科学和工程的许多领域具有重要的影响和实践意义，如统计学、信息论、编码等。
目前，研究此理论的的机构越来越多，美国Rice大学成功地研制出了单相素CS相机，它可利用单一的信号光子检测器采样得到比信号像素点数少得多的点，恢复得到信号，并且具有对信号波长自适应的能力，这种自适应能力是传统的CCD和CMOS 成像器件所不具备的。
11.1.2 一维小波分析进行压缩的原理
信号压缩的基本目标就是在不损失信号所携带的信息的前提下，尽可能地减少用于存储信号的开销。用于信号压缩的方法有很多，这里只介绍通过小波变换实现信号压缩的方法、其基本思想如下。
小波变换实现数据压缩的基本目标是使信号在时间—频率域的分解系数所占的存储空间尽可能小，同时还要保证压缩后的系数能精确地反映原信号所携带的信息。
一维小波分析进行信号压缩的一般步骤如下。
(1)信号的小波分解。
(2)将高频系数进行阈值量化处理。对从1到N的每一层高频系数都可以选择不同的阈值，并且用硬阈值进行系数的量化。
(3)对量化后的系数实行小波重构。
对比上一章中提到的小波阈值去噪方法，信号的压缩与去噪不同之处仅在第2步。
一般来说，有两种比较有效的信号压缩方法：第一种方法是对信号进行小波尺度的扩展，并且保留绝对值最大的系数，在这种情况下，可以用全局阈值来压缩信号，实现信号的压缩或相对均方误差规范的信号恢复，这时只需用一个输入参数；第二种方法是根据每一层分解后的效果来决定该层的阈值，而且每层的阈值是不同的。
类似于信号去噪，信号压缩也从两方面来度量某种压缩算法的性能，即压缩比和误差。在一般情况下，这两者在本质上是矛盾的，这就需要研究者在其间做比较。
R.A.DeVore提出了一个度量压缩性能的数学模型。设原始信号为X，压缩后的信号为
 ，则误差函数为：
在使用的压缩方法保留的系数不超过 N的情况下，err 度量了该压缩算法的误差。在这个意义下，我们通过eff的衰减速度来定义压缩算法的性能。如果eff满足：
即随着N的增加，eff以N的多项式的速度衰减，那么一个算法如果包含了其他算法的α类，则说明前一个算法较优。
11.1.3 小波压缩实现方法
在小波变换中，信号压缩可以归结为在信号的小波分解域对小波系数进行量化的过程。由于量化的过程可以看作是系数的内部扰动，所以在小波域中系数对扰动的稳定性是至关重要的。对于实现来讲，一个首要的问题就是用更少的位来存储小波系数并尽量保证重构信号的质量。
为了实现信息的分离，减少存储或传输的代价，压缩后用以存储或传输的数据应该是小波域经过量化以后的数据。一种简单的量化方法叫做归一化量化方法，其基本思想是把整个系数的值域分解为不同的区间，各区间所有值出现的概率之和相同，且对每个区间的值用相同的位数加以量化。
有了这些讨论，我们可以看一下信号压缩的具体流程，如图11-1所示。
图11-1 小波变换用于数据压缩的流程
信号f：被压缩的原始信号在小波分解域中是时频一致的。
零填充：是去除采样过程中的时间冗余，对采样率为N0
 的情况，作用如下算式：
编码器：对输入系数进行编码，去除其编码冗余，并根据需要进行量化。
存储介质：编码数据的存储容器。
解码器：将编码后的系数重新恢复到小波表示域，便于重构。
重构信号：恢复压缩后的信号。
11.2 MATLAB压缩函数
基于上一节介绍的原理，MATLAB小波工具箱主要针对阈值选取开发出了压缩函数，本节将作以介绍。
11.2.1 一维小波分析进行信号压缩
尽管同去噪过程有很大区别，但是压缩算法也是通过在信号的小波域表示作用阈值来实现的，因此其MATLAB函数大部分与去噪过程相同。这里先介绍压缩函数wdencmp。
函数wdencmp，其调用格式为：
[XC，CXC，LXC，PERF0，PERFL2]=
wdencmp('gbl'，X，'wname'，N，THR，SORH，KEEPAPP)
[XC，CXC，LXC，PERF0，PERFL2]=wdencmp('lvd'，X，'wname'，N，THR，SORH)
[XC，CXC，LXC，PERF0，PERFL2]=wdencmp('lvd'，C，L，'wname'，N，THR，SORH)
其中：gbl代表阈值是全局阈值；lvd代表阈值是在各层大小不同的数值；X是输入信号；wname是选取的小波基；N是分解层数；KEEPAPP=1 时将小波分解的低频系数不作任何处理；KEEPAPP=0时对小波分解的低频系数也进行阈值量化处理；XC返回的是经过原信号X进行压缩处理后的信号；CXC返回的是XC结构；LXC返回的是XC结构。
下面给出一个具体的例子，以便使读者对小波分析的这个重要应用有一个直观的认识。
【例11.1】现在利用小波分析，对给定信号进行压缩处理。MATLAB代码设置如下：
clc;
clear;
%加载信号
load leleccum;
%采样点
s=leleccum(2000:2800);
%画出原始信号
plot(s);
title('原始信号');
结果如图11-2所示。
%用db4小波对原始信号进行3层分解
[c，l]=wavedec(s，3，'db4');
%指定不同阈值对系数进行处理
sorh='h';
keepapp=1;
thr=33;
[s1，CXC，LXC，PERF01，PERFL21]=wdencmp('gbl'，c，l，'db4'，3，thr，sorh，keepapp);
thr=23;
[s2，CXC，LXC，PERF02，PERFL22]=wdencmp('gbl'，c，l，'db4'，3，thr，sorh，keepapp);
thr=13;
[s3，CXC，LXC，PERF03，PERFL23]=wdencmp('gbl'，c，l，'db4'，3，thr，sorh，keepapp);
%输出结果
figure(2)
subplot(3，1，1);
plot(s1);
grid;
title('不同阈值对信号进行压缩处理');
xlabel(['能量剩余'，num2str(PERFL21)，'零系数成分'，num2str(PERF01)]);
subplot(3，1，2);
plot(s2);
xlabel(['能量剩余'，num2str(PERFL22)，'零系数成分'，num2str(PERF02)]);
subplot(3，1，3);
plot(s3);
xlabel(['能量剩余'，num2str(PERFL23)，'零系数成分'，num2str(PERF03)]);
结果如图11-3所示。
图11-2 原始信号
图11-3 结果图
从图11-3可以看出，根据阈值不同，可以很灵活地根据需要确定应用所需要的压缩性能。
从此例还可以看出小波变换用于数据压缩的基本思想，就是将信号通过小波变换把低采样率的近似系数和少量需要保留的近似系数组合在一起，而信号其他的细节系数(不是信号的主要性质)，都通过阈值滤掉，最后完成压缩。
在本例中，阈值的作用方式使用的是全局阈值，可能对信号的局部性质刻画的不够精细，读者可以自己设计分层阈值的例子，并对压缩的性能与全局阈值方法的结果进行比较。
11.2.2 小波包分析进行信号压缩
在前面的章节中已经详细介绍了小波包分析，其信号压缩的思想和在小波分析中基本相同，不同之处在于：小波包提供了一种更为复杂、更为灵活的分析手段，因为小波包分析对上一层的低频部分和高频部分同时实行分解，具有更加精确的局部分析能力。
对信号进行小波包分解时，可以采用多种小波包基。通常根据分析信号的要求，从中选择最好的一种小波包基，即最优基。最优基的选择标准是熵标准。在MATLAB的小波工具箱中可通过besttree函数进行最优基的选择。
应用小波包分析对信号进行压缩处理是一个最基本的功能。小波包阈值压缩过程主要分为4个步骤，如下。
(1)信号的小波包分解。选择一种小波包基并确定所需分解的层次，然后对信号进行小波包分解。
(2)最优小波包基的选择。
(3)小波包分解系数的阈值化。对于每一个小波包分解系数，选择一个恰当的阈值对小波包分解后的系数进行阈值量化处理。
(4)信号的小波包重构。对低频系数和经过处理后的高频系数进行小波包重构。
利用小波包进行压缩首先要选取小波包基和分解的层数。对称性好的小波不产生相位畸变，正则性好的小波易于获得光滑的重构信号。
【例11.2】利用小波包分析对一个给定的信号进行压缩处理。MATLAB代码设置如下：
首先，加载信号：
clear
load noisdopp;
s=noisdopp(500:1000);
plot(s)
title('原始信号')
结果如图11-4所示。
[thr，sorh，deepapp，crit]=ddencmp(‘cmp’，'wp'，s);
[x，wpt，perf0，perfl2]=...
wpdencmp(s，sorh，3，'db9'，crit，thr，deepapp);
figure(2)
plot(x)
title('压缩信号1')
结果如图11-5所示。
图11-4 原始信号
图11-5 三层分解压缩
clear t x ;
[thr，sorh，deepapp，crit]=ddencmp(‘cmp’，'wp'，s);
[x，treed，perf0，perfl2]=...
wpdencmp(s，sorh，4，'db9'，crit，thr，deepapp);
figure(3)
plot(x)
title('压缩信号2')
结果如图11-6所示。
clear t x ;
[thr，sorh，deepapp，crit]=ddencmp(‘cmp’，'wp'，s);
[x，treed，perf0，perfl2]=...
wpdencmp(s，sorh，5，'db9'，crit，thr，deepapp);
figure(4)
plot(x)
title('压缩信号3')
结果如图11-7所示。
图11-6 四层分解压缩
图11-7 五层分解压缩
clear t x ;
[thr，sorh，deepapp，crit]=ddencmp(‘cmp’，'wp'，s);
[x，treed，perf0，perfl2]=...
wpdencmp(s，sorh，8，'db9'，crit，thr，deepapp);
figure(5)
plot(x)
title('压缩信号4')
结果如图11-8所示。
说明：本例选择对称性和正则性均较好的db9小波，并且分别选用3、4、5和8层小波包分解。
小波包分解层数是压缩处理过程的一个关键参数，通过上面的程序可以确定。从得到的不同分解层次后的压缩效果图可以看出，随着小波包分解层次的增加，压缩效果变好，但是分解层次增加到4层以后，压缩效果改善已经不明显，此时反而增加了计算量。
图11-8 八层分解压缩
11.3 MATLAB一维小波工具箱压缩
前面的章节曾经向读者介绍过如何利用小波工具箱的界面来进行对信号去噪的操作，本节主要介绍关于信号压缩的工具箱。
11.3.1 一维离散小波界面式压缩
本小节将向读者介绍利用小波工具箱的一维离散小波变换去除信号的噪声。
打开一维小波分析工具。首先，在窗口中键入wavemenu，弹出小波工具箱主界面。下面向读者介绍一维离散小波分析工具分析真实世界噪声信号的具体操作步骤，具体操作如下
(1)打开一维小波分析工具。首先，在窗口中键入wavemenu，弹出小波工具箱主界面，然后单击Wavelet 1-D按钮。
(2)下载信号源。选择File→Load→Signal命令，在弹出的对话框中选择noissin.mat文件，它的路径为toolbox/wavelet/wavedemo，单击OK按钮。这样信号就被下载到了此工具中，如图11-9所示。
图11-9 信号源
(3)执行3级小波系数分解。在如图11-9所示的界面中选择db4小波进行3级分解。同时单击Analyze按钮，经过短暂的计算后，工具显示分解系数图，如图11-10所示。
图11-10 分解系数图
(4)执行多级小波系数分解。在界面的右上方，从at levels下拉列表中选择3，并单击Analyze按钮，如图11-10所示。界面的右下方提供了6种分解模式选项，如图11-11所示。读者可以通过选择不同的选项来观看不同分解系数模式。
图11-11 分解模式选项
(5)压缩。单击Compress按钮(在Analyze按钮的下面)，一维小波压缩窗出现，如图11-12所示。
(6)选择阈值。在图11-12 右上方单击By level thresholdting按钮，如图11-13所示。显示每分解系数的阈值区间。
(7)信号压缩。单击View Compress Signal按钮，弹出压缩窗口，如图11-14所示。
图11-12 压缩图
图11-13 阈值区间
图11-14 压缩信号
(8)显示统计值。单击Residuals按钮，弹出压缩后的信号与原始信号的误差统计值，如图11-15所示。
图11-15 误差统计值
11.3.2 一维小波包界面式压缩
本小节将采用一维小波包图形工具来对一个噪声信号进行压缩。
(1)启动小波工具箱。方法同上一小节一样。
(2)启动一维小波包的主界面。单击Wavelet Packets 1-D按钮。
(3)下载信号源。选择 File→Load Signal命令，在弹出的对话框中选择 noissin.mat文件，它的路径为toolbox/wavelet/wavedemo，单击OK按钮。这样信号就被下载到了此工具中，如图11-16所示。计算时采用shannon熵标准。
图11-16 加载信号
(4)小波包分析。在图11-16的右边可以对信号进行设置。选择db4小波，并且分解层数设置为3，选择的熵类型为shannon。单击Analyze按钮，执行对信号分解，显示在窗口中，如图11-17所示。
图11-17 信号分析
(5)计算最佳树。在图11-17中的右上方单击Best Tree按钮，得到如图11-18所示的最佳树。
(6)选择阈值进行压缩。单击Compress按钮，弹出如图11-19所示的一维小波包压缩窗口。
图11-18 最佳树
图11-19 压缩
(7)选择Remove near 0选项，去掉附近的零，如图11-20所示。
图11-20 去零
(8)单击View Compress Signal按钮可以看到压缩信号，如图11-21所示。
图11-21 压缩信号图
说明：选中图11-21中的Origianl Signal复选框即可看出压缩信号(黑色)和原始信号(红色)在同一个显示框中。
(9)显示统计值。单击Residuals按钮，弹出压缩后的信号与原始信号的误差统计值，如图11-22所示。
通过MATLAB小波工具箱来进行信号压缩操作简单方便。同时，在这一节中读者可以对一种噪声信号通过小波变换和小波包变换两种小波方法来进行压缩处理，并可以比较发现两种方法的优劣。
图11-22 误差统计值
11.4 小波压缩综合实例
本节将通过实例来向读者演示小波方法压缩的算例，希望可以对这方面感兴趣的研究人员提供一些编程方法。
语音信号作为一维非平稳信号，含有很多尖峰和突变部分，高频信号占的比例很大，对其进行小波分解后每一层的幅值与白噪声容易出现混叠。
下面利用小波变换进行语音压缩。
【例11.3】小波变换进行语音压缩。
MATLAB代码设置如下：
%小波压缩MATLAB程序
clear
clc
%在噪声环境下语音信号的增强
%语音信号为读入的声音文件
%噪声为正态随机噪声
load handel.mat
hfile='handel.wav';
wavwrite(y，Fs，hfile)
plot(y);
title('原始语音信号');
结果如图11-23所示。
[c，l]=wavedec(y，3，'db1');
cA3=appcoef(c，l，'db1'，3);
[cD1，cD2，cD3]=detcoef(c，l，[1 2 3]);
figure(2)
subplot(411)
plot(cA3);
title('分解系数');
subplot(412);
plot(cD1);
subplot(413);
plot(cD2);
subplot(414);
plot(cD3);
结果如图11-24所示。
图11-23 原始信号
图11-24 分解系数
%quantize
[partition，codebook]=lloyds(cA3，2^4);
[index，quants]=quantiz(cA3，partition，codebook);
[partition1，codebook1]=lloyds(cD1，2^4);
[index1，quants1]=quantiz(cD1，partition1，codebook1);
[partition2，codebook2]=lloyds(cD2，2^4);
[index2，quants2]=quantiz(cD2，partition2，codebook2);
[partition3，codebook3]=lloyds(cD3，2^4);
[index3，quants3]=quantiz(cD3，partition3，codebook3);
%probability
g=(1/length(codebook))*ones(1，length(codebook));
g1=g;
g2=g;
g3=g;
%huffman编解码
[dict，avglen]=huffmandict(codebook，g);
comp=huffmanenco(quants，dict);% Encode the data.
dsig=huffmandeco(comp，dict);
isequal(quants，dsig)
%huffman编解码
[dict1，avglen1]=huffmandict(codebook1，g1);
comp1=huffmanenco(quants1，dict1);% Encode the data.
dsig1=huffmandeco(comp1，dict1);
isequal(quants1，dsig1)
%huffman编解码
[dict2，avglen2]=huffmandict(codebook2，g2);
comp2=huffmanenco(quants2，dict2);% Encode the data.
dsig2=huffmandeco(comp2，dict2);
isequal(quants2，dsig2)
%huffman编解码
[dict3，avglen3]=huffmandict(codebook3，g3);
comp3=huffmanenco(quants3，dict3);% Encode the data.
dsig3=huffmandeco(comp3，dict3);
isequal(quants3，dsig3)
%dsig3=dsig3';
%dsig2=dsig2';
%dsig1=dsig1';
%dsig=dsig';
%重构
cA2=idwt(dsig，dsig3，'db1');
icA2(1:13950)=cA2(1:13950);
cA1=idwt(icA2，dsig2(1:13950)，'db1');
icA1(1:13950)=cA1(1:13950);
XX=idwt(icA1，dsig1(1:13950)，'db1');
XX=XX';
figure(3)
plot(XX);
title('重构信号');
XX=XX(1:27000);
MSE=sum((y(1:27000)-XX(1:27000)).^2)/length(XX);
PSNR=10*log10((max(max(XX)))^2/MSE)
结果如图11-25所示。
注：为了让读者对语音信号有感观上的认识，我们将语音信号转化为数字信号处理，因此在运行此程序时会产生后缀名为wav的音频文件。
图11-25 重构信号
11.5 本章小结
随着小波变换(Wavelet Transform，WT)、分形几何理论、数学形态学等数学理论和相关学科(如模式识别、人工智能、神经网络、感知生理心理学等)的深入发展，相继出现了新颖高效的现代压缩方法。特别是，小波变换编码已经成为国际标准编码。对于应用小波来处理信号压缩与去噪既有区别又有联系，但方法类似，因此本章承接上一章的结构。
同时，基于小波变换对信号压缩效果较好的基础上，介绍了信号压缩的小波方法以及MATLAB压缩函数和算例，最后补充介绍了基于小波变换的语音信号压缩。希望对此感兴趣的爱好者可以从中有所收获！