多小波 nir.m

用于数据预处理的多小波分解与重构程序.其中多小波是采用CL4平衡多小波,所以分解前的预处理是采用等值预滤波器.如果要用别的非平衡多小波,只需替换为相应的预滤波器即可.

%调入待分解的数据
clear
load nirbeer;
xh=[Xcal;Xtest];
x=wextend('ac','sp1',xh,49,'b');

%预处理 
%平衡多小波的等值预滤波器
m=1;
p(:,:,1)=[1 0;0 1];
q(:,:,1)=[1 0;0 1];
%预处理(将原始一维数据x转化为二维数据c0)
n=length(x);
%对x前后延拓,前后各增加2*(m-1)个，共增加4*(m-1)个,这里共增加4个 8个。
z=zeros(60,(2*m-2));
x=[z x z];
N=(n+2*m-2)/2;%延拓数据后N=6 7 c0有2*6个  2*7个
c0=zeros(2,N,60);% 2*n(6)个0元素  2*n(7)个0元素
for j=1:60
for k=1:N
   for i=1:m
      c0(:,k,j)=c0(:,k,j)+p(:,:,i)*[x(j,(2*(k+i)-3));x(j,(2*(k+i)-2))];
   end
end
end

%分解
%平衡算子CL4多小波
M=4;
h(:,:,1)=sqrt(2)*[10-3*sqrt(10),5*sqrt(6)-2*sqrt(15);5*sqrt(6)-3*sqrt(15),5-3*sqrt(10)]/80;
h(:,:,2)=sqrt(2)*[30+3*sqrt(10),5*sqrt(6)-2*sqrt(15);-5*sqrt(6)-7*sqrt(15),15-3*sqrt(10)]/80;
s=[1 0;0 -1];
h(:,:,3)=s*h(:,:,2)*s;
h(:,:,4)=s*h(:,:,1)*s;
g(:,:,1)=sqrt(2)*[5*sqrt(6)-2*sqrt(15) -10+3*sqrt(10);-5+3*sqrt(10) 5*sqrt(6)-3*sqrt(15)]/80;
g(:,:,2)=sqrt(2)*[-5*sqrt(6)+2*sqrt(15) 30+3*sqrt(10);15-3*sqrt(10) 5*sqrt(6)+7*sqrt(15)]/80;
g(:,:,3)=s*g(:,:,2)*s;
g(:,:,4)=s*g(:,:,1)*s;
r=[1 -1;1 1]/sqrt(2);%平衡算子
h(:,:,1)=r*h(:,:,1)*r';
h(:,:,2)=r*h(:,:,2)*r';
h(:,:,3)=r*h(:,:,3)*r';
h(:,:,4)=r*h(:,:,4)*r';
g(:,:,1)=r*g(:,:,1)*r';
g(:,:,2)=r*g(:,:,2)*r';
g(:,:,3)=r*g(:,:,3)*r';
g(:,:,4)=r*g(:,:,4)*r';

%一级分解(由c0分解得到c1和d1)
N=length(c0(:,:,1));
yushu1=mod(N,2);
if yushu1~=0
c0_=[c0,c0(:,end,:)];
else c0_=c0;
end
N=length(c0_(:,:,1));
c0_=[c0_,c0_(:,1,:),c0_(:,2,:)];
c1=zeros(2,N/2,60);
for j=1:60
for i=1:(N/2)
    for k=1:M
         c1(:,i,j)= c1(:,i,j)+h(:,:,k)*c0_(:,(2*i+k-2),j);
     end
 end
end
for i=1:60
  c1ti(:,:,i)=wkeep(c1(:,:,i),[2,fix(length(xh)/4)]);
end

%二级分解
N=length(c1(:,:,1));
yushu1=mod(N,2);
if yushu1~=0
c1_=[c1,c1(:,end,:)];
else c1_=c1;
end
N=length(c1_(:,:,1));
c1_=[c1_,c1_(:,1,:),c1_(:,2,:)];
c2=zeros(2,N/2,60);
for j=1:60
for i=1:(N/2)
    for k=1:M
         c2(:,i,j)= c2(:,i,j)+h(:,:,k)*c1_(:,(2*i+k-2),j);
     end
 end
end
for i=1:60
  c2ti(:,:,i)=wkeep(c2(:,:,i),[2,fix(length(xh)/8)]);
end


%三级分解(由c2_分解得到c3和d3)
N=length(c2(:,:,1));
yushu1=mod(N,2);
if yushu1~=0
c2_=[c2,c2(:,end,:)];
else c2_=c2;
end
N=length(c2_(:,:,1));
c2_=[c2_,c2_(:,1,:),c2_(:,2,:)];
c3=zeros(2,N/2,60);
for j=1:60
for i=1:(N/2)
    for k=1:M
         c3(:,i,j)= c3(:,i,j)+h(:,:,k)*c2_(:,(2*i+k-2),j);
     end
 end
end
for i=1:60
  c3ti(:,:,i)=wkeep(c3(:,:,i),[2,fix(length(xh)/16)]);
end

%四级分解(由c3_分解得到c4和d4)
N=length(c3(:,:,1));
yushu4=mod(N,2);
if yushu4~=0
c3_=[c3,c3(:,end,:)];
else 
c3_=c3;    
end
N=length(c3_(:,:,1));
c3_=[c3_,c3_(:,1,:),c3_(:,2,:)];
c4=zeros(2,N/2,60);
for j=1:60
for i=1:(N/2)
    for k=1:M
         c4(:,i,j)= c4(:,i,j)+h(:,:,k)*c3_(:,(2*i+k-2),j);
     end
 end
end
for i=1:60
  c4ti(:,:,i)=wkeep(c4(:,:,i),[2,fix(length(xh)/32)]);
end
%5级分解(由c4_分解得到c5和d5)
N=length(c4(:,:,1));
yushu5=mod(N,2);
if yushu5~=0
c4_=[c4,c4(:,end,:)];
else 
c4_=c4;    
end
N=length(c4_(:,:,1));
c4_=[c4_,c4_(:,1,:),c4_(:,2,:)];
c5=zeros(2,N/2,60);
for j=1:60
for i=1:(N/2)
    for k=1:M
         c5(:,i,j)= c5(:,i,j)+h(:,:,k)*c4_(:,(2*i+k-2),j);
     end
 end
end
for i=1:60
  c5ti(:,:,i)=wkeep(c5(:,:,i),[2,fix(length(xh)/64)]);
end

%6级分解(由c5_分解得到c6和d6)
N=length(c5(:,:,1));
yushu6=mod(N,2);
if yushu6~=0
c5_=[c5,c5(:,end,:)];
else 
c5_=c5;    
end
N=length(c5_(:,:,1));
c5_=[c5_,c5_(:,1,:),c5_(:,2,:)];
c6=zeros(2,N/2,60);
for j=1:60
for i=1:(N/2)
    for k=1:M
         c6(:,i,j)= c6(:,i,j)+h(:,:,k)*c5_(:,(2*i+k-2),j);
     end
 end
end
for i=1:60
  c6ti(:,:,i)=wkeep(c6(:,:,i),[2,fix(length(xh)/128)]);
end

%7级分解(由c6_分解得到c7和d7)
N=length(c6(:,:,1));
yushu7=mod(N,2);
if yushu7~=0
c6_=[c6,c6(:,end,:)];
else 
c6_=c6;    
end
N=length(c6_(:,:,1));
c6_=[c6_,c6_(:,1,:),c6_(:,2,:)];
c7=zeros(2,N/2,60);
for j=1:60
for i=1:(N/2)
    for k=1:M
         c7(:,i,j)= c7(:,i,j)+h(:,:,k)*c6_(:,(2*i+k-2),j);
     end
 end
end
for i=1:60
  c7ti(:,:,i)=wkeep(c7(:,:,i),[2,fix(length(xh)/256)]);
end

