关于小波分析的matlab程序

合集下载

(参考)小波分析及其matlab实现

(参考)小波分析及其matlab实现

第8章 小波分析及其MATLAB 实现本文档节选自:Matlab 在数学建模中的应用,卓金武等编著, 北航出版社,2011年4月出版8.1小波分析基本理论8.1.1 Fourier 变换的局限性8.1.1.1变换的含义我们把那些定义域和因变域都不是数值或常量的函数称为变换或算子,它们是定义域和值域本身为函数集的函数,如傅里叶变换(Fourier Transform )和拉普拉斯变换(Laplace Transform ),其定义域是时间的函数,而因变域是频率的函数。

简单地说,变换的基本思想仍然是映射,变换是函数的函数。

8.1.1.2 Fourier 变换的局限性信号分析的主要目的是寻找一种简单有效的信号变换方法,使信号包含的重要特征能显示出来。

在小波变换兴起之前,Fourier 级数展开和Fourier 分析是刻画函数空间、求解微分方程、进行数学计算的有效方法和有效的数学工具。

从物理直观上看,一个周期性振动的量可以看成是具有简单频率的简谐振动的叠加。

Fourier 级数展开则是这一物理过程的数学描述,Fourier 变化和Fourier 逆变换公式如下:函数)()(1R L t f ∈的连续Fourier 变换定义为t t f ft d e )()(ˆi -⎰+∞∞-=ωω)(ˆωf 的Fourier 逆变换定义为 ωωωd e )(ˆπ21)(i ⎰+∞∞-=t f t f 从公式上看,Fourier 变换是域变换,它把时间域和频率域联系起来,在时间域上难以观察到的现象和规律,在频域往往能十分清楚地显示出来。

频谱分析的本质就是对)(ˆωf的加工、分析和滤波等处理。

Fourier 变换是平稳信号分析的最重要的工具。

然而在实际运用中,所遇到的信号大多数并不平稳,而是时变频率信号,这时人们需要知道信号在突变时刻所对应的频率成分,显然Fourier 变换的积分作用平滑了非平稳过程的突变部分,作为积分核tω-i e的幅值在任何情况下均为1,因此,在频谱)(ˆωf的任一频点值是由时间过程)(t f 在整个时间域上)(∞+-∞,上的贡献决定的;反之,时间过程)(t f 在某一时刻的状态也是由)(ˆωf在整个频率域)(∞+-∞,上贡献决定的。

关于小波分析的matlab程序

关于小波分析的matlab程序

关于小波分析的matlab程序(个人收集关于小波分析的matlab程序)小波滤波器构造和消噪程序 (2)小波谱分析mallat算法经典程序 (16)小波包变换分析信号的MATLAB程序 (20)利用小波变换实现对电能质量检测的算法实现 (31)基于小波变换的图象去噪Normalshrink算法 (35)小波滤波器构造和消噪程序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;elsecoefficient=1;endhigh_construct(1,i_high)=low_decompose(1,i_hi gh)*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(-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;elsecoefficient=1;endhigh_construct(1,i_high)=low_decompose(1,i_hi gh)*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.000001 h_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'); %原信号与消噪后信号的比较小波谱分析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.00 1)+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));endsubplot(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');%fftfigurey1=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('故障信号的功率谱') 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)); 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'); 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('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: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,l】中第六层低频部分进行重构a6=wrcoef('a',c,l,'db5',6);subplot(8,1,2);plot(a6);Ylabel('a6');%对分解结构【c,l】中各层高频部分进行重构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) % the condition of abnormal append.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)fork=N0:N0+99 %testin g of 1/4 period signals tosce=sce+s(k)*s(k)/10000;end;re=sqrt(sce*200) %re indicate the pike value of .sce=0;si=si+1;end;end;NlenN0n=1:300;figure(3)plot(n,rec);基于小波变换的图象去噪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:levelstemp_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%HLHL=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;%LHLH=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;%HHHH=img((temp_x+1):2*temp_x,(temp_y+1):2*te mp_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*tem p_y)=T_HH;Sub_T(3*(i-1)+3)=T_3;end。

利用dwt实现信号的小波分解 matlab代码

利用dwt实现信号的小波分解 matlab代码

利用dwt实现信号的小波分解 matlab代码
一、简介
小波变换是一种信号处理方法,它通过伸缩和平移小波函数来分析信号的局部特性。

离散小波变换(Discrete Wavelet Transform,DWT)是小波变换的一种形式,它可以用于信号的多尺度分析。

在MATLAB中,我们可以使用内置的小波分析工具箱来实现DWT。

二、代码实现
这段代码首先生成一个简单的正弦信号,然后使用Daubechies 1小波对其进行DWT分解。

最后,它绘制出原始信号以及各尺度的小波系数。

请注意,您需要先安装小波工具箱才能运行此代码。

三、说明与注意事项
此代码仅用于演示DWT的基本概念和实现。

在实际应用中,您可能需要根据具体的需求和信号特性选择合适的小波函数和分解层次。

此外,还需要考虑如何处理噪声和小波系数的阈值压缩等问题。

matlab小波分解程序

matlab小波分解程序

matlab小波分解程序小波分解是一种信号处理的方法,可以用于信号的分析和压缩。

在MATLAB中,可以使用内置的`wavedec`函数来进行小波分解。

下面是一个简单的MATLAB小波分解程序示例:matlab.% 创建一个示例信号。

x = randn(1,1024);% 选择小波基和分解级别。

wname = 'db4'; % 选择小波基,这里使用db4小波。

level = 3; % 选择分解级别。

% 进行小波分解。

[c, l] = wavedec(x, level, wname);% 从分解系数和长度信息中重构近似和细节系数。

appx = wrcoef('a',c,l,wname,level); % 近似系数。

det1 = wrcoef('d',c,l,wname,1); % 第一层细节系数。

det2 = wrcoef('d',c,l,wname,2); % 第二层细节系数。

det3 = wrcoef('d',c,l,wname,3); % 第三层细节系数。

% 绘制原始信号和重构的近似信号。

t = 1:1024;subplot(2,1,1);plot(t, x);title('Original Signal');subplot(2,1,2);plot(t, appx);title('Approximation Coefficients'); % 显示细节系数。

figure;subplot(3,1,1);plot(t, det1);title('Detail Coefficients Level 1'); subplot(3,1,2);plot(t, det2);title('Detail Coefficients Level 2'); subplot(3,1,3);plot(t, det3);title('Detail Coefficients Level 3');在这个示例中,我们首先生成了一个长度为1024的随机信号。

小波分析中的matlab使用

小波分析中的matlab使用

小波分析中的matlab使用Matlab主窗口File菜单File菜单,弹出如图1所示的菜单选项。

其中,各子菜单选项的功能如下:图 1New选项包含5个选项:M-File,Figure,Varible,Model和gui。

