解读matlab之小波库函数,解读matlab之小波库函数.doc

解读matlab之小波库函数

南京理工大学仪器科学与技术专业 谭彩铭

2010-4-2

使用的matlab软件版本为matlab7.1

1 dwt函数

dwt函数是单尺度一维小波变换函数。

dwt函数执行过程中调用了函数conv2，这个函数是运算的关键，需要首先明白conv2函数的执行过程。要明白conv2函数，需要先明白conv函数。

对w = conv(u,v)运算

Let m = length(u) and n = length(v). Then w is the vector of length m+n-1 whose kth element is

式(1)

假设h=[h(1) h(2) h(3) h(4)]，x=[x(1) x(2) x(3) x(4) x(5) x(6) x(7)]，为更直接地表达y=conv(h,x)的计算过程，作如下示意图。其中length(y)=7+4-1。

图1

对c=conv2(a,b)运算

这里，a和b为一维或二维矩阵，其计算过程可由下式表示

式(2)

The size of c in each dimension is equal to the sum of the corresponding dimensions of the input matrices, minus one. That is, if the size of a is [ma,na] and the size of b is [mb,nb], then the size of C is [ma+mb-1,na+nb-1].

其计算过程可以由下表表示

表1

c(1,:)

conv(a(1,:),b(1,:))

c(2,:)

conv(a(1,:),b(2,:))+ conv(a(2,:),b(1,:))

c(3,:)

conv(a(1,:),b(3,:))+ conv(a(2,:),b(2,:)) +conv(a(3,:),b(1,:))

……

……

下面研究一下conv2函数中的‘valid’参数的用法。

The formula is c=conv2(a,b,’valid’)

Valid：Returns only those parts of the convolution that are computed without the zero-padded edges. Using this option, c has size [ma-mb+1,na-nb+1] when all(size(a) >= size(b)). Otherwise conv2 returns [].

假设h=[h(1) h(2) h(3) h(4)]，x=[x(1) x(2) x(3) x(4) x(5) x(6) x(7)],计算y=conv2(x,h,’valid’)，它的计算过程可由下图表示，其中length(y)= 7-4+1。

图2

下面来看看dwt函数的工作过程

假设x=[x(1) x(2) x(3) x(4) x(5) x(6) x(7)],计算y=dwt(x,’db2’)。

其计算过程主要由两个部分组成：

第一部分：边缘延拓，它主要由函数wextend完成。

第二部分：卷积运算，它主要由函数conv2完成。

先看第一部分，仔细分析子程序部分，函数wextend的用法为y=wextend(1D,sym,x,3);

这样得到的y=[ x(3) x(2) x(1) x(1) x(2) x(3) x(4) x(5) x(6) x(7) x(7) x(6) x(5)]

在看第二部分，仔细分析子程序部分，核心语句有z=conv2(y,Lo_D,valid);

这里设Lo_D=[h(1) h(2) h(3) h(4)]。

结合图2所示对conv2函数用法的介绍，绘制下图表示该处卷积的计算过程。

图3

由此可见，小波库函数dwt中是如何处理边缘点的。

最后就是下采样，其下采样是按照式a = z(2:2:length(z))进行的，高频低频部分均如此，项数为floor((7+4-1)/2)。

2 idwt函数

idwt函数是单尺度一维离散小波逆变换函数。

小波重构函数的核心函数为上采样函数dyadup和卷积运算函数conv2。

下面先简要看一看dyadup函数。

dyadup implements a simple zero-padding scheme very useful in the wavelet reconstruction algorithm.

Y = dyadup(X,EVENODD), where X is a vector, returns an extended copy of vector X obtained by inserting zeros. Whether the zeros are inserted as even- or odd-indexed elements of Y depends on the value of positive integer EVENODD: If EVENODD is even, then Y(2k-1) = X(k), Y(2k) = 0. If EVENODD is odd, then Y(2k-1) = 0, Y(2k) = X(k).

Y = dyadup(X) is equivalent to Y = dyadup(X,1) (odd-indexed samples).

设x=[1 2 3 4 5 6 7],该函数的运算过程可由下表较直接地表示。对一维的情况，就表中所述几种。

