关于小波分析的matlab程序

1、集关于小波分析的 matlab程序)小波滤波器构造和消噪程序 1小波谱分析 mallat算法经典程序 7小波包变换分析信号的 MATLAB程序 9利用小波变换实现对电能质量检测的算法实现 15基于小波变换的图象去噪Normalshrink算法17小波滤波器构造和消噪程序1.重构% mallet_wavelet.m%此函数用于研究Mallet算法及滤波器设计%此函数仅用于消噪a=pi/8; %角度赋初值b=pi/8;%低通重构FIR滤波器h0(n)冲激响应赋值hO=cos(a)*cos(b);h1=s in( a)*cos(b);h2=-si n(a)*si n(b);h3=cos(a)*si

2、n(b);low_co nstruct=h0,h1,h2,h3;L_fre=4; %滤波器长度low_decompose=low_construct(end:-1:1); % 确定 h0(-n)，低通分解滤波器for i_high=1:L_fre; % 确定 h1(n)=(-1)An,高通重建滤波器if(mod(i_high,2)=0);coefficie nt=-1;else coefficient=1;endhigh_construct(1,i_high)=low_decompose(1,i_high)*coefficient; endh1(-n)high_decompose=high_c

3、onstruct(end:-1:1); % 高通分解滤波器L_signal=100; % 信号长度n=1:L_signal; % 信号赋值f=10;t=0.001;y=10*cos(2*pi*50*n*t).*exp(-20*n*t);figure(1);plot(y);title('原信号');check1=sum(high_decompose); %h0(n) 性质校验check2=sum(low_decompose);check3=norm(high_decompose);check4=norm(low_decompose);l_fre=conv(y,low_decomp

4、ose); % 卷积 l_fre_down=dyaddown(l_fre); % 抽取 ,得低频细节 h_fre=conv(y,high_decompose);h_fre_down=dyaddown(h_fre); % 信号高频细节figure(2);subplot(2,1,1)plot(l_fre_down);title(' 小波分解的低频系数 ');subplot(2,1,2);plot(h_fre_down);title(' 小波分解的高频系数 ');l_fre_pull=dyadup(l_fre_down); %0 差值 h_fre_pull=dyadu

5、p(h_fre_down);l_fre_denoise=conv(low_construct,l_fre_pull); h_fre_denoise=conv(high_construct,h_fre_pull); l_fre_keep=wkeep(l_fre_denoise,L_signal); % 取结果的中心部分 ,消除卷积影响 h_fre_keep=wkeep(h_fre_denoise,L_signal);sig_denoise=l_fre_keep+h_fre_keep; % 信号重构 compare=sig_denoise-y; % 与原信号比较 figure(3);subplot

6、(3,1,1)plot(y);ylabel('y'); % 原信号subplot(3,1,2);plot(sig_denoise);ylabel('sig_denoise'); % 重构信号subplot(3,1,3);plot(compare);ylabel('compare'); % 原信号与消噪后信号的比较2 消噪、% 此函数用于研究 Mallet 算法及滤波器设计% 此函数用于消噪处理%角度赋值%此处赋值使滤波器系数恰为db9%分解的高频系数采用db9 较好 ,即它的消失矩较大%分解的有用信号小波高频系数基本趋于零%对于噪声信号高频分解系

7、数很大，便于阈值消噪处理l,h=wfilters('db10','d');low_construct=l;L_fre=20; % 滤波器长度low_decompose=low_construct(end:-1:1); % 确定 h0(-n) ，低通分解滤波器for i_high=1:L_fre; % 确定 h1(n)=(-1)An,高通重建滤波器if(mod(i_high,2)=0);coefficient=-1;elsecoefficient=1;endhigh_construct(1,i_high)=low_decompose(1,i_high)*coeff

8、icient;endhigh_decompose=high_construct(end:-1:1); % 高通分解滤波器 h1(-n)L_signal=100; % 信号长度n=1:L_signal; % 原始信号赋值f=10;t=0.001;y=10*cos(2*pi*50*n*t).*exp(-30*n*t);zero1=zeros(1,60); % 信号加噪声信号产生zero2=zeros(1,30);noise=zero1,3*(randn(1,10)-0.5),zero2;y_noise=y+noise;figure(1);subplot(2,1,1);plot(y);title(&

9、#39;原信号');subplot(2,1,2);plot(y_noise);title(' 受噪声污染的信号 ');check1=sum(high_decompose); %h0(n), 性质校验check2=sum(low_decompose);check3=norm(high_decompose);check4=norm(low_decompose);l_fre=conv(y_noise,low_decompose); % 卷积l_fre_down=dyaddown(l_fre); % 抽取 ,得低频细节 h_fre=conv(y_noise,high_decom

10、pose);h_fre_down=dyaddown(h_fre); % 信号高频细节figure(2);subplot(2,1,1)plot(l_fre_down);title(' 小波分解的低频系数 ');subplot(2,1,2);plot(h_fre_down);title(' 小波分解的高频系数 ');% 消噪处理for i_decrease=31:44;if abs(h_fre_down(1,i_decrease)>=0.000001h_fre_dow n(1,i_decrease)=(10A-7);endendl_fre_pull=dyadu

11、p(l_fre_down); %0 差值h_fre_pull=dyadup(h_fre_down);l_fre_denoise=conv(low_construct,l_fre_pull);h_fre_denoise=conv(high_construct,h_fre_pull);l_fre_keep=wkeep(l_fre_denoise,L_signal); % 取结果的中心部分 ,消除卷积影响 h_fre_keep=wkeep(h_fre_denoise,L_signal);sig_denoise=l_fre_keep+h_fre_keep; % 消噪后信号重构%平滑处理for j=1:

12、2for i=60:70;sig_denoise(i)=sig_denoise(i-2)+sig_denoise(i+2)/2;end;end;compare=sig_denoise-y; % 与原信号比较figure(3);subplot(3,1,1)plot(y);ylabel('y'); % 原信号subplot(3,1,2);plot(sig_denoise);ylabel('sig_denoise'); % 消噪后信号subplot(3,1,3);plot(compare);ylabel('compare'); % 原信号与消噪后信号的

13、比较 小波谱分析 mallat 算法经典程序 clc;clear;% 1.正弦波定义f1=50; % 频率 1f2=100; % 频率 2fs=2*(f1+f2); % 采样频率 Ts=1/fs; % 采样间隔 N=120; % 采样点数 n=1:N;y=sin(2*pi*f1*n*Ts)+sin(2*pi*f2*n*Ts); % 正弦波混合 figure(1)plot(y);title(' 两个正弦信号 ')figure(2) stem(abs(fft(y); title(' 两信号频谱 ')% 2.小波滤波器谱分析 h=wfilters('db30&

14、#39;,'l'); % 低通 g=wfilters('db30','h'); % 高通 h=h,zeros(1,N-length(h); % 补零(圆周卷积，且增大分辨率变于观察) g=g,zeros(1,N-length(g); % 补零(圆周卷积，且增大分辨率变于观察) figure(3);stem(abs(fft(h);title(' 低通滤波器图 ')figure(4); stem(abs(fft(g);title(' 高通滤波器图 ')% 3.MALLET 分解算法 (圆周卷积的快速傅里叶变换实现)si

15、g1=ifft(fft(y).*fft(h); %低通 (低频分量 )sig2=ifft(fft(y).*fft(g); %高通 (高频分量 )figure(5); % 信号图 subplot(2,1,1)plot(real(sig1); title(' 分解信号 1')subplot(2,1,2) plot(real(sig2);title(' 分解信号 2') figure(6); % 频谱图 subplot(2,1,1) stem(abs(fft(sig1);title('分解信号1频谱')subplot(2,1,2) stem(abs(f

16、ft(sig2);title('分解信号2频谱')% 4.MALLET 重构算法 sig1=dyaddown(sig1); % 2 抽取 sig2=dyaddown(sig2); % 2 抽取 sig1=dyadup(sig1); % 2 插值 sig2=dyadup(sig2); % 2 插值sig1=sig1(1,1:N); % 去掉最后一个零 sig2=sig2(1,1:N); % 去掉最后一个零 hr=h(end:-1:1); % 重构低通 gr=g(end:-1:1); % 重构高通 hr=circshift(hr',1)' % 位置调整圆周右移一位

17、gr=circshift(gr',1)' % 位置调整圆周右移一位 sig1=ifft(fft(hr).*fft(sig1); %低频sig2=ifft(fft(gr).*fft(sig2); %高频sig=sig1+sig2; % 源信号 % 5. 比较 figure(7);subplot(2,1,1) plot(real(sig1); title(' 重构低频信号 ');subplot(2,1,2) plot(real(sig2); title(' 重构高频信号 '); figure(8);subplot(2,1,1) stem(abs(ff

18、t(sig1); title(' 重构低频信号频谱 ');subplot(2,1,2) stem(abs(fft(sig2); title(' 重构高频信号频谱 ');figure(9) plot(real(sig),'r','linewidth',2); hold on;plot(y);legend('重构信号，原始信号') title(' 重构信号与原始信号比较 ')小波包变换分析信号的 MATLA程序%t=0.001:0.001:1;t=1:1000;s1=sin(2*pi*50*t*0.00

19、1)+sin(2*pi*120*t*0.001)+rand(1，length(t);for t=1:500;s2(t)=sin(2*pi*50*t*0.001)+sin(2*pi*120*t*0.001)+rand(1，length(t);endfor t=501:1000;s2(t)=sin(2*pi*200*t*0.001)+sin(2*pi*120*t*0.001)+rand(1，length(t);end subplot(9,2,1)plot(s1)title(' 原始信号 ')ylabel('S1')subplot(9,2,2)plot(s2)titl

20、e(' 故障信号 ') ylabel('S2')wpt=wpdec(s1,3,'db1','shannon');%plot(wpt);s130=wprcoef(wpt,3,0); s131=wprcoef(wpt,3,1); s132=wprcoef(wpt,3,2); s133=wprcoef(wpt,3,3); s134=wprcoef(wpt,3,4); s135=wprcoef(wpt,3,5); s136=wprcoef(wpt,3,6); s137=wprcoef(wpt,3,7); s10=norm(s130); s

21、11=norm(s131);s12=norm(s132); s13=norm(s133); s14=norm(s134); s15=norm(s135); s16=norm(s136); s17=norm(s137);st10=std(s130);st11=std(s131);st12=std(s132);st13=std(s133);st14=std(s134);st15=std(s135);st16=std(s136);st17=std(s137);disp(' 正常信号的特征向量 ');snorm1=s10,s11,s12,s13,s14,s15,s16,s17 std1

22、=st10,st11,st12,st13,st14,st15,st16,st17 subplot(9,2,3);plot(s130);ylabel('S130');subplot(9,2,5);plot(s131);ylabel('S131');subplot(9,2,7);plot(s132);ylabel('S132');subplot(9,2,9);plot(s133);ylabel('S133');subplot(9,2,11);plot(s134);ylabel('S134');subplot(9,2,1

23、3);plot(s135);ylabel('S135');subplot(9,2,15);plot(s136);ylabel('S136');subplot(9,2,17);plot(s137);ylabel('S137');wpt=wpdec(s2,3,'db1','shannon');%plot(wpt);s230=wprcoef(wpt,3,0);s231=wprcoef(wpt,3,1);s232=wprcoef(wpt,3,2);s233=wprcoef(wpt,3,3);s234=wprcoef(wpt

24、,3,4);s235=wprcoef(wpt,3,5);s236=wprcoef(wpt,3,6);s237=wprcoef(wpt,3,7);s20=norm(s230);s21=norm(s231);s22=norm(s232);s23=norm(s233);s24=norm(s234);s25=norm(s235);s26=norm(s236);s27=norm(s237);st20=std(s230);st21=std(s231);st22=std(s232); st23=std(s233);st24=std(s234);st25=std(s235);st26=std(s236);st

25、27=std(s237);disp(' 故障信号的特征向量 '); snorm2=s20,s21,s22,s23,s24,s25,s26,s27 std2=st20,st21,st22,st23,st24,st25,st26,st27 subplot(9,2,4);plot(s230);ylabel('S230');subplot(9,2,6);plot(s231);ylabel('S231');subplot(9,2,8);plot(s232);ylabel('S232');subplot(9,2,10);plot(s233);

26、ylabel('S233');subplot(9,2,12);plot(s234);ylabel('S234');subplot(9,2,14);plot(s235);ylabel('S235');subplot(9,2,16);plot(s236);ylabel('S236');subplot(9,2,18);plot(s237);ylabel('S237');%fft figurey1=fft(s1,1024); py1=y1.*conj(y1)/1024;y2=fft(s2,1024);py2=y2.*con

27、j(y2)/1024;y130=fft(s130,1024); py130=y130.*conj(y130)/1024;y131=fft(s131,1024); py131=y131.*conj(y131)/1024;y132=fft(s132,1024); py132=y132.*conj(y132)/1024;y133=fft(s133,1024);py133=y133.*conj(y133)/1024; y134=fft(s134,1024); py134=y134.*conj(y134)/1024; y135=fft(s135,1024);py135=y135.*conj(y135)/

28、1024; y136=fft(s136,1024);py136=y136.*conj(y136)/1024; y137=fft(s137,1024);py137=y137.*conj(y137)/1024;y230=fft(s230,1024); py230=y230.*conj(y230)/1024;y231=fft(s231,1024); py231=y231.*conj(y231)/1024;y232=fft(s232,1024); py232=y232.*conj(y232)/1024;y233=fft(s233,1024); py233=y233.*conj(y233)/1024;y

29、234=fft(s234,1024); py234=y234.*conj(y234)/1024;y235=fft(s235,1024); py235=y235.*conj(y235)/1024;y236=fft(s236,1024); py236=y236.*conj(y236)/1024;y237=fft(s237,1024); py237=y237.*conj(y237)/1024;f=1000*(0:511)/1024; subplot(1,2,1);plot(f,py1(1:512);ylabel('P1');title(' 原始信号的功率谱 ') su

30、bplot(1,2,2);plot(f,py2(1:512);ylabel('P2');title(' 故障信号的功率谱 ')figuresubplot(4,2,1); plot(f,py130(1:512);ylabel('P130');title('S130 的功率谱 ') subplot(4,2,2);plot(f,py131(1:512);ylabel('P131'); title('S131 的功率谱 ') subplot(4,2,3); plot(f,py132(1:512);ylabe

31、l('P132');subplot(4,2,4); plot(f,py133(1:512);ylabel('P133');subplot(4,2,5); plot(f,py134(1:512);ylabel('P134');subplot(4,2,6); plot(f,py135(1:512);ylabel('P135');subplot(4,2,7); plot(f,py136(1:512);ylabel('P136');subplot(4,2,8); plot(f,py137(1:512);ylabel(

32、9;P137'); figuresubplot(4,2,1); plot(f,py230(1:512);ylabel('P230'); title('S230 的功率谱 ') subplot(4,2,2); plot(f,py231(1:512);ylabel('P231'); title('S231 的功率谱 ') subplot(4,2,3); plot(f,py232(1:512);ylabel('P232');subplot(4,2,4); plot(f,py233(1:512);ylabel(&#

33、39;P233');subplot(4,2,5); plot(f,py234(1:512);ylabel('P234');subplot(4,2,6); plot(f,py235(1:512);ylabel('P235');subplot(4,2,7); plot(f,py236(1:512);ylabel('P236');subplot(4,2,8);plot(f,py237(1:512);ylabel('P237');figure%plottree(wpt)利用小波变换实现对电能质量检测的算法实现N=10000;s=z

34、eros(1,N);for n=1:Nif n<0.4*N|n>0.8*Ns(n)=31.1*sin(2*pi*50/10000*n);elses(n)=22.5*sin(2*pi*50/10000*n);endendl=length(s);c,l=wavedec(s,6,'db5'); % 用 db5 小波分解信号到第六层subplot(8,1,1);plot(s);title(' 用 db5 小波分解六层： s=a6+d6+d5+d4+d3+d2+d1');Ylabel('s');%对分解结构【c, I】中第六层低频部分进行重构a

35、6=wrcoef('a',c,l,'db5',6);subplot(8,1,2);plot(a6);Ylabel('a6');%对分解结构【c, I】中各层高频部分进行重构for i=1:6decmp=wrcoef('d',c,l,'db5',7-i);subplot(8,1,i+2);plot(decmp);Ylabel('d',num2str(7-i);end%rec=zeros(1,300);rect=zeros(1,300);ke=1;u=0;d1=wrcoef('d',c,l,'db5',1);figure(2);plot(d1);si=0;N1=0;N0=0;sce=0;for n=20:N-30rect(ke)=s(n);ke=ke+1;if(ke>=301)if(si=2)rec=rect; u=2;end;si=0;ke=1;end;if(d1(n)>0.01)N1=n;if(N0=0)N0=n; si=si+1;end;if(N1>N0