1)M-File选项:打开m文件编辑器;2)Figure选项:将打开一个空白的图形窗口;3)Variable选项:可变因素;4)Model选项:用于创建新模型的窗口;5)Gui选项:创建新的图形用户界面的对话框。

Open选项:打开一个open对话框,可以在对话框中选择相应的文件,然后matlab将用相应的编辑器打开该文件。

Close…选项:跟随某个打开的视窗名。

单击该选项,将关闭该视窗。

Importdata…选项:打开一个import对话框,用户可以选择相应的数据文件,然后将该数据文件中的数据导入到matlab工作空间。

Saveworkspaceas…选项:打开一个savetomat-File对话框,用户需要为保存的工作空间命名。

Setpath…选项:打开设置路径对话框。

通过该对话框可以更改matlab执行命令时搜索的路径。

Preferences:首选参数。

Pagesetup选项:用于设置页面布局,页面的页眉,页面所用的文字。

Print…选项:用于打印预定义好的页面内容,也可以设置一些参数。

Printselection…选项:当选中命令窗口内的一部分内容后,该选项将处于激活状态,此时单击该选项,将打印对话框中选中的内容。

Exitmatlab选项:关闭matlab。

也可以通过快捷键ctrl+O来关闭。

Edit菜单单击edit菜单,会弹出如图2所示的菜单选项。

其中,各子菜单选项的功能如下:Undo选项:取消上一次的操作。

Redo选项:重复上一次的操作。

Cut选项:剪切所选中的部分。

Copy选项:选复制被选中的部分。

Paste选项:把存放在缓冲区中的内容粘贴到光标所在的位置。

Pastespecial选项:打开导入数据向导,该向导引导用户把存放在缓冲区中的内容以特定格式存放到剪贴板变量中。

小波变换的原理及matlab仿真程序

小波变换的原理及matlab仿真程序

基于小波变换的信号降噪研究2 小波分析基本理论设Ψ(t)∈L 2( R) ( L 2( R) 表示平方可积的实数空间,即能量有限的信号空间) , 其傅立叶变换为Ψ(t)。

当Ψ(t)满足条件[4,7]:2()Rt dw wCψψ=<∞⎰(1)时,我们称Ψ(t)为一个基本小波或母小波,将母小波函数Ψ(t)经伸缩和平移后,就可以得到一个小波序列:,()()a bt bt aψ-=,,0a b R a ∈≠ (2) 其中a 为伸缩因子,b 为平移因子。

对于任意的函数f(t)∈L 2( R)的连续小波变换为:,(,),()()f a b Rt bW a b f f t dt aψψ-=<>=⎰(3) 其逆变换为:211()(,)()fR R t b f t W a b dadb C a aψψ+-=⎰⎰ (4) 小波变换的时频窗是可以由伸缩因子a 和平移因子b 来调节的,平移因子b,可以改变窗口在相平面时间轴上的位置,而伸缩因子b 的大小不仅能影响窗口在频率轴上的位置,还能改变窗口的形状。

小波变换对不同的频率在时域上的取样步长是可调节的,在低频时,小波变换的时间分辨率较低,频率分辨率较高:在高频时,小波变换的时间分辨率较高,而频率分辨率较低。

使用小波变换处理信号时,首先选取适当的小波函数对信号进行分解,其次对分解出的参数进行阈值处理,选取合适的阈值进行分析,最后利用处理后的参数进行逆小波变换,对信号进行重构。

3 小波降噪的原理和方法3.1 小波降噪原理从信号学的角度看 ,小波去噪是一个信号滤波的问题。

尽管在很大程度上小波去噪可以看成是低通滤波 ,但由于在去噪后 ,还能成功地保留信号特征 ,所以在这一点上又优于传统的低通滤波器。

由此可见 ,小波去噪实际上是特征提取和低通滤波的综合 ,其流程框图如图所示[6]:小波分析的重要应用之一就是用于信号消噪 ,一个含噪的一维信号模型可表示为如下形式:(k)()()S f k e k ε=+* k=0.1…….n-1其中 ,f( k)为有用信号,s(k)为含噪声信号,e(k)为噪声,ε为噪声系数的标准偏差。

小波分析MATLAB实例

小波分析MATLAB实例

小波分析MATLAB实例小波分析是一种信号处理方法,可以用于信号的时频分析和多尺度分析。

在MATLAB中,可以使用Wavelet Toolbox实现小波分析。

这个工具箱提供了丰富的函数和工具,可以方便地进行小波分析的计算和可视化。

小波分析的核心是小波变换,它将信号分解成一组不同尺度和频率的小波基函数。

在MATLAB中,可以使用`cwt`函数进行连续小波变换。

以下是一个小波分析的MATLAB实例,用于分析一个心电图信号的时频特性。

首先,导入心电图信号数据。

假设心电图数据保存在一个名为`ecg_signal.mat`的文件中,包含一个名为`ecg`的变量。

可以使用`load`函数加载这个数据。

```MATLABload('ecg_signal.mat');```接下来,设置小波变换的参数。

选择一个小波基函数和一组尺度。

这里选择Morlet小波作为小波基函数,选择一组从1到64的尺度。

可以使用`wavelet`函数创建一个小波对象,并使用`scal2frq`函数将尺度转换为频率。

```MATLABwavelet_name = 'morl'; % 选择Morlet小波作为小波基函数scales = 1:64; % 选择1到64的尺度wavelet_obj = wavelet(wavelet_name);scales_freq = scal2frq(scales, wavelet_name, 1);```然后,使用`cwt`函数进行小波变换,得到信号在不同尺度和频率下的小波系数。

将小波系数的幅度平方得到信号的能量谱密度。

```MATLAB[wt, f] = cwt(ecg, scales, wavelet_name);energy = abs(wt).^2;``````MATLABimagesc(1:length(ecg), scales_freq, energy);colormap('jet');xlabel('时间(样本)');ylabel('频率(Hz)');```运行整个脚本之后,就可以得到心电图信号的时频图。

第六章基于MATLAB的小波分析

第六章基于MATLAB的小波分析

第六章基于MATLAB的小波分析小波分析是一种用来分析和处理信号的数学方法,其基本原理是通过将信号分解成不同频率范围的小波基函数来揭示信号的特征。

MATLAB是一种功能强大的科学计算和数据分析软件,提供了丰富的工具箱和函数,可以方便地进行小波分析。

在MATLAB中,小波分析可以通过使用Wavelet Toolbox来实现。

该工具箱提供了几种常用的小波基函数,如Daubechies、Coiflets、Symlets等,同时还包括了一系列小波分析的函数。

下面将介绍基于MATLAB的小波分析的基本步骤。

首先,需要导入待分析的信号数据。

可以使用MATLAB的数据导入和处理工具来加载信号数据,如load函数、importdata函数等。

加载数据后,可以使用plot函数将信号数据可视化,以便直观地了解信号的特点。