表2

dyadup(x)

[0 1 0 2 0 3 0 4 0 5 0 6 0 7 0]

dyadup(x,1)

[0 1 0 2 0 3 0 4 0 5 0 6 0 7 0]

dyadup(x,3)

[0 1 0 2 0 3 0 4 0 5 0 6 0 7 0]

dyadup(x,0)

[1 0 2 0 3 0 4 0 5 0 6 0 7]

dyadup(x,0,1)

[1 0 2 0 3 0 4 0 5 0 6 0 7 0]

dyadup(x,2)

[1 0 2 0 3 0 4 0 5 0 6 0 7]

至于conv2函数，之前已作介绍。这里要关心的是边界点的处理问题。

idwt函数中，对于数值的取舍问题用到的函数为wkeep1，下面先研究下函数wkeep1。(好像在help中搜寻不到，可在命令窗口中输入help wkeep1命令，得到下列使用说明)

Y = WKEEP1(X,L,OPT) extracts the vector Y

from the vector X. The length of Y is L.

If OPT = c (l , r, respectively), Y is the central(left, right, respectively) part of X.

Y = WKEEP1(X,L,FIRST) returns the vector X(FIRST:FIRST+L-1).

Y = WKEEP1(X,L) is equivalent to Y = WKEEP1(X,L,c).

下表中列出了iwt函数中用到的两种情况。

设x=[1 2 3 4 5 6 7]。

表3

wkeep1(x,3,c)

[3 4 5]

wkeep1(x,2,c)

[3 4]

下面仔细分析一下边缘点的处理问题。之前在文档中对完全重构滤波器的分解与重构过程已作分析。下列根据理论自己编写的程序便能完成完全重构过程。

clear;

load noissin;

v=noissin(1:6);

[ca1,cd1,tip1]=funbreakupindb2(v);

c=funreconstructindb2(ca1,cd1,tip1);

plot(v);

hold on;

plot(c,r.);

hold off;

图4

function [ca1,cd1,tip]=funbreakupindb2(v)

[h0,h1,h2,h3] = wfilters(db2);

w0=conv(h0,v);

x0=conv(h1,v);

N2=length(w0);

tip=mod(N2,2);

if(tip == 1)

w0(N2+1)=0;

x0(N2+1)=0;

end

N=floor(length(w0)/2);

w0=reshape(w0,2,N);

x0=reshape(x0,2,N);

w1=w0(1,:);

x1=x0(1,:);

ca1=w1;

cd1=x1;

图5

function c=funreconstructindb2(ca1,cd1,tip)

[h0,h1,h2,h3] = wfilters(db2);

N=length(ca1);

w1=ca1;

x1=cd1;

w2=[w1;zeros(1,N)];

w3=w2(:);

x2=[x1;zeros(1,N)];

x3=x2(:);

y=conv(h2,w3)+conv(h3,x3);

if(tip==1)

N1=2*N-4;

else

N1=2*N-3;

end

c=y(4:N1+3);

图6

图7

但是你会发现边缘处的尺度函数系数和小波函数系数偏离的很远，这在阈值滤波时便会在边缘点处产生很大的误差。因为边缘点处没有完全用到4个滤波系数所致。

matlab小波库函数dwt中对这个问题的解决方案是边缘延拓，之前已作介绍，如图3所示。这里的描述借助图3的描述，依然假设x=[x(1) x(2) x(3) x(4) x(5) x(6) x(7)],计算y=dwt(x,’db2’)。作边缘延拓后，得到图3所示结果。如图3所示，作卷积运算后，得到的项数共有16项，记为v(1)~v(16),但是作边缘取舍和下采样后得到项为z(2)，z(4)，z(6)，z(8)，z(10)。

接着的是重构过程，重构过程也是作卷积的过程，我们先对v(1)~v(16)项作下采样并作重构卷积运算，就可以清晰地知道z项的重构卷积运算是否达到重构效果。下图显示了重构过程。

图8

这里输入序列x的项数为奇数7，对于偶数项有稍许不同。