c7r=c7;
%7级重构(由c7r重构c6r)
N=length(c7r(:,:,1));
c7r_=[c7r(:,N,:),c7r];
c6r=zeros(2,(2*N),60);
for j=1:60
for i=1:(2*N)
    for k=(2-mod(i,2)):2:M
    c6r(:,i,j)=c6r(:,i,j)+h(:,:,k)'*c7r_(:,((i-k+4)/2),j);
    end
end
end
c6r=c6;
%6级重构(由c6r重构c5r)
N=length(c6r(:,:,1));
c6r_=[c6r(:,N,:),c6r];
c5r=zeros(2,(2*N),60);
for j=1:60
for i=1:(2*N)
    for k=(2-mod(i,2)):2:M
    c5r(:,i,j)=c5r(:,i,j)+h(:,:,k)'*c6r_(:,((i-k+4)/2),j);
    end
end
end

c5r=c5;
%5级重构(由c5r重构c4r)
N=length(c5r(:,:,1));
c5r_=[c5r(:,N,:),c5r];
c4r=zeros(2,(2*N),60);
for j=1:60
for i=1:(2*N)
    for k=(2-mod(i,2)):2:M
    c4r(:,i,j)=c4r(:,i,j)+h(:,:,k)'*c5r_(:,((i-k+4)/2),j);
    end
end
end

c4r=c4;
%四级重构(由c4r重构c3r)
N=length(c4r(:,:,1));
c4r_=[c4r(:,N,:),c4r];
c3r=zeros(2,(2*N),60);
for j=1:60
for i=1:(2*N)
    for k=(2-mod(i,2)):2:M
    c3r(:,i,j)=c3r(:,i,j)+h(:,:,k)'*c4r_(:,((i-k+4)/2),j);
    end
end
end

c3r=c3;
%三级重构(由c3r重构c2r)
N=length(c3r(:,:,1));
c3r_=[c3r(:,N,:),c3r];
c2r=zeros(2,(2*N),60);
for j=1:60
for i=1:(2*N)
    for k=(2-mod(i,2)):2:M
    c2r(:,i,j)=c2r(:,i,j)+h(:,:,k)'*c3r_(:,((i-k+4)/2),j);
    end
end
end

c2r=c2;
%二级重构(由c2r重构c1r)
N=length(c2r(:,:,1));
c2r_=[c2r(:,N,:),c2r];
c1r=zeros(2,(2*N),60);
for j=1:60
for i=1:(2*N)
    for k=(2-mod(i,2)):2:M
    c1r(:,i,j)=c1r(:,i,j)+h(:,:,k)'*c2r_(:,((i-k+4)/2),j);
    end
end
end
c1r=c1;
%一级重构(由c1r重构c0r)
N=length(c1r(:,:,1));
c1r_=[c1r(:,N,:),c1r];
c0r=zeros(2,(2*N),60);
for j=1:60
for i=1:(2*N)
    for k=(2-mod(i,2)):2:M
    c0r(:,i,j)=c0r(:,i,j)+h(:,:,k)'*c1r_(:,((i-k+4)/2),j);
    end
end
end

% 4
%后处理(由c0r恢复原始一维数据)
xr=zeros(2,(length(x)/2),60);
for j=1:60
 c0r_=c0r(:,:,j);   
for k=1:(length(x)/2-1)
  for i=1:m
      xr(:,k,j)=xr(:,k,j)+q(:,:,i)*[c0r_((2*(k+i)-3));c0r_((2*(k+i)-2))];
  end
end
end
 %将复原信号xr转化为一维数据xnr,y应该等于x
 xnr=zeros(60,length(x));
for j=1:60
for i=1:2:(length(x))
    xnr(j,i)=xr(1,((i+1)/2),j);
    xnr(j,(i+1))=xr(2,((i+1)/2),j);
end
end
for i=1:60
xnrti(i,:)=wkeep(xnr(i,:),length(xh));
end
xcal=xnrti((1:40),:);
xtest=xnrti((41:60),:);