接下来,需要选择合适的小波基函数进行分析。

小波基函数的选择与信号的特征和分析目标相关。

可以使用waveinfo函数来查看Wavelet Toolbox提供的小波基函数的特性和参数,并选择适合的小波基函数。

然后,使用wavedec函数对信号进行小波分解。

wavedec函数可以将信号分解成多个尺度的小波系数。

分解得到的小波系数包括近似系数和细节系数,近似系数反映了信号在低频范围的特征,而细节系数则反映了信号在高频范围的细节特征。

分解后,可以使用可视化函数如plot、imshow等来展示小波系数的分布和变化情况。

通过观察小波系数的变化,可以得到信号的频率特征和局部特征。

除了观察小波系数,还可以根据需要进行小波系数的处理和分析。

例如,可以使用细节系数来提取信号中的细节特征,如边缘、尖峰等,也可以使用近似系数来提取信号的整体趋势。

最后,可以使用waverec函数将处理过的小波系数重构成原始信号。

重构得到的信号可以与原始信号进行对比,以验证分析的结果和提取的特征。

综上所述,MATLAB提供了丰富的工具和函数来实现小波分析,可以方便地进行信号的频率分析和特征提取。

Daubechies小波分析的部分Matlab实现

Daubechies小波分析的部分Matlab实现

实验5:Daubechies 小波分析的部分Matlab 实现注:本实验在原本操作过程中在重构部分存在着问题,经思考检验调试之后,排除了隐藏的问题,实验得到成功进行同时达到了预期的结果。

第一部分Daubechies 尺度序列,二进点上的尺度函数图像及小波函数图像第一部分实验的目标有两个:1,给出各阶Daubechies 小波的低通滤波器,即尺度函数两尺度序列。

具体的程序算法由《小波导论》中7.3节中的引理7.16,定理7.17给现,其中定理7.17的证明过程实际上给出了一个显式的计算过程,这是程序实现的重要依据。

2,根据两尺度关系,由第一步骤中得到的两尺度序列逐步计算得到Daubechies 尺度函数及小波函数的二进点值,进而给出近似图像。

第一步骤相关函数程序解释:1、sinj(n),实现将引理7.16中的2(sin)2n ω转化为形如cos k ω,n 即表示(sin)2ω的次数。

2、sinj(n)中涉及函数cosn(n), cosn(n)实现将cos kω转化为cos()kkw ∑,n 表示三角函数次数。

程序编写的思路就是先将2(sin)2n ω转化为cos k ω再转化为cos()kkw ∑,然后根据定理7.17的证明过程计算出两尺度序列,其中用于提取多项式系数的函数利用了三角函数的正交性,使用了积分函数int()。

daucoef(n)给出n 阶daubechies 尺度函数的两尺度序列。

具体的原程序代码见附件,部分程序运行结果如下:各阶(2--12)Daubechies 尺度函数的两尺度序列如下:)第二步骤:给出n 阶daubechies 尺度函数图像,由函数dauinterp(n,max1)实现,即给出n 阶daubechies 尺度函数的max1次迭代的函数值,所得到的点都是函数的二进点。

Dauinterpwa(n,max1)则给出小波函数的二进点上的函数值。

这一部分内容中,比较难以实现的地方在于下标的变换,这常常引起一些数据的混乱或错误,相关的数学推导总是难以避免的,由于过程过于繁复,这里不给出详细的推导过程及解释。

实验三 MATLAB小波实现

实验三 MATLAB小波实现

