关于小波分析的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;
else
coefficient=1;
end
high_construct(1,i_high)=low_decompose(1,i_hi gh)*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_hi gh)*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.00 1)+rand(1,length(t));
end
for t=501:1000;
s2(t)=sin(2*pi*200*t*0.001)+sin(2*pi*120*t*0.0
01)+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) % 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 %testin g 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);
基于小波变换的图象去噪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*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。