matlab仿真程序.doc

4搜集整理小波分析的matlab程序4.1小波滤波器构造和消噪程序重构% mallet_wavelet.m % 此函数用于研究Mallet算法及滤波器设计 % 此函数仅用于消噪 a=pi/8; %角度赋初值 b=pi/8; %低通重构FIR滤波器h0(n)冲激响应赋值 h0=cos(a)*cos(b); h1=sin(a)*cos(b); h2=-sin(a)*sin(b); h3=cos(a)*sin(b); low_construct=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)n,高通重建滤波器 if(mod(i_high,2)=0); coefficient=-1; else coefficient=1; end high_construct(1,i_high)=low_decompose(1,i_high)*coefficient; end high_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(-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_decompose); %卷积 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=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; %信号重构 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); %原信号与消噪后信号的比较消噪、% 此函数用于研究Mallet算法及滤波器设计% 此函数用于消噪处理%角度赋值%此处赋值使滤波器系数恰为db9%分解的高频系数采用db9较好,即它的消失矩较大%分解的有用信号小波高频系数基本趋于零%对于噪声信号高频分解系数很大，便于阈值消噪处理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)n,高通重建滤波器if(mod(i_high,2)=0);coefficient=-1;elsecoefficient=1;endhigh_construct(1,i_high)=low_decompose(1,i_high)*coefficient;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(原信号);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_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(小波分解的高频系数);% 消噪处理for i_decrease=31:44;if abs(h_fre_down(1,i_decrease)=0.000001h_fre_down(1,i_decrease)=(10-7);endendl_fre_pull=dyadup(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: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); %原信号与消噪后信号的比较4.2小波谱分析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,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分解算法(圆周卷积的快速傅里叶变换实现)sig1=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(fft(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); % 位置调整圆周右移一位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(fft(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(重构信号与原始信号比较)小波包变换分析信号的MATLAB程序%t=0.001:0.001:1;t=1:1000;s1=sin(2*pi*50*t*0.001)+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);end for 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)title(故障信号)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);s11=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,s17std1=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,13);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,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);st27=std(s237); disp(故障信号的特征向量);snorm2=s20,s21,s22,s23,s24,s25,s26,s27std2=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);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.*conj(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)/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;y234=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(原始信号的功率谱)subplot(1,2,2);plot(f,py2(1:512);ylabel(P2);title(故障信号的功率谱) figure subplot(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);ylabel(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(P137);figure subplot(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(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)4.3利用小波变换实现对电能质量检测的算法实现N=10000;s=zeros(1,N);for n=1:N if n0.8*N s(n)=31.1*sin(2*pi*50/10000*n); else s(n)=22.5*sin(2*pi*50/10000*n); end end l=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，l】中第六层低频部分进行重构a6=wrcoef(a,c,l,db5,6);subplot(8,1,2);plot(a6);Ylabel(a6);%对分解结构【c，l】中各层高频部分进行重构for i=1:6 decmp=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-30 rect(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) % the condition of abnormal append. N1=n; if(N0=0) N0=n; si=si+1; end; if(N1N0+30) Nlen=N1-N0; Tab=Nlen/10000; end; end; if(si=1) for k=N0:N0+99 %testing of 1/4 period signals to sce=sce+s(k)*s(k)/10000; end; re=sqrt(sce*200) %re indicate the pike value of . sce=0; si=si+1; end;end; Nlen N0 n=1:300; figure(3) plot(n,rec); 4.4基于小波变换的图象去噪 Normalshrink算法function T_img,Sub_T=threshold_2_N(img,levels) % reference :image denoising using wavelet thresholding xx,yy=size(img);HH=img(xx/2+1):xx,(yy/2+1):yy); delt_2=(std(HH(:)2;%(median(abs(HH(:)/0.6745)2;% T_img=img; for i=1:levels temp_x=xx/2i; temp_y=yy/2i;% belt=1.0*(log(temp_x/(2*levels)0.5; belt=1.0*(log(temp_x/(2*levels)0.5; %2.5 0.8 %HL HL=img(1:temp_x,(temp_y+1):2*temp_y); delt_y=std(HL(:); T_1=belt*delt_2/delt_y; % T_HL=sign(HL).*max(0,abs(HL)-T_1); % T_img(1:temp_x,(temp_y+1):2*temp_y)=T_HL; Sub_T(3*(i-1)+1)=T_1; %LH LH=img(temp_x+1):2*temp_x,1:temp_y); delt_y=std(LH(:); T_2=belt*delt_2/delt_y; % T_LH=sign(LH).*max(0,abs(LH)-T_2); % T_img(temp_x+1):2*temp_x,1:temp_y)=T_LH; Sub_T(3*(i-1)+2)=T_2; %HH HH=img(temp_x+1):2*temp_x,(temp_y+1):2*temp_y); delt_y=std(HH(:); T_3=belt*delt_2/delt_y; % T_HH=sign(HH).*max(0,abs(HH)-T_3); % T_img(temp_x+1):2*temp_x,(temp_y+1):2*temp_y)=T_HH; Sub_T(3*(i-1)+3)=T_3;end 4.5小波去噪matlab程序 * clear clc %在噪声环境下语音信号的增强 %语音信号为读入的声音文件 %噪声为正态随机噪声 sound=wavread(c12345.wav); count1=length(sound); noise=0.05*randn(1,count1); for i=1:count1 signal(i)=sound(i); end for i=1:count1 y(i)=signal(i)+noise(i); end %在小波基db3下进行一维离散小波变换 coefs1,coefs2=dwt(y,db3); %低频 高频 count2=length(coefs1); count3=length(coefs2); energy1=sum(abs(coefs1).2); energy2=sum(abs(coefs2).2); energy3=energy1+energy2; for i=1:count2 recoefs1(i)=coefs1(i)/energy3; end for i=1:count3 recoefs2(i)=coefs2(i)/energy3; end %低频系数进行语音信号清浊音的判别 zhen=160; count4=fix(count2/zhen); for i=1:count4 n=160*(i-1)+1:160+160*(i-1); s=sound(n); w=hamming(160); sw=s.*w; a=aryule(sw,10); sw=filter(a,1,sw); sw=sw/sum(sw); r=xcorr(sw,biased); corr=max(r); %为清音（unvoice）时，输出为1；为浊音（voice）时，输出为0 if corr=0.8 output1(i)=0; elseif corr=0.1 output1(i)=1; end end for i=1:count4 n=160*(i-1)+1:160+160*(i-1); if output1(i)=1 switch abs(recoefs1(i) case abs(recoefs1(i)0.002 & abs(recoefs1(i)=0.8 output2(i)=0; elseif corr=0.1 output2(i)=1; end end for i=1:count5 n=160*(i-1)+1:160+160*(i-1); if output2(i)=1 switch abs(recoefs2(i) case abs(recoefs2(i)0.002 & abs(recoefs2(i)=0.003 recoefs2(i)=sgn(recoefs2(i)*(0.003*abs(recoefs2(i)-0.000003)/0.002; otherwise recoefs2(i)=recoefs2(i); end elseif output2(i)=0 recoefs2(i)=recoefs2(i); end end %在小波基db3下进行一维离散小波反变换 output3=idwt(recoefs1, recoefs2,db3); %对输出信号抽样点值进行归一化处理 maxdata=max(output3); output4=output3/maxdata; %读出带噪语音信号，存为101.wav wavwrite(y,5500,16,c101); %读出处理后语音信号，存为102.wav wavwrite(output4,5500,16,c102);4.6 神经网络小波Matlab实现%clear allload dataM1=20;epo=15;A=4;B=18;B2=B/2+1N=500;M=(A+1)*(B+1);for a0=1:A+1;for b0=1:B+1;i=(B+1)*(a0-1)+b0;b_init(i)=(b0-B2)/10)/(2(-A); a_init(i)=1/(2(-A);c_init(i)=(20-A)/2;end end w0=ones(1,M); for i=1:N for j=1:M t=x(i); t= a_init(j)*t-b_init(j); %P0(i,j)= (cos(1.75*t)*exp(-t*t/2)/2c_init(j); P0(i,j)= (1-t*t)*exp(-t*t/2)/2c_init(j); end end%calculation of output of network for i=1: