小波分析的理解

小波变换是克服其他信号处理技术缺陷的一种分析信号的方法。小波由一族小波基函数构成，它可以描述信号时间（空间）和频率（尺度）域的局部特性。采用小波分析最大优点是可对信号进行实施局部分析，可在任意的时间或空间域中分析信号。小波分析具有发现其他信号分析方法所不能识别的、隐藏于数据之中的表现结构特性的信息，而这些特性对机械故障和材料的损伤等识别是尤为重要的。如何选择小波基函数目前还没有一个理论标准，常用的小波函数有Haar、Daubechies(dbN)、Morlet、Meryer、Symlet、Coiflet、Biorthogonal小波等15种。但是小波变换的小波系数为如何选择小波基函数提供了依据。小波变换后的系数比较大，就表明了小波和信号的波形相似程度较大；反之则比较小。另外还要根据信号处理的目的来决定尺度的大小。如果小波变换仅仅反映信号整体的近似特征，往往选用较大的尺度；反映信号细节的变换则选用尺度不大的小波。由于小波函数家族成员较多，进行小波变换目的各异，目前没有一个通用的标准。根据实际运用的经验，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各自对应的尺度了呢？我的意思是，是不是有一个频率和尺度的换算公式？实际频率＝小波中心频率×采样频率/尺度在小波分解中，若将信号中的最高频率成分看作是1，则各层小波小波分解便是带通或低通滤波器，且各层所占的具体频带为（三层分解）a1:0~0.5d1:0.5~1;a2:0~0.25d2:0.25~0.5;a3:0~0.125;d3:0.125~0.25可以这样理解吗？如果我要得到频率为0.125~0.25的信号信息，是不是直接对d3的分解系数直接重构之后就是时域信息了？这样感觉把多层分解纯粹当作滤波器来用了，又怎么是多分辨分析？？怎样把时频信息同时表达出来？？这个问题非常好，我刚开始的时候也是被这个问题困惑住了，咱们确实是把它当成了滤波器来用了，也就是说我们只看重了小波分析的频域局部化的特性。但是很多人都忽略其时域局部化特性，因为小波是变时频分析的方法，根据测不准原理如果带宽大，则时窗宽度就要小。那么也就意味着如果我们要利用其时域局部化特性就得在时宽小的分解层数下研究，也就是低尺度下。这样我们就可以更容易看出信号在该段时间内的细微变化，但是就产生一个问题，这一段的频率带很宽，频率局部化就体现不出来了。对d3进行单支重构就可以得到0.125－0.25的信号了，当然频域信息可能保存的比较好，但如果小波基不是对称的话，其相位信息会失真。小波变换主要也是用在高频特征提取上。层数不是尺度，小波包分解中，N应该是层数，个人理解对应尺度应该是2^N小波分解的尺度为a，分解层次为j。如果是连续小波分解尺度即为a。离散小波分解尺度严格意义上来说为a＝2^j,在很多书上就直接将j称为尺度，因为一个j就对应者一个尺度a。其实两者是统一的。小波基：一般从线性相位，消失矩，相似性，紧支撑等来选择。Daubechies小波基的构造%此程序实现构造小波基%periodic_wavelet.mfunctionss=periodic_wavelet;clear;clc;%globalMOMENT;%消失矩阶数%globalLEFT_SCALET;%尺度函数左支撑区间%globalRIGHT_SCALET;%尺度函数右支撑区间%globalLEFT_BASIS;%小波基函数左支撑区间%globalRIGHT_BASIS;%小波基函数右支撑区间%globalMIN_STEP;%最小离散步长%globalLEVEL;%计算需要的层数(离散精度）%globalMAX_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.mMOMENT=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;fori=LEFT_SCALET+1:RIGHT_SCALET-1forj=LEFT_SCALET+1:RIGHT_SCALET-1k=2*i-j+1;if(k=1&k=RIGHT_SCALET+1)a(i,j)=h(k);%矩阵系数矩阵elsea(i,j)=0;endendend[s,w]=eig(a);%求特征向量，解的基s1=s(:,1);s2=[0;s1/sum(s1);0];%根据条件SUM(FI(T))=1,求解;%本函数实现尺度函数经伸缩后的离散值%scalet_stretch.mfunction[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;%周期小波最大计算层数forj=1:LEVEL%需要计算到尺度函数的层数t=0;fori=1:2:2*length(s2)-3%需要计算的离散点取值（0，1，2，3-1/2,3/2,5/2)t=t+1;fi(t)=0;forn=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);%反复应用双尺度方程求解endendendclearsn1=length(s2);n2=length(fi);fori=1:length(s2)+length(fi)%变换后的矩阵长度if(mod(i,2)==1)s(i)=s2((i+1)/2);%矩阵奇数下标为小波上一层(0,1,2,3)离散值elses(i)=fi(i/2);%矩阵偶数下标为小波下一层(1/2,3/2,5/2)(经过伸缩变换后)的离散值endends2=s;end%采用双尺度方程求解小波基函数PSI(T)%wavelet.mfunctionwave_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;fort=LEFT_BASIS:MIN_STEP:RIGHT_BASIS;%小波基支撑长度s=0;forn=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);%计算任意精度的小波基函数值endendi=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(