小波分析的理解

但是小波变换的小波系数为如何选择小波基函数提供了依据。小波变换后的系数比较大，就表明了小波和信号的波形相似程度较大；反之则比较小。    另外还要根据信号处理的目的来决定尺度的大小。如果小波变换仅仅反映信号整体的近似特征，往往选用较 大的尺度；反映信号细节的变换则选用尺度不大的小波。由于小波函数家族成员较多，进行小波变换目的各异，目前没有一个通用的标准。

根据实际运用的经验， Morlet 小波应用领域较广，可以用于信号表示和分类、图像识别 特征提取；墨西哥草帽小波用于系统识别；样条小波用于材料探伤； Shannon 正交基用于差

分方程求解。

是不是小波以一个尺度分解一次就是小波进行一层的分解？

比如： [C,L]=wavedec(X,N,'wname') 中， N 为尺度，若为 1 ，就是进行单尺度分解，也就是分解一层。      但是 W=CWT(X,[2:2:128],'wname','plot') 的分解尺度又是从 2 ～ 128 以 2 为步进的，这里的“分解尺度”跟上面那个“尺度”的意思一样吗？

[C,L]=wavedec(X,N,'wname') 中的 N 为分解层数 ,  不是尺度 ,' 以 wname' 是 DB 小波为例 ,  如 DB4, 4 为消失矩 , 则一般滤波器长度为 8,  阶数为 7.

wavedec 针对于离散 ,CWT 是连续的。

多尺度又是怎么理解的呢？

多尺度的理解 :   如将 0-pi 定义为空间 V0,   经过一级分解之后 V0 被分成 0-pi/2 的低频子空间 V1 和 pi/2-pi 的高频子空间 W1,   然后一直分下去 .... 得到  VJ+WJ+....W2+W1.     因为 VJ 和 WJ 是正交的空间 ,  且各 W 子空间也是相互正交的 .   所以分解得到了是相互不包含的多个频域区间 , 这就是多分辩率分析 ,  即多尺度分析 .

.

是不是说分解到 W3 、 W2 、 W1 、 V3 就是三尺度分解？

简单的说尺度就是频率，不过是反比的关系．确定尺度关键还要考虑你要分析信号的采样频率大小，因为根据采样频率大小才能确定你的分析频率是多少．（采样定理）．然后再确定你到底分多少层．

假如我这有一个 10hz 和 50hz 的正弦混合信号，采样频率是 500hz ，是不是就可以推断出 10hz 和 50hz 各自对应的尺度了呢？我的意思是，是不是有一个频率和尺度的换算公式？

尺度

可以这样理解吗？如果我要得到频率为 0.125~0.25 的信号信息，是不是直接对 d3 的分解系数直接重构之后就是时域信息了？这样感觉把多层分解纯粹当作滤波器来用了，又怎么是多分辨分析？？怎样把时频信息同时表达出来？？

这个问题非常好，我刚开始的时候也是被这个问题困惑住了，咱们确实是把它当成了滤波器来用了，也就是说我们只看重了小波分析的频域局部化的特性。但是很多人都忽略其时域局部化特性，因为小波是变时频分析的方法，根据测不准原理如果带宽大，则时窗宽度就要小。那么也就意味着如果我们要利用其时域局部化特性就得在时宽小的分解层数下研究，也就是低尺度下。这样我们就可以更容易看出信号在该段时间内的细微变化，但是就产生一个问题，这一段的频率带很宽，频率局部化就体现不出来了。

对 d3 进行单支重构就可以得到 0.125 － 0.25 的信号了，当然频域信息可能保存的比较好，

但如果小波基不是对称的话，其相位信息会失真。

Daubechies 小波基的构造

%   此程序实现构造小波基

%  periodic_wavelet.m

function ss=periodic_wavelet;

clear;clc;

% global MOMENT;  %   消失矩阶数

% global LEFT_SCALET;  %   尺度函数左支撑区间

% global RIGHT_SCALET;  %   尺度函数右支撑区间

% global LEFT_BASIS;  %   小波基函数左支撑区间

% global RIGHT_BASIS;  %   小波基函数右支撑区间

% global MIN_STEP;  %   最小离散步长

% global LEVEL;  %   计算需要的层数 ( 离散精度）

% global MAX_LEVEL;  %  

周期小波最大计算层数

[s2,h]=scale_integer;

[test,h]=scalet_stretch(s2,h);

wave_base=wavelet(test,h);

ss=periodic_waveletbasis(wave_base);

function [s2,h]=scale_integer;

%   本函数实现求解小波尺度函数离散整数点的值

%  sacle_integer.m

MOMENT=10;  %   消失矩阶数

LEFT_SCALET=0;  %   尺度函数左支撑区间

RIGHT_SCALET=2*MOMENT-1;  %   尺度函数右支撑区间

LEFT_BASIS=1-MOMENT;    %   小波基函数左支撑区间

RIGHT_BASIS=MOMENT;     %   小波基函数右支撑区间

MIN_STEP=1/512;          %   最小离散步长

LEVEL=-log2(MIN_STEP);  %   计算需要的层数 ( 离散精度）

MAX_LEVEL=8;  %  

周期小波最大计算层数

h=wfilters('db10','r');  % 滤波器系数

h=h*sqrt(2); % FI(T)=SQRT(2)*SUM(H(N)*FI(2T-N)) N=0:2*MOMENT-1;

for i=LEFT_SCALET+1:RIGHT_SCALET-1

for j=LEFT_SCALET+1:RIGHT_SCALET-1

k=2*i-j+1;

if (k>=1&k<=RIGHT_SCALET+1)

a(i,j)=h(k);  %   矩阵系数矩阵

else

a(i,j)=0;

end

end

end

[s,w]=eig(a);  %   求特征向量，解的基

s1=s(:,1);

s2=[0;s1/sum(s1);0]; %   根据条件 SUM(FI(T))=1, 求解 ;

%   本函数实现尺度函数经伸缩后的离散值

%  scalet_stretch.m

function [s2,h]=scalet_stretch(s2,h);

MOMENT=10;  %   消失矩阶数

LEFT_SCALET=0;  %   尺度函数左支撑区间

RIGHT_SCALET=2*MOMENT-1;  %   尺度函数右支撑区间

LEFT_BASIS=1-MOMENT;  %   小波基函数左支撑区间

RIGHT_BASIS=MOMENT;  %   小波基函数右支撑区间

MIN_STEP=1/512;  %   最小离散步长

LEVEL=-log2(MIN_STEP);  %   计算需要的层数 ( 离散精度）

MAX_LEVEL=8;  %  

周期小波最大计算层数

for j=1:LEVEL  % 需要计算到尺度函数的层数

t=0;

for i=1:2:2*length(s2)-3  %   需要计算的离散点取值（ 0 ， 1 ， 2 ， 3 -> 1/2, 3/2, 5/2)

t=t+1;

fi(t)=0;

for n=LEFT_SCALET:RIGHT_SCALET;  %  低通滤波器冲击响应紧支撑判断

if ((i/2^(j-1)-n)>=LEFT_SCALET&(i/2^(j-1)-n)<=RIGHT_SCALET) %   小波尺度函数紧支撑判断

fi(t)=fi(t)+h(n+1)*s2(i-n*2^(j-1)+1);  %   反复应用双尺度方程求解

end

end

end

clear s

n1=length(s2);

n2=length(fi);

for i=1:length(s2)+length(fi)  %   变换后的矩阵长度

if (mod(i,2)==1)

s(i)=s2((i+1)/2);  %   矩阵奇数下标为小波上一层 (0,1,2,3) 离散值

else

s(i)=fi(i/2);  %   矩阵偶数下标为小波下一层 (1/2,3/2,5/2)( 经过伸缩变换后 ) 的离散值

end

end

s2=s;

end

% 采用双尺度方程求解小波基函数  PSI(T)

%  wavelet.m

function wave_base=wavelet(test,h);

MOMENT=10;  %   消失矩阶数

LEFT_SCALET=0;  %   尺度函数左支撑区间

RIGHT_SCALET=2*MOMENT-1;  %   尺度函数右支撑区间

LEFT_BASIS=1-MOMENT;  %   小波基函数左支撑区间

RIGHT_BASIS=MOMENT;  %   小波基函数右支撑区间

MIN_STEP=1/512;  %   最小离散步长

LEVEL=-log2(MIN_STEP);  %   计算需要的层数 ( 离散精度）

MAX_LEVEL=8;  %  

周期小波最大计算层数

i=0;

for t=LEFT_BASIS:MIN_STEP:RIGHT_BASIS;  %   小波基支撑长度

s=0;

for n=1-RIGHT_SCALET:1-LEFT_SCALET  %  g(n) 取值范围

if((2*t-n)>=LEFT_SCALET&(2*t-n)<=RIGHT_SCALET)  %   尺度函数判断

s=s+h(1-n+1)*(-1)^(n)*test((2*t-n)/MIN_STEP+1);  %   计算任意精度的小波基函数值  

end

end

i=i+1;

wave_base(i)=s;

end

一维数字滤波器 filter():

Y=filter(B, A, X)    由传递函数模型向量 B 、 A 描述的滤波器对向量 X 中的元素进行滤波，并将结果数据存放在向量 Y 中。

[Y, Zf]=filter(B, A, X, Zi)    给出了滤波器延时的初始和终止条件 Zf 和 Zi 。

例子：

         人体心电信号在测量过程中往往受到工业高频干扰，所以必须经过低通滤波处理后，才能判断心脏功能的有用信息。下面给出一实际心电图信号采样序列样本 x(n) ，其中存在高频干扰。在试验中以 x(n) 作为输入序列，滤除其中的干扰成分 。

{ x(n) } = {-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,

-2,-4,8,12,12,10,6,6,4,0,0,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0}

Matlab程序设计如下：

X=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,4,0,0,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0];

figure;

plot(X);

xlabel(' 时间 ');

ylabel(' 幅值 ');

wp=40;   ws=50;  rp=0.5;  rs=40;   Fs=200;

[N, Wn] = buttord(wp/(Fs/2), ws/(Fs/2), rp, rs);

[b, a]=butter(N, Wn);

figure;

[H, W]=freqz(b, a);

plot(W*Fs/(2*pi), abs(H));   grid;

xlabel(' 频率 /Hz');

ylabel(' 幅值 ');

Y=filter(b,a,X);

figure

plot(Y)

xlabel(' 时间 ');

ylabel(' 幅值 ');

figure

psd(X, [ ],200);

figure

psd(Y, [ ],200);

分析这段程序可知包括以下几部分：

（ 1 ）首先绘制原始数据的图形；

（ 2 ）设计一个 Butterworth 低通滤波器并绘制出它的幅频响应曲线；

（ 3 ）用设计的滤波器对原数据进行滤波；

（ 4 ）绘制滤波以后的数据图形；

（ 5 ）绘制原数据功率谱图形；

）绘制滤波以后数据功率谱图形。

1. 监测信号的自相似性

直观上讲，小波分解系数表示了信号与小波之间的“相似指数”，如果相似程度越高，则相似指数越大。因此如果一个信号的不同的尺度之间相似，则小波系数在不同的尺度上也应该相似。因此可以通过小波分解检测信号的自相似性，即检测信号的分形特征。实践表明，通过小波分解可以很好地研究信号或图像的分形特征。

下面通过一个简单的例子来说明小波分析在检测信号自相似性中的应用，待检测的信号是经过反复迭代生成的信号，因此具有自相似性。

程序代码如下：

load vonkoch;

x=vonkoch;

subplot(211);

plot(x);

title(' 原始信号 ');

subplot(212);

% 进行一维连续小波变换

f=cwt(x,[2:2:128],'coif3','plot');

从图中可以看出分解后的小波系数在许多尺度上看上去都非常相似。

2. 信号的奇异性检测

信号的突变点和奇异点等不规则部分通常包含重要信息。

一般信号的奇异性分为两种情况：（ 1 ）信号在某一时刻其幅值发生突变，引起信号的非连续，这种类型的突变称为第一类型的间断点；     （ 2 ）信号在外观上很光滑，幅值没有发生突变，但是信号的一阶微分有突变发生且一阶微分不连续，这种类型的突变称为第二类型的间断点。

应用小波分析可以检测出信号中的突变点的位置、类型以及变化的幅度。下面介绍小波分析在信号奇异性检测中的应用。

（ 1 ）第一类型间断点的检测

下面举例说明小波分析用于检测第一类型的间断点。

在本例中，信号的不连续是由于低频特征的正弦信号在后半部分突然有高频特征的正弦信号加入，首先利用傅里叶变换分析对信号在频域进行分析，发现无检测突变点，接着利用小波分析进行分析，结果证明它能够准确地检测出了信号幅值突变的位置，即高频信号加入的时间点。

程序代码如下：

load  freqbrk;

x=freqbrk;

% 对信号进行傅里叶变换

f=fft(x,1024);

f=abs(f);

figure;

subplot(211);

plot(x);

subplot(212);

plot(f);

% 使用 db6 小波进行 6 层分解

[c,l]=wavedec(x,6,'db6');

figure(2);

subplot(811);

plot(x);

ylabel('x');

% 对分解的第六层低频系数进行重构

a=wrcoef('a',c,l,'db6',6);

subplot(812);

plot(a);

ylabel('a6');

for i=1:6

% 对分解的第 6 层到第 1 层的高频系数分别进行重构

d=wrcoef('d',c,l,'db6',7-i);

subplot(8,1,i+2);

plot(d);

ylabel(['d',num2str(7-i)]);

end

由图中可以看出，由于傅里叶变换不具有时间分辨力，因此无法检测信号的间断点。而在小波分析的图中，在信号的小波分解的第一层高频系数 d1 和第二层高频系数 d2 中，可以非常清楚地观察到信号的不连续点，用 db1 小波比用 db6 小波要好。

这个例子也表明小波分析在检测信号的奇异点时具有傅里叶变换无法比拟的优越性，利用小波分析可以精确地检测出信号的突变点。

程序代码如下：

load sumsin;

x=sumsin;

figure;

subplot(611);

plot(x);

ylabel('x');

title(' 原始信号以及各层近似信号 ');

% 使用 db3 小波进行 5 层分解

[c,l]=wavedec(x,5,'db3');

for i=1:5

% 对分解的第 5 层到第 1 层的低频系数分别进行重构

a=wrcoef('a',c,l,'db3',6-i);

subplot(6,1,i+1);

plot(a);

ylabel(['a',num2str(6-i)]);

end

figure;

subplot(611)

plot(x);

ylabel('x')

for i=1:5

% 对分解的第 5 层到第 1 层的高频系数进行重构

d=wrcoef('d',c,l,'db3',6-i);

subplot(6,1,i+1);

plot(d);

ylabel(['d',num2str(6-i)]);

end

分析：

。