当前位置:首页 >> 电力/水利 >>

关于小波分析的matlab程序


http://www.chinavib.com/forum-viewthread-tid-9141-extra-page%3D1-page-1.html(个人收 集关于小波分析的 matlab 程序) 小波滤波器构造和消噪程序 .............................................................................

.............................. 1 小波谱分析 mallat 算法经典程序 ................................................................................................... 7 小波包变换分析信号的 MATLAB 程序 ........................................................................................ 9 利用小波变换实现对电能质量检测的算法实现 ......................................................................... 15 基于小波变换的图象去噪 Normalshrink 算法............................................................................ 17

小波滤波器构造和消噪程序
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'); %原信号与消噪后信号的比较 2. 消噪、 % 此函数用于研究 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; 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(-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.000001 h_fre_down(1,i_decrease)=(10^-7); end end 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; %消噪后信号重构 %平滑处理 for j=1:2 for 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'); %原信号与消噪后信号的比较

小波谱分析 mallat 算法经典程序
clc;clear; %% 1.正弦波定义 f1=50; % 频率 1 f2=100; % 频率 2

fs=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,s17]

std1=[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,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); 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 figure y1=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)

利用小波变换实现对电能质量检测的算法 实现
N=10000; s=zeros(1,N); for n=1:N if n<0.4*N||n>0.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) N1=n; if(N0==0) N0=n; si=si+1; end; if(N1>N0+30) Nlen=N1-N0; Tab=Nlen/10000; end; end; if(si==1) for k=N0:N0+99 sce=sce+s(k)*s(k)/10000; end; re=sqrt(sce*200) sce=0; si=si+1; end; end; Nlen N0 n=1:300; figure(3) plot(n,rec);

% the condition of abnormal append.

%testing of 1/4 period signals to

%re indicate the pike value of .

基于小波变换的图象去噪 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/2^i; temp_y=yy/2^i; % 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


相关文章:
小波分析的应用及其MATLAB程序的实现
小波分析的应用及其 MATLAB 程序的实现摘要:在简单介绍小波分析的发展的基础上,对傅立叶变换和小波变换比较 摘要 分析,介绍了小波分析在实际生活中的应用,重点阐述...
matlab小波周期分析程序
matlab小波周期分析程序_IT/计算机_专业资料。小波分析程序 clear s=xlsread('1') for i=1:length(s); for j=1:length(s); wreal(i,j)=0.0; wimage...
小波变换的原理及matlab仿真程序_图文
小波变换的原理及matlab仿真程序_计算机软件及应用_IT/计算机_专业资料。巨有用,初学者全靠这份文档 基于小波变换的信号降噪研究 2 小波分析基本理论设 Ψ(t)∈L...
小波分析MATLAB实例
4.2 小波去噪的模型建立 4.2.1 去噪的 Matlab 程序 局部放电试验所采集的信号...关于小波分析的matlab程... 18页 免费 小波分析中的matlab使用 9页 5下载券 ...
matlab小波去噪详解超全超全有程序
matlab小波去噪详解超全超全有程序_数学_自然科学_专业资料。小波去噪 [xd,...基于小波去噪matlab程序... 3页 1下载券 小波去噪方法分析与MATL... 2页...
小波分析中MATLAB阈值获取函数及其应用附程序代码
1、小波分析MATLAB 阈值获取函数 MATLAB实现阈值获取的函数有 ddencmp、 thselect、 wbmpen 和 wwdcbm, 下面对它们的用法进行简单的说明。 一、ddencmp 的...
小波变换MATLAB程序
小波变换MATLAB程序_工学_高等教育_教育专区。故障诊断 小波变换的matlab程序 这个程序首先画出原函数 f=sin(0.03*t)的图像,然后画出原函数加上噪音的图像,在 ...
五种常见小波基函数及其matlab实现
五种常见小波基函数及其matlab实现_信息与通信_工程科技_专业资料。与标准的傅里叶变换相比,小波分析中使用到的小波函数具有不唯一性,即小波函数 具有多样性。小波...
用matlab小波分析的实例
文中给出了详细的程序范例,用 MATLAB 实现了基于小波变换的图 像处理。 1.2...域才发挥它独特的优势,解决了各类问题,为人们提供了更多的关于时域分析的信息。...
小波信号分解与重构的Matlab程序
以下是一维小波分解的程序Matlab 小波分析工具箱丰富的函数和强大的仿真功能为我们学习小波、 用好 小波提供了方便、快捷的途径,但是,如果我们要深入掌握小波分析的...
更多相关标签: