mymallat.m

[原创]自己写的mallat算法matlab源程序

clc
clear
f1=1;f2=20;f3=40;fs=200;N=400;
omiga=2*pi/fs;
t=0:N-1;
s1=sin(f1*omiga*t);
s2=sin(f2*omiga*t);
s3=sin(f3*omiga*t);
s=s1+s2+s3;
wavelet='sym8';
[Lo_D,Hi_D,Lo_R,Hi_R] = wfilters(wavelet);
figure(2);
stem(abs(fft(Lo_D)));
title('低通滤波器图')
figure(3);
stem(abs(fft(Hi_D)));
title('高通滤波器图')
% A0=zeros(1,N+16);
% A0(8:1:N+7)=s;
%A0 = s;
figure(1)
subplot(521);plot(s);Ylabel('x(t)')
%%%%%%%%%%%%%%%%%%%%%%%%分解%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      Dtemp = conv(Hi_D,s);
      Atemp = conv(Lo_D,s);
      Atemp= wkeep(Atemp,N);
      Dtemp= wkeep(Dtemp,N);
      subplot(522)
      plot(t,Atemp);Ylabel('a1')
      subplot(523)
      plot(t,Dtemp);Ylabel('d1')
      d1=Dtemp(2:2:N);
      a1=Atemp(2:2:N);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      t=0:N/2-1;
      Dtemp = conv(Hi_D,a1);     %卷积
      Atemp = conv(Lo_D,a1);
      Atemp= wkeep(Atemp,N/2);   %将卷积得到的数据进行截断处理
      Dtemp= wkeep(Dtemp,N/2);
      subplot(524)
      plot(t,Atemp);Ylabel('a2')
      subplot(525)
      plot(t,Dtemp);Ylabel('d2')
      d2=Dtemp(2:2:N/2);
      a2=Atemp(2:2:N/2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      t=0:N/4-1;
      Dtemp = conv(Hi_D,a2);
      Atemp = conv(Lo_D,a2);
      Atemp= wkeep(Atemp,N/4);
      Dtemp= wkeep(Dtemp,N/4);
      subplot(526)
      plot(t,Atemp);Ylabel('a3')
      subplot(527)
      plot(t,Dtemp);Ylabel('d3')
      d3=Dtemp(2:2:N/4);
      a3=Atemp(2:2:N/4);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      t=0:N/8-1;
      Dtemp = conv(Hi_D,a3);
      Atemp = conv(Lo_D,a3);
      Atemp= wkeep(Atemp,N/8);
      Dtemp= wkeep(Dtemp,N/8);
      subplot(528)
      plot(t,Atemp);Ylabel('a4')
      subplot(529)
      plot(t,Dtemp);Ylabel('d4')
      d4=Dtemp(2:2:N/8);
      a4=Atemp(2:2:N/8);
%%%%%%%%%%%%%%%%%%重构%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      a=zeros(1,N/8);
      d=zeros(1,N/8);
      a(1:2:N/8)=a4;
      d(1:2:N/8)=d4;
      Dtemp = conv(Hi_R,d);
      Atemp = conv(Lo_R,a);
      w= (Dtemp+Atemp);
      a3= wkeep(w,N/8);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      a=zeros(1,N/4);
      d=zeros(1,N/4);
      a(1:2:N/4)=a3;
      d(1:2:N/4)=d3;
      Dtemp = conv(Hi_R,d);
      Atemp = conv(Lo_R,a);
      w= (Dtemp+Atemp);
      a2= wkeep(w,N/4);
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      a=zeros(1,N/2);
      d=zeros(1,N/2);
      a(1:2:N/2)=a2;
      d(1:2:N/2)=d2;
      Dtemp = conv(Hi_R,d);
      Atemp = conv(Lo_R,a);
      w= (Dtemp+Atemp);
      a1= wkeep(w,N/2);
      a=zeros(1,N);
      d=zeros(1,N);
      a(1:2:N)=a1;
      d(1:2:N)=d1;
      Dtemp = conv(Hi_R,d);
      Atemp = conv(Lo_R,a);
      w= (Dtemp+Atemp);
      a0= wkeep(w,N);
      t=0:N-1;
      figure(4)
      subplot(211); plot(t,s);
      subplot(212); plot(t,a0); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%重构d3，频率为20的信号
% d=zeros(1,N/4);
% d(1:2:N/4)=d3;
% Dtemp = conv(Hi_R,d);
% d2= wkeep(Dtemp,N/4);
% d=zeros(1,N/2);
% d(1:2:N/2)=d2;
% Dtemp = conv(Lo_R,d);
% d1= wkeep(Dtemp,N/2);
% d=zeros(1,N);
% d(1:2:N)=d1;
% Dtemp = conv(Lo_R,d);
% d0= wkeep(Dtemp,N);
% figure(4)
% plot(d0);ylabel('d3');title('20Hz重构信号');
% figure(5);plot(s2,'r');hold on;plot(d0,'b');title('20Hz重构信号与原始信号s2比较')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%重构a3，频率为1的信号
% a=zeros(1,N/4);
% a(1:2:N/4)=a3;
% Atemp = conv(Lo_R,a);
% a2= wkeep(Atemp,N/4);
% a=zeros(1,N/2);
% a(1:2:N/2)=a2;
% Atemp = conv(Lo_R,a);
% a1= wkeep(Atemp,N/2);
% a=zeros(1,N);
% a(1:2:N)=a1;
% Atemp = conv(Lo_R,a);
% a0= wkeep(Atemp,N);
% figure(4)
% plot(a0); ylabel('a3')
% figure(5);plot(s1,'r');hold on;plot(a0,'b');title('重构信号与原始信号比较')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%重构d2，频率为40的信号
% d=zeros(1,N/2);
% d(1:2:N/2)=d2;
% Dtemp = conv(Hi_R,d);
% d1= wkeep(Dtemp,N/2);
% d=zeros(1,N);
% d(1:2:N)=d1;
% Dtemp = conv(Lo_R,d);
% d0= wkeep(Dtemp,N);
% figure(4)
% plot(d0);ylabel('d2');title('20Hz重构信号');
% figure(5);plot(s3,'r');hold on;plot(d0,'b');title('40Hz重构信号与原始信号s3比较');