下面设x=[x(1) x(2) x(3) x(4) x(5) x(6)]，作边缘延拓并作卷积运算得到的项数为15项，记为v(1)~v(15)，再作边缘取舍和下采样后得到项为z(2)，z(4)，z(6)，z(8)。下图显示了重构过程。

图9

从而得证。

下图程序是dwt和idwt函数的简单应用。

clear;

load noissin;

s=noissin(1:6);

[ca1,cd1] = dwt(s,db2);

ss = idwt(ca1,cd1,db2);

plot(s);

hold on;

plot(ss,r.);

hold off;

图10

图11

对图10所示程序，如果输入数据是奇数，重构后的点或多一个，不过最后两个点值一样，这也许是idwt函数考虑得不全面所致。如下图所示。

clear;

load noissin;

s=noissin(1:7);

[ca1,cd1] = dwt(s,db2);

ss = idwt(ca1,cd1,db2);

plot(s);

hold on;

plot(ss,r.);

hold off;

图12

图13

3 wavedec函数

wavedec函数是多尺度一维小波离散分解函数。

wavedec performs a multilevel one-dimensional wavelet analysis using either a specific wavelet (wname) or a specific wavelet decomposition filters (Lo_D and Hi_D, see wfilters).

[C,L] = wavedec(X,N,wname) returns the wavelet decomposition of the signal X at level N, using wname. N must be a strictly positive integer (see wmaxlev for more information). The output decomposition structure contains the wavelet decomposition vector C and the bookkeeping vector L. The structure is organized as in this level-3 decomposition example:

图14

wavedec函数中的核心函数为dwt函数，作多少尺度的分解就调用几次dwt函数。

4 waverec函数

waverec函数是多尺度一维离散小波重构函数。

waverec performs a multilevel one-dimensional wavelet reconstruction using either a specific wavelet (wname, see wfilters) or specific reconstruction filters (Lo_R and Hi_R). waverec is the inverse function of wavedec in the sense that the abstract statement waverec(wavedec(X,N,wname),wname) returns X.

X = waverec(C,L,wname) reconstructs the signal X based on the multilevel wavelet decomposition structure [C,L] and wavelet wname. (For information about the decomposition structure, see wavedec.)

X = waverec(C,L,Lo_R,Hi_R) reconstructs the signal X as above, using the reconstruction filters you specify. Lo_R is the reconstruction low-pass filter and Hi_R is the reconstruction high-pass filter.

Note that X = waverec(C,L,wname) is equivalent to X = appcoef(C,L,wname,0).

waverec函数中的核心函数是appcoef函数，故先要研究下appcoef函数函数的工作过程。

5 appcoef函数

appcoef函数为提取一维离散小波近似分量函数。

appcoef is a one-dimensional wavelet analysis function.

appcoef computes the approximation coefficients of a one-dimensional signal.

A = appcoef(C,L,wname,N) computes the approximation coefficients at level N using the wavelet decomposition structure [C,L] (see wavedec for more information).

wname is a string containing the wavelet name. Level N must be an integer such that 0≤N≤length(L)-2.

A = appcoef(C,L,wname) extracts the approximation coefficients at the last level: length(L)-2.

Instead of giving the wavelet name, you can give the filters.

For A = appcoef(C,L,Lo_R,Hi_R) or A = appcoef(C,L,Lo_R,Hi_R,N), Lo_R is the reconstruction low-pass filter and Hi_R is the reconstruction high-pass filter (see wfilters for more information).

因为小波重构的过程实际上就是得到‘0尺度’下的近似分量。故waverec函数可以说是appcoef函数功能的一部分。

appcoef函数中的两个主要函数是detcoef函数和idwt函数。

下面先看一下detcoef函数。

detcoef函数为提取一维离散小波细节分量函数。我们看一下wavedec函数，wavedec函数的返回值为c,l。 detcoef函数是为了更方便地取出c中的细节分量数据，实现过程简单。

appcoef函数可以得到任意尺度下的近似分量数据。

下图是wavedec函数和waverec函数的简单应用程序。

clear;

load noissin;

x=noissin(1:7);

[c,l]=wavedec(x,1,db2)

x1= clear;

load noissin;

s=noissin(1:7);

[ca1,cd1] = dwt(s,db2);

ss = idwt(ca1,cd1,db2);

plot(s);

hold on;

plot(ss,r.);

hold off;

(c,l,db2)

plot(x,.-);

hold on;

plot(x1,.r);

hold off;

图15

图16

观察图13，它出现了多出一个输出点的情况。而图16中就不存在了。原因简单，idwt函数返回的细节分量和近似分量的项数均为floor((N+4-1)/2)(设小波为db2小波，N为原始数据长度)。易知N=7和N=8两种情况下的项数均为5，这时如图12中所示程序调用idwt函数ss=idwt(ca1,cd1,db2)，idwt显然无法判别原始数据长度。故会出现多一点的情况。实际上前人早已考虑到这一点，idwt函数中可以加入数据长度这个参数。在图15中的wavedec函数内部调用idwt函数时即考虑了这个问题。

6 upwlev函数

upwlev函数是单尺度一维离散小波分析的重构函数。该函数是在wavedec函数的分解作用下向上一层重构。没有多大作用。下列程序可以比较直观地说明它的作用。

load sumsin; s = sumsin;

[c,l] = wavedec(s,3,db1);

subplot(311); plot(s);

subplot(312); plot(c);

[nc,nl] = upwlev(c,l,db1);

subplot(313); plot(nc);

图17

图18

7 wrcoef函数

wrcoef函数用于对一维小波变换进行单支重构，这包括两种情况。

第一种情况：令细节分量为0，进行重构。

第二种情况：令近似分量为0，进行重构。

两种情况均可以选择在那个尺度下开始。此函数貌似用处也不大。

8 upcoef函数

该函数可用于得到各个尺度下的尺度函数和小波函数。

9 ddencmp函数

自动生成小波消噪或压缩处理的阈值处理方案。

调用方式：[THR,SORH,KEEPAPP,CRIT] = ddencmp(IN1,IN2,X)自动生成信号x的小波(或小波包)消噪或压缩的阈值选取方案。

输入参数x为一维或二维的信号向量或矩阵。

输入参数IN1指定处理的目的是压缩还是消噪，可选值为IN1=’den’，为信号消噪；IN1=’cmp’，为信号压缩。

输入参数IN2可选值为IN2=’wv’,使用小波分解；IN2=’wp’,使用小波包分解。

输出参数THR为函数选择的阈值。

输出参数SORH为函数选择的阈值使用方式，可选值为SORH=s，为软阈值；SORH=h，为硬阈值。

输出参数CRIT为使用小波包分解时选取的熵函数类型。

下面以下图所示程序为例分析这个函数。

clear;

load noissin;

x=noissin;

[thr,sorh,keepapp] = ddencmp(den,wv,x)

图19

仔细分析ddencmp函数的代码(if isequal(dorc,den) & isequal(worwp,wv) , sorh = s; else ,sorh = h; end)，发现只要输入参数为’den’和’wv’，那么输出参数sorh一定为’s’，即使用软阈值法。

ddencmp函数中，阈值thr的计算是用下式计算的

式(3)

其中，1000为图19所示程序中输入信号x的数据长度，c为x经’db1’小波一次分解后得到的细节分量。下面需要研究下median函数，它应该不同于mean函数。

If A is a vector, median(A) returns the median value of A.

Median value是指什么？是指平均值还是从小到大排列后序列的中间值？

经分析median函数代码后得，对于一维向量A，如果该向量长度为奇数，那么median(A)为A经从小到大重新排列后的中间值；如果该向量长度为偶数，那么median(A)为A经从小到大重新排列后的中间两个值的平均值。如下表所示。

表4

A

Median(A)

[1 2 3 5 9]

3

[2 5 1 9 3]

3

[1 2 3 9]

2.5

[3 1 9 2]

2.5

10 thselect函数

选取用于小波消噪处理的阈值。

thr=thselect(x,tptr)：根据信号x和阈值选择标准tptr来确定一个消噪处理过程中所采用的自适应的阈值。

阈值选择标准tptr有下表所示的4种取值方式。

表5

tptr

取值方式说明

‘rigrsure’

使用stein的无偏似然估计原理所得到的自适应阈值