实验三 MATLAB 小波实现实验题目:利用MATLAB 实现一些小波函数,并绘图实验内容:1. Haar 小波函数(1) 定义函数源程序:function Haarfor x=0:0.001:2if x>=0&x<=0.5H=1;elseif x>=0.5&x<=1H=-1;elseH=0;endplot(x,H,'k')axis([0,2,-2,2]) title('Harr 定义函数')hold onendhold offend运行程序得图:⎪⎪⎩⎪⎪⎨⎧<≤-≤≤=其它012112/101x x H ψ(2) 尺度函数源程序:function harrfor x=0:0.001:2if x>=0&x<=1H=1;elseH=0;endplot(x,H,'k')axis([0,2,-2,2]) title('Harr 尺度函数')hold onendhold offend00.20.40.60.81 1.2 1.4 1.6 1.82Harr 定义函数⎩⎨⎧≤≤=其它0101)(x x φ运行程序得图:2. Mexican Hat 小波函数函数定义为:源程序:function MexicanHat[psi,x]=mexihat(-5,5,1000);plot(x,psi,'k') title('Mexican Hat 小波')end运行程序得图:00.20.40.60.81 1.2 1.4 1.6 1.82Harr 尺度函数2/24/12e )1(π32)(x x x ---=ψ3. Meyer 小波函数(1)定义函数v (a )为构造Meyer 小波的辅助函数,且有v (a )=a 4(35-84a +70a 2-20a 3),a ∈[0,1](2)尺度函数-5-4-3-2-1012345-0.4-0.20.20.40.60.81Mexican Hat 小波⎪⎪⎪⎩⎪⎪⎪⎨⎧∉≤≤-≤≤-=--]3π8,3π2[,03π83π4)),1π23(2πcos(e π)2(3π43π2)),1π23(2πsin(e π)2()(ˆ2/i 2/12/i 2/1ωωωνωωνωψωω⎪⎪⎪⎩⎪⎪⎪⎨⎧>≤≤-≤=--3π403π43π2))1π23(2πcos(π)2(3π2π)2()(2/12/1ωωωνωωΦ源程序:function Meyer[phi,psi,x]=meyer(-8,8,1024);figure(1);plot(x,psi,'k'),title('Meyer 定义函数') figure(2);plot(x,phi,'k'),title('Meyer 尺度函数') end运行程序得图:-8-6-4-202468-1-0.50.511.5Meyer 定义函数4. Morlet 小波函数函数定义为:源程序:function Morlet[psi,x]=morlet(-4,4,1000);plot(x,psi,'k')title('MorletС²¨')end运行程序得图:-8-6-4-202468-0.4-0.20.20.40.60.811.2Meyer 尺度函数)5cos(e )(2/2x C x x -=ψ5. Gauss 函数族 同一坐标系中ααπα4221)(x e x g -= 161,41,1=α 及其频域形式 源程序:function Gaussx=-5:0.001:5;g1=1/(2*pi^(1/2))*exp(-x.^2/4);g2=1/(pi^(1/2))*exp(-x.^2);g3=2/(pi^(1/2))*exp(-4*x.^2);plot(x,g1,x,g2,x,g3);title('Gauss 函数族')end运行程序得图:-4-3-2-101234-1-0.8-0.6-0.4-0.20.20.40.60.81Morlet 小波-5-4-3-2-101234500.20.40.60.811.21.4Gauss 函数族。

matlab wsst 实现方法

matlab wsst 实现方法

matlab wsst 实现方法使用MATLAB实现WSST方法引言:小波分析是一种在时间和频率域上进行信号分析的有效工具。

小波分析可以将信号分解成不同频率的成分,使得对信号的分析更加全面和准确。

其中,WSST(Wavelet Synchrosqueezed Transform)方法是一种基于小波分析的信号处理方法,可以用于时频分析、频谱估计和信号特征提取等领域。

本文将介绍如何使用MATLAB实现WSST方法,并通过一个示例来展示其应用。

一、MATLAB中的小波分析工具MATLAB提供了丰富的小波分析工具箱,可以方便地进行小波变换、小波重构和小波分析等操作。

在实现WSST方法之前,我们首先需要了解MATLAB中的小波分析工具。

1. 小波变换小波变换是一种将信号分解成不同频率的成分的方法。

MATLAB中的小波变换函数为“wavetrans”。

通过选择不同的小波基函数和尺度参数,可以得到不同频率的小波系数。

2. 小波重构小波重构是一种将小波系数合成为原始信号的方法。

MATLAB中的小波重构函数为“iwavetrans”。

通过将不同频率的小波系数进行合成,可以得到原始信号的近似重构。

3. 小波分析工具箱MATLAB提供了丰富的小波分析工具箱,包括小波变换、小波重构、小波包分析、小波阈值去噪等功能。

通过使用这些工具,可以方便地进行小波分析和信号处理。

二、WSST方法的原理WSST方法是一种基于小波分析的信号处理方法,可以将信号在时频域上进行分析。

其原理是通过对信号进行小波变换,然后对小波系数进行重构,得到信号在时频域上的表示。

WSST方法可以提取信号的时频特征,从而实现对信号的分析和处理。

1. 小波变换我们需要对信号进行小波变换。

通过选择合适的小波基函数和尺度参数,可以将信号分解成不同频率的小波系数。

MATLAB中的小波变换函数为“wavetrans”。

2. 小波重构然后,我们需要对小波系数进行重构,得到信号在时频域上的表示。

个人收集的一些关于小波分析的matlab程序

个人收集的一些关于小波分析的matlab程序

个人收集的一些关于小波分析的matlab程序都是从网上收集来的,在文档结构图里查看,每个标题对应一个程序,难免重复,请大家理解。

标题: 提升法97经典程序1.%% 本程序实现任意偶数大小图像第二代双正交97提升小波变换%% 注1:采用标准正交方法,对行列采用不同矩阵(和matlab里不同)2.%% 注2:为了保证正交,所有边界处理,全部采用循环处理3.%% 注3:正交性验证,将单位阵带入函数,输出仍是单位阵(matlab不具有此性质)4.%% 注4:此程序是矩阵实现,所以图像水平分量和垂直分量估计被交换位置5.%% 注5:此程序实现的是类小波(wavelet-like)变换,是介于小波包变换与小波变换之间的变换6.%% 注6:此程序每层变换相对原图像矩阵,产生的矩阵都是正交阵,这和小波包一致7.%% 注7:但小波变换每层产生的矩阵,是相对每个待分解子块的正交矩阵,而不是原图像的正交矩阵8.%% 注8:且小波变换产生的正交矩阵维数,随分解层数2分减少9.%% 注9: 提升系数可以在MATLAB7.0以上版本,用liftwave('9.7')获取,这里直接给出,考虑兼容性10.%% 注10:由于MATLAB数组下标从1开始,所以注意奇偶序列的变化11.%% 注11:d为对偶上升,即预测;p为原上升,即更新%%编程人沙威安徽大学12.%% 编程时间2004年12月18日%% x输入图像,y输出图像13.%% flag_trans为正变换或反变换标志,0执行正变换,1执行反变换14.%% flag_max,是否最大层数变换标志,0执行用户设定层数,1执行最大层数变换15.%% layer,用户层数设置(小于最大层)functiony=db97(x,flag_trans,flag_max,layer); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%16.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%17.18.% 1.输入参数检查% 矩阵维数判断19.[sa,sb]=size(x); if (sa~=sb) % 防止非图像数据20.errordlg('非图像数据!');21.error('非图像数据!');22.end; % 变换标志判断23.[sa,sb]=size(flag_trans);24.if ((sa~=1) | (sb~=1)) % 变换标志错误25.errordlg('变换标志错误!');26.error('变换标志错误!');27.end; if ((flag_trans~=1) & (flag_trans~=0)) % 变换标志错误28.errordlg('变换标志错误!');29.error('变换标志错误!');30.end; % 最大层数标志判断31.[sa,sb]=size(flag_max);32.if ((sa~=1) | (sb~=1)) % 最大层数标志错误33.errordlg('最大层数标志错误!');34.error('最大层数标志错误!');35.end; if ((flag_max~=1) & (flag_max~=0)) % 最大层数标志错误36.errordlg('最大层数标志错误!');37.error('最大层数标志错误!');38.end; % 用户设置层数判断39.if (flag_max~=1) [sa,sb]=size(layer);40.if ((sa~=1) | (sb~=1)) % 层数设置错误41.errordlg('层数设置错误!');42.error('层数设置错误!');43.end; if (flag_max<0) % 层数设置错误44.errordlg('层数设置错误!');45.error('层数设置错误!');46.end;47.end;48.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 49.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%50.51.% 2.提升系数确定52.% t1=liftwave('9.7'); % 获取提升系数(MATLAB7.0以后)d1=[-1.586100000000000e+000,-1.586134342069360e+00 0];53.p1=[1.079600000000000e+000,-5.298011857188560e-002];54.d2=[-8.829110755411875e-001,-8.829110755411875e-001];55.p2=[4.435068*********e-001,1.576123746148364e+000];56.d3=-8.698644516247808e-001;57.p3=-1.149604398860242e+000;58.59.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 60.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%61.62.% 3.分解层数确定63.% 采用用户输入和自动给出最大层数两种方法N=length(x); % 矩阵大小64.S=N; % 变量65.s=log2(N); % 最大循环次数66.n1=N/2; % 初始一半矩阵大小67.n2=N; % 初始矩阵大小68.u=0; % 初始值% 对非2的整数幂大小图像确定最大分解层数69.for ss=1:s70.if (mod(S,2)==0)71.u=u+1;72.S=S/2;73.end;74.end;75.u=u-1; % 分解最大层数减1(后面的边界处理造成) % 最大层数确定76.if (flag_max==0) % 手动输入77.T=layer; % 用户输入值78.else % 自动确定最大层数79.T=u; % 分解最大层数80.end81.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 82.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%83.84.% 4.最大层数和图像大小检查if (T>u) % 防止用户层数越界85.errordlg('已超过最大分解层数!或者非偶数大小图像!');86.error('已超过最大分解层数!或者非偶数大小图像!');87.end; if (mod(N,2)~=0) % 防止图像大小错误88.errordlg('非偶数大小图像!');89.error('非偶数大小图像!');90.end;91.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 92.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%93.94.% 5.提升法正变换if (flag_trans==0)95.for time=1:T; % 行正变换96.97.% d;98.x1(n1,:)=x(n2,:)+d1(2)*x(n2-1,:)+d1(1)*x(1,:);99.x1([1:n1-1],:)=x([2:2:n2-2],:)+d1(2)*x([1:2:n2-3],:)+d1(1)*x([3:2:n2-1],:);100.101.% p;102.x(1,:)=x(1,:)+p1(2)*x1(n1,:)+p1(1)*x1(1,:);103.x([2:n1],:)=x([3:2:n2-1],:)+p1(2)*x1([1:n1-1],:)+p1(1)* x1([2:n1],:);104.x([n1+1:n2],:)=x1([1:n1],:);105.106.% d;107.x(n1+1,:)=x(n1+1,:)+d2(2)*x(n1,:)+d2(1)*x(1,:); 108.x([n1+2:n2],:)=x([n1+2:n2],:)+d2(2)*x([1:n1-1],:)+d2(1 )*x([2:n1],:);109.110.% p;111.x(n1,:)=x(n1,:)+p2(2)*x(n1+1,:)+p2(1)*x(n1+2,:); 112.x(n1-1,:)=x(n1-1,:)+p2(2)*x(n2,:)+p2(1)*x(n1+1,:); 113.x([1:n1-2],:)=x([1:n1-2],:)+p2(2)*x([n1+2:n2-1],:)+p2(1 )*x([n1+3:n2],:);114.115.% 归一116.x([1:n1],:)=p3*x([1:n1],:);117.x([n1+1:n2],:)=d3*x([n1+1:n2],:); clear x1;118.119.% 列正变换120.121.% d;122.x1(:,[1:n1])=x(:,[2:2:n2]);123.124.% p;125.x(:,1)=x(:,1)-d1(1)*x1(:,n1)-d1(2)*x1(:,1);126.x(:,[2:n1])=x(:,[3:2:n2-1])-d1(1)*x1(:,[1:n1-1])-d1(2)*x1 (:,[2:n1]);127.x(:,[n1+1:n2])=x1(:,[1:n1]);128.129.% d;130.x(:,n2)=x(:,n2)-p1(1)*x(:,n1)-p1(2)*x(:,1);131.x(:,[n1+1:n2-1])=x(:,[n1+1:n2-1])-p1(1)*x(:,[1:n1-1])-p 1(2)*x(:,[2:n1]);132.133.% p;134.x(:,n1,:)=x(:,n1)-d2(1)*x(:,n2)-d2(2)*x(:,n1+1);135.x(:,[1:n1-1])=x(:,[1:n1-1])-d2(1)*x(:,[n1+1:n2-1])-d2(2) *x(:,[n1+2:n2]);136.137.% d;138.x(:,n1+1)=x(:,n1+1)-p2(1)*x(:,n1-1)-p2(2)*x(:,n1); 139.x(:,n1+2)=x(:,n1+2)-p2(1)*x(:,n1)-p2(2)*x(:,1);140.x(:,[n1+3:n2])=x(:,[n1+3:n2])-p2(1)*x(:,[1:n1-2])-p2(2) *x(:,[2:n1-1]);141.142.% 归一143.x(:,[1:n1])=d3*x(:,[1:n1]);144.x(:,[n1+1:n2])=p3*x(:,[n1+1:n2]); clear x1;145.146.n2=n2/2; % 原大小147.n1=n2/2; % 一半大小148.end; 149.%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%% 150.%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%151.152.% 6.提升法反变换else153.n2=N/(2.^(T-1)); % 分解最小子块维数154.n1=n2/2;155.for time=1:T; % 行反变换156.157.% 去归一158.x([1:n1],:)=x([1:n1],:)/p3;159.x([n1+1:n2],:)=x([n1+1:n2],:)/d3; % 反p;160.x(n1,:)=x(n1,:)-p2(2)*x(n1+1,:)-p2(1)*x(n1+2,:); 161.x(n1-1,:)=x(n1-1,:)-p2(2)*x(n2,:)-p2(1)*x(n1+1,:);162.x([1:n1-2],:)=x([1:n1-2],:)-p2(2)*x([n1+2:n2-1],:)-p2(1) *x([n1+3:n2],:);163.164.% 反d;165.x(n1+1,:)=x(n1+1,:)-d2(2)*x(n1,:)-d2(1)*x(1,:);166.x([n1+2:n2],:)=x([n1+2:n2],:)-d2(2)*x([1:n1-1],:)-d2(1) *x([2:n1],:);167.168.% 反p;169.x1(1,:)=x(1,:)-p1(2)*x(n2,:)-p1(1)*x(n1+1,:);170.x1([2:n1],:)=x([2:n1],:)-p1(2)*x([n1+1:n2-1],:)-p1(1)*x( [n1+2:n2],:);171.172.% 反d;173.x(n2,:)=x(n2,:)-d1(2)*x1(n1,:)-d1(1)*x1(1,:);174.x([2:2:n2-2],:)=x([n1+1:n2-1],:)-d1(2)*x1([1:n1-1],:)-d1(1)*x1([2:n1],:);175.176.% 偶数177.x([1:2:n2-1],:)=x1([1:n1],:);178.179.clear x1;180.181.% 列反变换182.183.% 归一184.x(:,[1:n1])=x(:,[1:n1])/d3;185.x(:,[n1+1:n2])=x(:,[n1+1:n2])/p3; % 反d;186.x(:,n1+1)=x(:,n1+1)+p2(1)*x(:,n1-1)+p2(2)*x(:,n1); 187.x(:,n1+2)=x(:,n1+2)+p2(1)*x(:,n1)+p2(2)*x(:,1); 188.x(:,[n1+3:n2])=x(:,[n1+3:n2])+p2(1)*x(:,[1:n1-2])+p2(2 )*x(:,[2:n1-1]);189.190.% 反p;191.x(:,n1,:)=x(:,n1)+d2(1)*x(:,n2)+d2(2)*x(:,n1+1); 192.x(:,[1:n1-1])=x(:,[1:n1-1])+d2(1)*x(:,[n1+1:n2-1])+d2(2 )*x(:,[n1+2:n2]);193.194.% 反d;195.x(:,n2)=x(:,n2)+p1(1)*x(:,n1)+p1(2)*x(:,1);196.x(:,[n1+1:n2-1])=x(:,[n1+1:n2-1])+p1(1)*x(:,[1:n1-1])+ p1(2)*x(:,[2:n1]);197.198.% 反p;199.x1(:,1)=x(:,1)+d1(1)*x(:,n2)+d1(2)*x(:,n1+1);200.x1(:,[2:n1])=x(:,[2:n1])+d1(1)*x(:,[n1+1:n2-1])+d1(2)* x(:,[n1+2:n2]); % 奇偶201.x(:,[2:2:n2])=x(:,[n1+1:n2]);202.x(:,[1:2:n2-1])=x1(:,[1:n1]); clear x1;203.204.n2=n2*2; % 原大小205.n1=n2/2; % 一半大小end;206.end;207. 208.%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%% 209.%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%210.211.% 7.结果输出y=x;212.% 传输最后结果%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%% 213.%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%214.215.% 8.内存清理clear x;216.clear flag_max;217.clear layer;218.clear flag_trans;219.clear N;220.clear n1;221.clear n2;222.clear s;223.clear ss;224.clear u;225.clear d1;226.clear d2;227.clear d3;228.clear p1;229.clear p2;230.clear p3;231.clear sa;232.clear sb;233. 234.%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%复制代码标题: 2代小波示意程序1.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%2.%% 此程序用提升法实现第二代小波变换3.%% 我用的是非整数阶小波变换4.%% 采用时域实现,步骤先列后行5.%% 正变换:分裂,预测,更新;6.%% 反变换:更新,预测,合并7.%% 只做一层(可以多层,而且每层的预测和更新方程不同)clear;clc; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%8.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%9.%%%%10.11.% 1.调原始图像矩阵load wbarb; % 下载图像12.f=X; % 原始图像13.% f=[0 0 0 0 0 0 0 0 ;...14.% 0 0 0 1 1 0 0 0 ;...15.% 0 0 2 4 4 2 0 0 ;...16.% 0 1 4 8 8 4 1 0 ;...17.% 0 1 4 8 8 4 1 0 ;...18.% 0 0 2 4 4 2 0 0 ;...19.% 0 0 0 1 1 0 0 0 ;...20.% 0 0 0 0 0 0 0 0 ;]; % 原始图像矩阵N=length(f); % 图像维数21.T=N/2;22.23.% 子图像维数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%24.%%%%%%%%%%%%%%%%%正变换%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%% 26.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%27.%%%% 1.列变换28.29.% A.分裂(奇偶分开)f1=f([1:2:N-1],:); % 奇数30.f2=f([2:2:N],:); % 偶数% f1(:,T+1)=f1(:,1); % 补列31.% f2(T+1,:)=f2(1,:); % 补行% B.预测for i_hc=1:T;32.high_frequency_column(i_hc,:)=f1(i_hc,:)-f2(i_hc,:);33.end; %high_frequency_column(T+1,:)=high_frequency_column(1, :); % 补行% C.更新for i_lc=1:T;34.low_frequency_column(i_lc,:)=f2(i_lc,:)+1/2*high_frequency_column(i_lc,:);35.end; % D.合并36.f_column([1:1:T],:)=low_frequency_column([1:T],:);37.f_column([T+1:1:N],:)=high_frequency_column([1:T],:);38.39.40.figure(1)41.colormap(map);42.image(f); figure(2)43.colormap(map);44.image(f_column);45.46.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 47.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%48.%%%% 2.行变换49.50.% A.分裂(奇偶分开)f1=f_column(:,[1:2:N-1]); % 奇数51.f2=f_column(:,[2:2:N]); % 偶数52.% f2(:,T+1)=f2(:,1); % 补行% B.预测for i_hr=1:T;53.high_frequency_row(:,i_hr)=f1(:,i_hr)-f2(:,i_hr);54.end; %high_frequency_row(:,T+1)=high_frequency_row(:,1); % 补行% C.更新for i_lr=1:T;55.low_frequency_row(:,i_lr)=f2(:,i_lr)+1/2*high_frequency_row(:,i_lr);56.end; % D.合并57.f_row(:,[1:1:T])=low_frequency_row(:,[1:T]);58.f_row(:,[T+1:1:N])=high_frequency_row(:,[1:T]);59.figure(3)60.colormap(map);61.image(f_row);62.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 63.%%%%%%%%%%%%%%%%%%%反变换%%%%%%%%%%%%%%%%%%%%%%%%% 64.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 65.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 66.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%67.68.%%%% 1.行变换69.% A.提取(低频高频分开)f1=f_row(:,[T+1:1:N]); % 奇数70.f2=f_row(:,[1:1:T]); % 偶数71.% f2(:,T+1)=f2(:,1); % 补行% B.更新for i_lr=1:T;72.low_frequency_row(:,i_lr)=f2(:,i_lr)-1/2*f1(:,i_lr);73.end; % C.预测for i_hr=1:T;74.high_frequency_row(:,i_hr)=f1(:,i_hr)+low_frequency_row(:,i_hr);75.end; %high_frequency_row(:,T+1)=high_frequency_row(:,1); % 补行76.% D.合并(奇偶分开合并)77.f_row(:,[2:2:N])=low_frequency_row(:,[1:T]);78.f_row(:,[1:2:N-1])=high_frequency_row(:,[1:T]);79.figure(4)80.colormap(map);81.image(f_row);82.83.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 84.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%85.%%%% 2.列变换86.87.% A.提取(低频高频分开)f1=f_row([T+1:1:N],:); % 奇数88.f2=f_row([1:1:T],:); % 偶数% f1(:,T+1)=f1(:,1); % 补列89.% f2(T+1,:)=f2(1,:); % 补行% B.更新for i_lc=1:T;90.low_frequency_column(i_lc,:)=f2(i_lc,:)-1/2*f1(i_lc,:);91.end; % C.预测for i_hc=1:T;92.high_frequency_column(i_hc,:)=f1(i_hc,:)+low_frequency_column(i_hc,:);93.end; %high_frequency_column(T+1,:)=high_frequency_column(1,:); % 补行% D.合并(奇偶分开合并)94.f_column([2:2:N],:)=low_frequency_column([1:T],:);95.f_column([1:2:N-1],:)=high_frequency_column([1:T],:);96.97.98.figure(5)99.colormap(map);100.image(f_column);复制代码标题: 二代小波漫谈现在我就举例,对一个8点序列,怎样实现第二代小波变换。

Matlab中的小波分析工具箱

Matlab中的小波分析工具箱

欢迎下载 可修改
13
其他的一维函数:
抽样: dyaddow 补零插值:dyaup 滤波器生成:qmf,orthfilt,wfilters 反变换:idwt,idwtper, 重构: upwlev,waverec,wrcoef,
欢迎下载 可2
[cA,cH,cV,cD]=dwt2(X,’wname’)
Ten lectures on wavelets,
CBMS, SIAM, 61, 1994, 117-119, 137, 152.
欢迎下载 可修改
5
计算小波滤波器系数的函数:
参数表示
morlet mexihat meyer meyeraux dbwavf dbaux symwavf coifwavf biowavf
欢迎下载 可修改
21
LastWave
小波信号和图像处理软件,用C语言编写, 可在Unix和Macintosh上运行。
下载地址:
欢迎下载 可修改
22
值得关注的几个发展方向:
提升小波变换(Lifting scheme wavelet transform)
多小波变换(Multiwavelet transform)
小波基的名称
计算Morlet小波滤波器系数
计算墨西哥草帽小波滤波器系数
计算Meyer小波与尺度滤波器系数
计算Meyer小波辅助函数
计算紧支集双正交小波滤波器系数
计算紧支集双正交小波尺度滤波器系数
计算近似对称的紧支集双正交小波滤波器系数
计算Coifmant小波尺度滤波器系数
计算双正交样条小波尺度滤波器参数
DWT
possible but without FWT

小波变换matlab程序

小波变换matlab程序

小波变换matlab程序小波变换是一种信号处理技术,它可以将信号分解成不同频率的成分,并且可以在不同时间尺度上进行分析。

在Matlab中,可以使用内置的小波变换函数来实现这一技术。

下面是一个简单的小波变换Matlab程序示例:matlab.% 生成一个示例信号。

t = 0:0.001:1; % 时间范围。

f1 = 10; % 信号频率。

f2 = 50; % 信号频率。

y = sin(2pif1t) + sin(2pif2t); % 信号。

% 进行小波变换。

[c, l] = wavedec(y, 3, 'db1'); % 进行3层小波分解,使用db1小波基函数。

% 重构信号。

yrec = waverec(c, l, 'db1'); % 使用小波系数和长度进行信号重构。

% 绘制原始信号和重构信号。

subplot(2,1,1);plot(t, y);title('原始信号');subplot(2,1,2);plot(t, yrec);title('重构信号');这个程序首先生成了一个包含两个频率成分的示例信号,然后使用`wavedec`函数对信号进行小波分解,得到小波系数和长度。

接着使用`waverec`函数对小波系数和长度进行信号重构,最后绘制了原始信号和重构信号的对比图。

小波变换在信号处理、图像处理等领域有着广泛的应用,可以用于信号去噪、特征提取、压缩等方面。

通过Matlab中的小波变换函数,我们可以方便地进行小波分析和处理,从而更好地理解和利用信号的特性。

小波分析Matlab程序

小波分析Matlab程序

1 绪论1.1概述小波分析是近15年来发展起来的一种新的时频分析方法。

其典型应用包括齿轮变速控制,起重机的非正常噪声,自动目标所顶,物理中的间断现象等。

而频域分析的着眼点在于区分突发信号和稳定信号以及定量分析其能量,典型应用包括细胞膜的识别,金属表面的探伤,金融学中快变量的检测,INTERNET 的流量控制等。

从以上的信号分析的典型应用可以看出,时频分析应用非常广泛,涵盖了物理学,工程技术,生物科学,经济学等众多领域,而且在很多情况下单单分析其时域或频域的性质是不够的,比如在电力监测系统中,即要监控稳定信号的成分,又要准确定位故障信号。

这就需要引入新的时频分析方法,小波分析正是由于这类需求发展起来的。

在传统的傅立叶分析中,信号完全是在频域展开的,不包含任何时频的信息,这对于某些应用来说是很恰当的,因为信号的频率的信息对其是非常重要的。

但其丢弃的时域信息可能对某些应用同样非常重要,所以人们对傅立叶分析进行了推广,提出了很多能表征时域和频域信息的信号分析方法,如短时傅立叶变换,Gabor 变换,时频分析,小波变换等。

其中短时傅立叶变换是在傅立叶分析基础上引入时域信息的最初尝试,其基本假定在于在一定的时间窗内信号是平稳的,那么通过分割时间窗,在每个时间窗内把信号展开到频域就可以获得局部的频域信息,但是它的时域区分度只能依赖于大小不变的时间窗,对某些瞬态信号来说还是粒度太大。

换言之,短时傅立叶分析只能在一个分辨率上进行。

所以对很多应用来说不够精确,存在很大的缺陷。

而小波分析则克服了短时傅立叶变换在单分辨率上的缺陷,具有多分辨率分析的特点,在时域和频域都有表征信号局部信息的能力,时间窗和频率窗都可以根据信号的具体形态动态调整,在一般情况下,在低频部分(信号较平稳)可以采用较低的时间分辨率,而提高频率的分辨率,在高频情况下(频率变化不大)可以用较低的频率分辨率来换取精确的时间定位。

因为这些特定,小波分析可以探测正常信号中的瞬态,并展示其频率成分,被称为数学显微镜,广泛应用于各个时频分析领域。

小波分析的应用及其MATLAB程序的实现

小波分析的应用及其MATLAB程序的实现

小波分析的应用及其MATLAB 程序的实现 摘要:在简单介绍小波分析的发展的基础上,对傅立叶变换和小波变换比较分析,介绍了小波分析在实际生活中的应用,重点阐述了MA 的应用研究现存的几个TLAB 小波分析信号处理的方法.分析了小波分析在故障诊断中问题,并对解决这些问题和未来的发展进行了探讨。

关键词:小波分析;信号处理;MATLAB1.引言故障诊断中的首要问题就是对观测信号的故障特征提取,即对观测信号进行信号处理,从中获取反映故障信息的特征。

由于故障诊断中所遇到的信号绝大多数都是非平稳信号,而特别适用于非平稳信号处理的工具就是小波分析,所以小波分析在故障诊断中的应用越来越受到人们的青睬。

小波变换的基本思想类似于傅立叶变换,小波分析优于博立叶之处在于它能够实现时域和频域的局部分析,即通过伸缩和平移等运算功能对函数或信号进行多尺度细化分析,从而可以聚焦到信号的任意细节。

因此,小波变换被誉为分析信号的微镜。

现在,小波分析技术在信号处理、图像处理、语音分析、模识别、量子物理、生物医学工程、计算机视觉、故障诊断及众多非线性科学领域都有广泛的应用。

2、从傅立叶变换到小波变换小波分析属于时频分析的一种,传统的信号分析是建立在傅立叶变换的基础上的,由于傅立叶分析使用的是一种全局的变换,要么完全在时域,要么完全在时域,要么完全在频域,因此无法表述信号的时频局域性质,而这种性质恰恰是非平稳信号最根本和最关键的性质。

为了分析和处理非平稳信号,人们对傅立叶分析进行了推广乃至根本性的革命,提出并发展了一系列新的信号分析理论:短时傅立叶变换、Gabor 变换、时频分析、小波变换、分数阶傅立叶变换、线调频小波变换、循环统计量理论和调幅-调频信号分析等。

其中,短时傅立叶变换和小波变换也是应传统的傅立叶变换不能够满足信号处理的要求而产生的。

短时傅立叶变换分析的基本思想是:假定非平稳信号在分析窗函数g (t )的一个短时间间隔内是平稳(伪平稳)的,并移动分析窗函数,使)()(τ-t g t f 在不同的有限时间宽度内是平稳信号,从而计算出各个不同时刻的功率谱。

整数小波变换的matlab程序

整数小波变换的matlab程序

整数小波变换的matlab程序整数小波变换的matlab程序是一种用于对信号进行分析和处理的工具。

它可以将连续的信号转化为离散的信号,以实现对信号的不同频率组分的提取和处理。

在本文中,我们将一步一步地介绍整数小波变换的Matlab 程序。

首先,我们需要明确整数小波变换的概念。

整数小波变换是一种将信号分解为不同频率下的子信号和低频近似部分的方法。

通过将信号与特定的小波基函数进行卷积,可以得到不同频率分量的系数。

整数小波变换不仅能够提取信号的频率信息,还能提供时间和频率的局部化。

在Matlab中实现整数小波变换,首先需要加载信号。

在这个例子中,我们将使用经典的ECG信号作为输入信号。

可以使用以下代码将ECG信号加载到Matlab环境中:load('ecg_data.mat');接下来,我们可以定义小波基函数。

在整数小波变换中,通常使用的小波基函数有Haar小波、Symlet小波和Daubechies小波。

在这个例子中,我们选择使用Daubechies小波,可以使用以下代码定义:wname = 'db4';接下来,我们需要将信号进行分解。

Matlab提供了一个名为`wavedec`的函数,可以用于对信号进行离散小波变换。

以下是分解信号的代码:[coeff, l] = wavedec(ecg_data, N, wname);在这个代码中,`wavedec`函数的第一个参数是输入信号,第二个参数是分解的层数N,第三个参数是小波基函数的名称。

`wavedec`函数将返回分解系数和信号的长度l。

接下来,我们可以提取不同频率分量的系数。

通过对分解系数进行适当的索引,可以得到不同频率分量的系数。

以下是提取不同频率系数的代码:cA = appcoef(coeff, l, wname, 1);cD1 = detcoef(coeff, l, 1);cD2 = detcoef(coeff, l, 2);在这个代码中,`appcoef`函数是用于提取低频近似系数的函数,`detcoef`函数是用于提取详细系数的函数。

关于小波分析的matlab程序

关于小波分析的matlab程序

关于小波分析的matlab程序小波分析是一种在信号处理和数据分析领域中广泛应用的方法。

它可以帮助我们更好地理解信号的时域和频域特性,并提供一种有效的信号处理工具。

在本文中,我将介绍小波分析的基本原理和如何使用MATLAB编写小波分析程序。

一、小波分析的基本原理小波分析是一种基于窗口函数的信号分析方法。

它使用一组称为小波函数的基函数,将信号分解成不同频率和不同时间尺度的成分。

与傅里叶分析相比,小波分析具有更好的时频局部化性质,可以更好地捕捉信号的瞬时特征。

小波函数是一种具有局部化特性的函数,它在时域上具有有限长度,并且在频域上具有有限带宽。

常用的小波函数有Morlet小波、Haar小波、Daubechies小波等。

这些小波函数可以通过数学运算得到,也可以通过MATLAB的小波函数库直接调用。

小波分析的基本步骤如下:1. 选择合适的小波函数作为基函数。

2. 将信号与小波函数进行卷积运算,得到小波系数。

3. 根据小波系数的大小和位置,可以分析信号的时频特性。

4. 根据需要,可以对小波系数进行阈值处理,实现信号的去噪和压缩。

二、MATLAB中的小波分析工具MATLAB提供了丰富的小波分析工具箱,可以方便地进行小波分析的计算和可视化。

下面介绍几个常用的MATLAB函数和工具箱:1. `waveinfo`函数:用于查看和了解MATLAB中可用的小波函数的信息,如小波函数的名称、支持的尺度范围等。

2. `wavedec`函数:用于对信号进行小波分解,得到小波系数。

3. `waverec`函数:用于根据小波系数重构原始信号。

4. `wdenoise`函数:用于对小波系数进行阈值处理,实现信号的去噪。

5. 小波分析工具箱(Wavelet Toolbox):提供了更多的小波分析函数和工具,如小波变换、小波包分析、小波阈值处理等。

可以通过`help wavelet`命令查看工具箱中的函数列表。

三、编写小波分析程序在MATLAB中编写小波分析程序可以按照以下步骤进行:1. 导入信号数据:首先需要导入待分析的信号数据。

使用MATLAB小波工具箱进行小波分析:

使用MATLAB小波工具箱进行小波分析:

使用MATLAB小波工具箱进行小波分析:如上图所示的小波分解过程,可以调用wfilters 来获得指定小波的分解和综合滤波器系数,例如:% Set wavelet name.wname = 'db5';% Compute the four filters associated with wavelet name given% by the input string wname.[Lo_D,Hi_D,Lo_R,Hi_R] = wfilters(wname);subplot(221); stem(Lo_D);title('Decomposition low-pass filter');subplot(222); stem(Hi_D);title('Decomposition high-pass filter');subplot(223); stem(Lo_R);title('Reconstruction low-pass filter');subplot(224); stem(Hi_R);title('Reconstruction high-pass filter');xlabel('The four filters for db5')% Editing some graphical properties,% the following figure is generated.以上例子,得到’db5’小波的分解和综合滤波器系数,并显示出来。

下面是wfilters的具体用法:Wname 可指定为列表中的任意一种小波,直接调用[Lo_D,Hi_D,Lo_R,Hi_R] = wfilters('wname')会返回分解和综合滤波器系数。

如果只想返回其中的一些而不是全部,可以调用[F1,F2] = wfilters('wname','type')其中’type’可指定为4种类型,每种类型的具体意义详见matlab wfilters帮助。

相关主题
  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

/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;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(-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;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'); %原信号与消噪后信号的比较小波谱分析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));endfor t=501:1000;s2(t)=sin(2*pi*200*t*0.001)+sin(2*pi*120*t*0.001)+rand(1,length(t));endsubplot(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');%fftfigurey1=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('故障信号的功率谱') figuresubplot(4,2,1);plot(f,py130(1:512));ylabel('P130');title('S130的功率谱') subplot(4,2,2);plot(f,py131(1:512));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'); 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('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));subplot(4,2,8);plot(f,py237(1:512));ylabel('P237');figure%plottree(wpt)利用小波变换实现对电能质量检测的算法实现N=10000;s=zeros(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,l】中第六层低频部分进行重构a6=wrcoef('a',c,l,'db5',6);subplot(8,1,2);plot(a6);Ylabel('a6');%对分解结构【c,l】中各层高频部分进行重构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) % the condition of abnormal append.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 %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;NlenN0n=1:300;figure(3)plot(n,rec);基于小波变换的图象去噪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:levelstemp_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%HLHL=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;%LHLH=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;%HHHH=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。

相关文档
最新文档