《小波分析基础》

小波分析的基本知识屠2001.8.2A.基本知识A1.小波(WAVELET)分类1.原始小波:(1).高斯gaus, (2).莫来特morlet, (3).墨西哥帽mexihat2.无限正则小波浪:(4).梅耶meyr (5).离散梅耶dmey3.正交和紧支集小波:(6).达比切斯dbN(Daubechies), (7).对称symN(symlets), (8).coifN4.双正交和紧支集小波: (9).双正交biorNr (10). 逆双正交rbioNr.Nd5.复小波: (11).复高斯cgauN (12)复莫来特 cmor Fb-Fc (13)复香农shan Fb-Fc(14).复频率B样条 fbspM Fb-Fc注1:db1小波也称哈尔Harr小波,也是原始小波注2.symlet小波是Daubechies小波的改进,由不对称改成近似对称注3.紧支集即函数在有限区域内不为零A2.小波函数和尺度(SCALE)函数1.小波函数(psi)--由高通滤波器确定,产生小波分解的细节 D(detail,)2.尺度函数(phi)--由低通正交镜象滤波器确定, 产生小波分解的逼近 A(approximation)A3.小波分解:S(SIGNAL)=A1+D1=(A2+D2)+D1=(A3+D3)+D2+D1=(A4+D4)+D3+D2+D1=...A4.小波包(WP=Wavelet Packet)分解:S=A1+D1=(AA2+DA2)+(AD2+DD2)=(AAA3+DAA3)+(ADA3+DDA3)+(ADA3+DDA3)+...A5.WAVEMENU: 开始图象用户界面GUI工具A6.WAVEDEMO: 小波工具箱演示B小波变换B1.一维连续小波变换:CWT coefs=cwt(S,scales,"wname')coefs=cwt(S,scales,'wname','plot')coefs=cwt(S,scales,'wname','plotmode')scales--正实数,如1:32,[64 32 16:-2:2],...COLORMENU,COLORBARB2.单级一维离散小波变换:DWT,UPCOEF[Ca1,Cd1]=dwt(x,'wname'), Ca1--逼近系数 Cd1--细节系数[Ca1,Cd1]=dwt(x,Lo_D,Hi_D)a1=upcoef('a',Ca1,'wname',1,L); a1--逼近 L--length(x)d1=upcoef('d',Cd1,'wname',1,L); d1--细节 L--length(x)B3.单级一维逆离散小波变换:IDWT, x=idwt(Ca1,Cd1,'wname')B4.多级一维离散小波分解:WAVEDEC,APPCOEF,DETCOEF,WRCOEF[C,L]=wavedec(x,N,'wname'),N--级(LEVEL)数 C--分解(DECOMPOSITION)矢量L--辅助操作(Bookkeeping)矢量B5.APPCOEF:提取一维小波逼近系数,A=appcoef(C,L,'wname',N)B6.DETCOEF:提取一维小波细节系数,A=detcoef(C,L,'wname',N)B7.WRCOEF:X=wrcoef('type',C,L,'wname',N).type=a,逼近;type=d,细节B8.WAVEREC(多级一维离散小波重构) 重构--RECONSTRUCTIONx=waverec(C,L,'wname')x=waverec(C,L,Lo_R,Hi_R)B9.WFILTERS--小波滤波器[Lo_D,Hi_D,Lo_R,Hi_R]=wfilters('wname'),'wname'=db,coif,sym,bior,rbioB10.DYADDOWN:二进(Dyadic)降采样 Y=dyaddown(x,EVENODD)EVENODD--even,y(k)=x(2k), --odd,y(k)=x(2k-1)B11. DYADUP:二进增采样(填零), y=dyadup(x,EVENODD)B12. WKEEP:保留矢量或矩阵的一部分C.小波包变换C1. WPDEC一维离散小波包分解:[T,D]=wpdec(x,N,'wname',E,P), T--树结构Tree structure, D--数据结构E-熵 Entropy E='shannon','threshold','norm','log energy','user'P-附加参数'threshold' 'sure':P=threshold(0<=P)'norm':P=power,1<=P<2)C2. WPREC一维离散小波包重构x=wprec(T,D) T--小波包树(TREE) N—节点(NODE)C3. WPCOEF小波包系数x=wpcoef(S,D,N)D.MALLAT算法(FWT)E.一维试验信号(b1(t): b2(t): )oislop(ramp+color noise):1=<t<=499,(t/500)+b2(t);500=<t<=1000,1+b2(t)2.freqbrk:1=<t<=500,sin(0.03t);501=<t<=1000,sin(0.3t)3. heavysin4.nelec(2000 电力消耗)5.leleccum(4320分(72小时)电力消耗6.linchirp(线性快扫)7.mfreqbrk8.mishmash 9.nearbrk(1~499,511~1500)10.noisbloc 11.noisbump12.noischir 13.noisdopp 噪声多普勒14.noismima 15.noispol: 在[1 1000]间 t^2-t+1+b1(t) 16.noissin:sin(0.03t)+b1(t) 17.qdchirp18.quachip19.scddvbrk:二阶导数不连续,t<0,exp(-4t^2);t>=0,exp(-t^2),t=[-0.5 0.5]20.sinfract 21.sinper8 22.sumlichr23.sumsin:sin(3t)+sin(0.3t)+sin(0.03t)24.trsin:1=<t<=500,((t-1)/500)+sin(0.3t);501=<t<=1000,((1000-t)/500)+sin(0.3t)25.vonkoch:分形,科克雪花26.warma:AR(3),b2(t)=-1.5*b2(t-1)-0.75*b2(t-2)-0.125*b2(t-3)+b1(t)+0.527.wcantor:分形,康托(三分取一)曲线28.whitnois:在[-0.5 0.5]间的均匀白噪声29.wnoislop:1=<t<=499,(3t/500)+b1(t);500=<t<=1000,3+b1(t)30.wntrsin:1=<t<=500,((t-1)/500)+sin(0.3t)+b1(t);501=<t<=1000,((1000-t)/500)+sin(0.3t)+b1(t)31.wstep:1=<t<=500,s=0;501=<t<=1000,s=20.32.cuspamax(1024):x=linspace(0,1.1024);y=exp(-128*((X-0.3).^2))-3*(abs(x-0.7).^0.433.brkintri:顶端折线三角34.wcantsym(2188):对称康托集disp('******)*************MALLAT算法示例***********************************************')x=[1.8 1.0 -1.0 -1.8];[Lo_D,Hi_D]=wfilters('db1','d');tmpo1=conv(x,Lo_D); [1.8 1.0 -1.0 -1.8]*[0.7071 0.7071]tmpo2=conv(x,Hi_D);Ca1=dyaddown(tmpo1);Cd1=dyaddown(tmpo2);disp('低通分解滤波器系数Lo_D 高通分解滤波器系数Hi_D');disp( [(Lo_D)' (Hi_D)'] ),disp('卷积conv(x,Lo_D 卷积conv(x,Hi_D)');disp( [(tmpo1)’ (tmpo2)’] ),disp('一级逼近系数Ca1 一级细节系数Cd1');disp( [(Ca1)’ (Cd1)’] ),% Ca1=1.9799 -1.9799 Cd1= 0.5657 0.5657[Lo_R,Hi_R]=wfilters('db1','r');disp('低通重构滤波器系数Lo_R=');disp(Lo_R),disp('高通重构滤波器系数Hi_R=');disp(Hi_R),tmp1=dyadup(Cd1);tmpo3=conv(tmp1,Hi_R);d1=wkeep(tmpo3,4);tmp2=dyadup(Ca1);tmpo4=conv(tmp2,Lo_R);a1=wkeep(tmpo4,4);disp( '一级逼近a1 一级细节d1 ');DISP( [(a1)’ (d1)’] ),% 一级逼近a1= 1.4000 1.4000 -1.4000 -1.4000% 一级细节d1= 0.4000 -0.4000 0.4000 -0.4000figure(1),a0=a1+d1;subplot(521),bar(x,0.1),title('原始波形x=[1.8 1.0 -1.0 -1.8]'), grid,axis([0 5 -2 2]),subplot(522),bar(a0,0.1),title('分解后重构波形s=a1+d1'),grid,axis([0 5 -2 2])subplot(523),bar(Ca1,0.1),title(' 逼近系数Ca1=[1.98 -1.98]'),grid,axis([0 5 -2 2])subplot(524),bar(Cd1,0.1),title(' 细节系数Cd1=[0.566 0.566]'),grid,axis([0 5 0 1])subplot(525),bar(a1,0.1),title(' 一级逼近a1=[1.4 1.4 -1.4 -1.4]'),grid,axis([0 5 -2 2])subplot(526),bar(d1,0.1),title(' 一级细节d1=[0.4 -0.4 0.4 -0.4]'),grid,axis([0 5 -1 1])subplot(527),bar(Lo_D,0.1),title('低通分解滤波器系数Lo_D'),grid,axis([0 5 0 1])subplot(528),bar(Hi_D,0.1),title('高通分解滤波器系数Hi_D'),grid,axis([0 5 -1 1])subplot(5,2,9),bar(Lo_R,0.1), title('低通重构滤波器系数Lo_R'),grid,axis([0 5 0 1])subplot(5,2,10),bar(Hi_R,0.1),title('高通重构滤波器系数Hi_R'),grid,axis([0 5 -1 1])%******以上为MALLAT算法原理,实际上用简单命令DWT,UPCOEF计算如下************************** x=[1.8 1.0 -1.0 -1.8];length(x);[Ca1,Cd1]=dwt(x,'db1');a1=upcoef('a',Ca1,'db1',1,4);d1=upcoef('d',Cd1,'db1',1,4);x1=a1+d1;a0=idwt(Ca1,Cd1,'db1',4);------------------------------------------------------------------------------------ x=[1.8 -1.8 1.8 -1.8];x=[1.8 1.0 -1.0 -1.8];[Lo_D,Hi_D]=wfilters('db40','d');tmpo1=conv(x,Lo_D);tmpo2=conv(x,Hi_D);Ca1=dyaddown(tmpo1);Cd1=dyaddown(tmpo2);disp('低通分解滤波器系数Lo_D 高通分解滤波器系数Hi_D');disp( [(Lo_D)'(Hi_D)'] )disp('卷积conv(x,Lo_D)');disp(tmpo1),disp('卷积conv(x,Hi_D)');disp(tmpo2),disp('一级逼近系数Ca1=');disp(Ca1),disp('一级细节系数Cd1=');disp(Cd1),[Lo_R,Hi_R]=wfilters('db40','r');disp('低通重构滤波器系数Lo_R=');disp(Lo_R),disp('高通重构滤波器系数Hi_R=');disp(Hi_R),tmp1=dyadup(Cd1);tmpo3=conv(tmp1,Hi_R);d1=wkeep(tmpo3,4);tmp2=dyadup(Ca1);tmpo4=conv(tmp2,Lo_R);a1=wkeep(tmpo4,4);disp('一级逼近a1');disp(a1),disp('一级细节d1');disp(d1),figure(2),a0=a1+d1;subplot(521),bar(x,0.1),title('原始波形x=[1.8 1.0 -1.0 -1.8]'),subplot(521),bar(x,0.1),title('原始波形x=[1.8 -1.8 1.8 -1.8]'),grid,axis([0 5 -2 2]),subplot(522),bar(a0,0.1),title('分解后重构波形s=a1+d1'),grid,axis([0 5 -2 2]) subplot(523),bar(Ca1,0.1),title(' 逼近系数Ca1'),grid,axlimdlg, axis([0 5 -2 2]) subplot(524),bar(Cd1,0.1),title(' 细节系数Cd1'),grid,axlimdlg, axis([0 5 0 1])subplot(525),bar(a1,0.1),title(' 一级逼近a1=[1.296,0.911,-0.6502,-1.5585]'),subplot(525),bar(a1,0.1),title(' 一级逼近a1=[ ]'), axlimdlg,grid, axis([0 5 -2 2])subplot(526),bar(d1,0.1),title(' 一级细节d1=[0.504,0.089,-0.3498,-0.2415'), subplot(526),bar(d1,0.1),title(' 一级细节d1=[ ]'),axlimdlg,grid, axis([0 5 -1 1])subplot(527),bar(Lo_D,0.1),title('低通分解滤波器系数Lo_D'),grid, axis([0 5 0 1]) subplot(528),bar(Hi_D,0.1),title('高通分解滤波器系数Hi_D'),grid, axis([0 5 -1 1]) subplot(5,2,9),bar(Lo_R,0.1), title('低通重构滤波器系数Lo_R'),grid, axis([0 5 0 1]), axlimdlg,subplot(5,2,10),bar(Hi_R,0.1),title('高通重构滤波器系数Hi_R'),grid, axis([0 5 -1 1]) axlimdlg,k=[1.8 1.0 -1.0 -1.8];s=[1.296 0.911 -0.6502 -1.5585];t=[0.504 0.089 -0.3498 -0.2415];ss=abs(fft(s,21));tt=abs(fft(t,21));kk=abs(fft(k,21));subplot(311),plot(kk),grid,axlimdlg,subplot(312),plot(ss),grid,axlimdlg,subplot(313),plot(tt),grid,axlimdlg,k1=[1.8 1.0 -1.0 -1.8];s1=[1.4 1.4 -1.4 -1.4];t1=[0.4 -0.4 0.4 -0.4];k2=[1.8 1.0 -1.0 -1.8];s2=[1.296 0.911 -0.6502 -1.5585];t2=[0.504 0.089 -0.3498 -0.2415];S1=abs(fft(s1,21));T1=abs(fft(t1,21));K1=abs(fft(k1,21));S2=abs(fft(s2,21));T2=abs(fft(t2,21));K2=abs(fft(k2,21));subplot(321),plot(K1),grid,axis([1 11 0 6]),title('Harr')subplot(323),plot(S1),grid,axis([1 11 0 6]),title('Harr')subplot(325),plot(T1),grid,axis([1 11 0 2]),title('Harr')subplot(322),plot(K2),grid,axis([1 11 0 6]),title('db40')subplot(324),plot(S2),grid,axis([1 11 0 6]),title('db40')subplot(326),plot(T2),grid,axis([1 11 0 2]),title('db40')disp('**********MALLAT算法可用简单命令DWT,UPCOEF重算如下*******************')x=[1.8 1.0 -1.0 -1.8];length(x); =4[Ca1,Cd1]=dwt(x,'db1');a1=upcoef('a',Ca1,'db1',1,4);d1=upcoef('d',Cd1,'db1',1,4);disp('一级逼近系数Ca1=');disp(Ca1), disp('一级细节系数Cd1=');disp(Cd1),disp('一级逼近a1=');disp(a1), disp('一级细节d1=');disp(d1),x1=a1+d1;a0=idwt(Ca1,Cd1,'db1',4);figure(1),subplot(321),bar(x,0.1),title('x=a1+d1=[1.8 1.0 -1.0 -1.8]'),grid,axis([0 5 -2 2]),subplot(322),bar(a0,0.1),title('a0=idwt=x'),grid,axis([0 5 -2 2])subplot(323),bar(Ca1,0.1),title(' 逼近系数Ca1=[1.98 -1.98]'),grid,axis([0 5 -2 2])subplot(324),bar(Cd1,0.1),title(' 细节系数Cd1=[0.566 0.566]'),grid,axis([0 5 0 1]) subplot(325),bar(a1,0.1),title(' 一级逼近a1=[1.4 1.4 -1.4 -1.4]'),grid,axis([0 5 -1.5 1.5])subplot(326),bar(d1,0.1),title(' 一级细节d1=[0.4 -0.4 0.4 -0.4]'),grid,axis([0 5 -1 1])pausedisp(' *******************************************************************'), disp(' * *'), disp(' * *'), disp(' * 低通滤波器减低通滤波器等于带通滤波器 *'), disp(' * *'), disp(' *******************************************************************'), pause,f=-10:0.01:10;t=-50:1/20:50;y1=cos(2*pi*100*f);y2=cos(2*pi*100*t);y1(1:50)=zeros(1,50);y1(1952:2001)=zeros(1,50);y2(1:250)=zeros(1,250);y2(1752:2001)=zeros(1,250);yy=y1-y2;u=cos(2*pi*7*t);v=sinc(t);r=u.*v;U=fft(u);V=fft(v);R=fft(r);x1=real(ifft(y1));x2=real(ifft(y2));xx=real(ifft(yy));figure(2),subplot(331),plot(f,y1),axis([-12,12,0,1.1]),...title('低通(尺度) Y1(f),Fc=9.5Hz.'),subplot(332),plot(f,y2),axis([-12,12,0,1.1]),...title('低通(尺度) Y2(f),Fc=7.5Hz.'),subplot(333),plot(f,yy),axis([-12,12,0,1.1]),...title('带通(小波) YY(f),BW=2Hz.'),subplot(334),plot(t,ifftshift(x1)),axis([-5 5 -0.1 1.0]),...title('X1(t)=IFFT(Y1)'),xlabel('t(s)'),...subplot(335),plot(t,ifftshift(x2)),axis([-5 5 -0.3 0.8]),...title('X2(t)=IFFT(Y2)'),xlabel('t(s)'),subplot(336),plot(t,ifftshift(xx)),axis([-5 5 -0.1 0.16]),...title('X3(t)=IFFT(YY)'),xlabel('t(s)'),pausedisp(' ******************************************************************'), disp(' * *'), disp(' * 调制引起频移,低通变成带通 *'), disp(' * *'), disp(' ******************************************************************'), pause, figure(3),subplot(331),plot(t,u),axis([-2 2 -1.1 1.1]),title('u=cos(2pi*7t),t=-50~50'),subplot(332),plot(t,v),axis([-4 4 -0.3 1.1]),title('v=sinc(t),t=-50~50')subplot(333),plot(t,r),axis([-4 4 -0.9 1.1]),title('r=uv,t=-50~50')subplot(334),plot(f,abs(U)),axis([-5 5 0 900]),title('FFT(u),F=3Hz'),xlabel('Hz') subplot(335),plot(f,fftshift(abs(V))),axis([-5 5 0 23]),...title('V=FFT(v)),低通:Fc=0.5Hz'),xlabel('Hz')subplot(336),plot(f,abs(R)),axis([-5 5 0 11]),title('FFT(r),带通:BW=1Hz'),...xlabel('Hz'),pause,********************************************************************************* t1=-10:0.02:10;f1=0:0.05:50;ta=-20:0.02:20;tb=0:0.02:40;f1=0:1/40:50;x1=cos(2*pi*50*t1);x2=cos(2*pi*50*[0:0.02:20]);x3=cos(2*pi*50*[20:0.02:40]);xa=[zeros(1,500) x1 zeros(1,500)];xb=[x2 zeros(1,1000)];xc=[zeros(1,1000) x3]; fxa=fft(xa);fxb=fft(xb);fxc=fft(xc);subplot(531),plot(ta,xa),axis([-20 20 0 1.1]),title(''),subplot(532),plot(tb,xb),axis([0 40 0 1.1]),title(''),subplot(533),plot(tb,xc),axis([0 40 0 1.1]),title(''),subplot(534),plot(f1,abs(fxa)),grid,axlimdlg,title(''),subplot(535),plot(f1,abs(fxb)),grid,axlimdlg,title(''),subplot(536),plot(f1,abs(fxc)),grid,axlimdlg,title(''),subplot(537),plot(f1,angle(fxa)),grid,axlimdlg,title(''),subplot(538),plot(f1,angle(fxb)),grid,axlimdlg,title(''),subplot(539),plot(f1,angle(fxc)),grid,axlimdlg,title(''),subplot(5,3,10),plot(f1,unwrap(angle(fxa))),grid,axlimdlg,title(''),subplot(5,3,11),plot(f1,unwrap(angle(fxb))),grid,axlimdlg,title(''),subplot(5,3,12),plot(f1,unwrap(angle(fxc))),grid,axlimdlg,title(''),disp('**************************END***********************************'),。