‘heursure’

启发式阈值选择，为最优预测变量阈值选择

‘sqtwolog’

固定阈值形式，大小为sqrt(2*log(length(x)))

‘minimaxi’

采用极大极小值原理选择阈值

阈值的选择规则是基于基本模型y=f(t)+e,其中e是白噪声N(0,1)。对于方差未知的噪声或非白噪声可以重新调整输出阈值。

以下图所示程序为例分析这几种阈值的计算方法。

clear;

load noissin;

x=noissin;

thr1 = thselect(x,rigrsure)

thr2 = thselect(x,sqtwolog)

thr3 = thselect(x,heursure)

thr4 = thselect(x,minimaxi)

图20

当tptr=rigrsure时，thr1的计算过程相对繁琐，用公式表达不如还是用编程语句表达，如下图所示。

sx2=sort(abs(x)).^2;

risks=(n-(2*(1:n))+(cumsum(sx2)+(n-1:-1:0).*sx2))/n;

[risk,best] = min(risks);

thr1 = sqrt(sx2(best));

图21

上图所示程序中，sort命令用于排序，这里还需要解释的命令有cumsum函数和min函数。

对cumsum函数，设向量x，y=cumsum(x),那么y可由下式表示。

式(4)

对min函数，设a为向量，[y,v] = min(a)，y为最小值，v为该最小值所在的位置。

当tptr=sqtwolog时，thr = sqrt(2*log(length(x)))。

当tptr=heursure时，thr3的计算过程可用下列编程语句表示。

hthr = sqrt(2*log(n));

eta = (norm(x).^2-n)/n;

crit = (log(n)/log(2))^(1.5)/sqrt(n);

if eta < crit

thr = hthr;

else

thr3 = min(thselect(x,rigrsure),hthr);

end

图22

上图所示程序中，norm为求模命令，设x为向量，y=norm(x),则有下式成立。

式(5)

当tptr=minimaxi时，thr4的计算过程可用下列编程语句表示

if n <= 32,thr = 0;

else thr4 = 0.3936 + 0.1829*(log(n)/log(2));end

图23

11 wden函数

调用方式[xd,cxd,lxd]=wden(x,tptr,sorh,scal,n,’wname’)或[xd,cxd,lxd]=wden(c,l,tptr,sorh,scal,n,’wname’)，返回经过小波消噪处理后的信号xd，及其小波分解结构[cxd,lxd]。

输入参数tptr为阈值选择标准，它的参数可以选择四种，’rigrsure’，’heursure’，’sqtwolog’和’minimaxi’。

输入参数sorh为函数选择的阈值使用方式，即sorh=’s’，为软阈值；sorh=’h’，为硬阈值。

输入参数scal规定了阈值求解随噪声水平的变化，其可选值为scal=’one’，不随噪声水平变化；scal=’sln’，根据第一层小波分解的噪声水平估计进行调整；scal=’mln’，根据每一层小波分解的噪声水平进行调整。

以下图所示程序为例仔细分析这个问题。

clear;

snr = 3; init = 2055615866;

[xref,x] = wnoise(3,11,snr,init);

lev = 5;

xd = wden(x,heursure,s,one,lev,sym8);

subplot(611), plot(xref), axis([1 2048 -10 10]);

title(Original signal);

subplot(612), plot(x), axis([1 2048 -10 10]);

title([Noisy signal - Signal to noise ratio = ,...

num2str(fix(snr))]);

subplot(613), plot(xd), axis([1 2048 -10 10]);

title(De-noised signal - heuristic SURE);

xd = wden(x,rigrsure,s,one,lev,sym8);

subplot(614), plot(xd), axis([1 2048 -10 10]);

title(De-noised signal - SURE);

xd = wden(x,sqtwolog,s,sln,lev,sym8);

subplot(615), plot(xd), axis([1 2048 -10 10]);

title(De-noised signal - Fixed form threshold);

xd = wden(x,minimaxi,s,sln,lev,sym8);

subplot(616), plot(xd), axis([1 2048 -10 10]);

title(De-noised signal - Minimax);

图24

上图所示程序中，wnoise函数用于产生信噪比为3的特定信号。

下面研究下程序中的wden函数，该函数由以下四个环节组成。

小波分解(wavedec)

求阈值(thselect)

软硬阈值处理，处理细节分量(wthresh)

)

小波重构(waverec)

图25

12 wdencmp函数

一维信号的小波消噪或压缩处理。在一维信号降噪方面，wdencmp函数与wden函数相比，少了求阈值环节，该函数的阈值需要由输入变量thr输入。wdencmp函数主要有三个主要环节。

小波分解(wavedec)

软硬阈值处理，处理细节分量(wthresh)

)

小波重构(waverec)

图26

在一维信号消噪方面，该函数的输入变量的不同组合只对上图所示三个环节中的第二个环节产生影响，如下表所示，共有三种不同的阈值处理方式。

表6

第一个输入变量

KEEPAPP

阈值处理方式

‘gbl’

KEEPAPP=1

只对信号的细节部分进行统一阈值处理

KEEPAPP≠1

对信号的细节部分和近似部分均进行统一阈值处理

‘lvd’

无此输入变量

对信号的各层细节部分用不同的阈值处理

其对应的示例程序如下图所示。

clear;

load leleccum; indx = 2600:3100;

x = leleccum(indx);

thr=35;

[xd,cxd,lxd] = wdencmp(gbl,x,db3,3,thr,h,0);

[xd1,cxd,lxd] = wdencmp(gbl,x,db3,3,thr,h,1);

thr=[35 32 30];

[xd2,cxd,lxd] = wdencmp(lvd,x,db3,3,thr,h);

plot(x);

hold on;

plot(xd,g);

plot(xd1,y.);

plot(xd2,r-.);

hold off;

图27

图28

13 wnoisest函数

估计小波分量中细节分量的标准差。由于函数wavedec的输出变量[c,l]的数据存放形式，不方便直接计算各个细节分量的标准差，故有此函数，就像提取细节分量函数detcoef一样。可想而知，该函数的语句逻辑是很简单的。

14 wthcoef函数

该函数相当于阈值处理函数，功能上比软硬阈值处理函数wthresh更加多些。下面以下图所示程序说明这个函数的用法。

clear;

load noissin;

x=noissin(1:1000);

[c,l]=wavedec(x,3,db4);

N=[1,2,3];

P=[98,99,97];

nc=wthcoef(d,c,l,N,P);

snc1=waverec(nc,l,db4);

nc=wthcoef(d,c,l,N);

snc2=waverec(nc,l,db4);

nc=wthcoef(a,c,l);

snc3=waverec(nc,l,db4);

T=[0.4 0.45 0.42];

nc=wthcoef(t,c,l,N,T,s);

snc4=waverec(nc,l,db4);

subplot(321);plot(x);

subplot(322);plot(snc1,r);

subplot(323);plot(snc2,g);

subplot(324);plot(snc3,y);

subplot(325);plot(snc4,k);

图29

如上图程序所示，该函数共有四种不同的用法。输入信号x经db4小波3层分解后得到CA3,CD3,CD2和CD1。下表对此作了总结。

表7

第一个输入变量

N

P

T

SORH

对应功能

‘d’

[1,2,3]

[98,99,97]

无此输入变量

无此输入变量

CD1按绝对值从小到大排列，前98%的部分置0；CD2，CD3依此类推

无此输入变量

CD3,CD2和CD1全部置0

‘a’

无此输入变量

无此输入变量

CA3全部置0

‘t’

[1,2,3]

无此输入变量

[0.4 0.45 0.42]

‘s’

CD1按阈值0.4作软阈值处理；CD2,CD3依此类推

15 wthresh函数

软硬阈值处理函数。下图所示程序和运行结果可以比较清晰地看出该程序的执行过程。

clear;

y = linspace(-1,1,100);

thr = 0.4;

ythard = wthresh(y,h,thr);

ytsoft = wthresh(y,s,thr);

subplot(311);plot(y);

subplot(312);plot(ythard);

subplot(313);plot(ytsoft);

图30

图31

16 wbmpen函数

计算阈值函数。

17 wdcbm函数

计算阈值函数。

19

展开阅读全文