hmt.m.txt

小波域的树模型去噪采用ＥＭ算法训练模型参数由Ｍatlab编程实现仿真

clc;
clear;
tic;
wname='haar';  % 采用的小波

[s,noisyx] = wnoise(4,10,5);     %%%%%% 测试信号：bumps 2 'heavy sine' 3,Doppler 4 
e1=noisyx-s;
SNR1=10*log10((s*s')/(e1*e1'))
%%%%%%%    
L0=length(noisyx);       %%% 信号长度
N=5;                     %%% 分解层数    
[c,l] = wavedec(noisyx,N,wname);   % 用wavedec函数对信号进行小波分解
noise = wnoisest(c,l,1:N);        % wnoisest 函数是估计一维小波系数的标准偏差的函数
L=l(2:N+2);               %%%% 各层高频系数的长度和信号长度：从最粗到最细 
D=zeros(N,L0/2);
L1=L(1)+1;
for n=1:N
    L2=L(n)+L1-1;
    x=c(L1:L2);
    D(n,1:L(n))=x;      %%%%% 各层高频系数，从最粗到最细
    L1=L2+1;
end        
%%%%%%%%%%%%%%%%%%%%%   建树   
node.value=0;     %%%%  1节点的值
node.p1=[0.4,0.6];    %%%%  2节点状态概率
node.p2=[0.45,0.1;0.05,0.4];    %%%%  3父子联合概率
node.beta1=[2,7];     %%%%  4以节点为根的小波树分布函数
node.beta2=[3,8];     %%%%  5知道父状态时，节点小波树分布函数 
node.beta3=[4,9];     %%%%  6父树去除节点树后的分布函数
node.alfa=[0.4,0.6];      %%%%  7整棵树去除节点树后分布
node.u=[0,0];         %%%%  8分布均值
node.ro=[noise(1)/2,2*noise(1)];        %%%%  9分布方差 
node.divert=[0.6,0.1;0.4,0.9];  %%%10 状态转移阵
node.likehood=0;

K=L0/(2^N);    %%%% K棵小波树
for n=1:N
    for m=1:2^(n-1)*K
        Tree(n,m)={node};
        Tree{n,m}.value=D(n,m);
        Tree{n,m}.ro=[noise(N+1-n)/2,2*noise(N+1-n)];
    end
end
%%%%%%%%%%%%%%%%%%%%%  EM 算法  
KKK=10;
for testn=1:80000  %%%% while  循环寻优
%%%%%%%%%%%%%%% 向上算法 
for n=1:N
    for m=1:2^(n-1)*K      %%%%  初始化：每个点算一次guass
        u=Tree{n,m}.u;
        ro=Tree{n,m}.ro;
        w=Tree{n,m}.value;
        temp=((w-u).^2)./(2*ro);      
        temp=exp(-temp)./(sqrt(2*pi*ro));         
        Tree{n,m}.beta1=temp*KKK;        
    end
end 

for k=1:K     %%%% K棵小波树
   for n=N:-1:2   %%%  每棵小波树分别使用向上向下算法
     for m=(2^(n-1)*(k-1)+1):(2^(n-1)*k)
        p=floor((m+1)/2);
        Tree{n,m}.beta2=(Tree{n,m}.beta1)*(Tree{n,m}.divert);
        temp=Tree{n-1,p}.beta1;                %%%   取出已算出父节点guass              
        Tree{n-1,p}.beta1=temp.*Tree{n,m}.beta2;
     end
     for m=(2^(n-1)*(k-1)+1):(2^(n-1)*k)
         p=floor((m+1)/2);
         Tree{n,m}.beta3=(Tree{n-1,p}.beta1)./(Tree{n,m}.beta2); %%% 
     end 
   end
   %%%%%%%%%%%%%%  向下算法 
   Tree{1,k}.alfa=Tree{1,k}.p1;       
   for n=2:N
     for m=(2^(n-1)*(k-1)+1):(2^(n-1)*k)
       p=floor((m+1)/2);
       temp=(Tree{n-1,p}.alfa).*(Tree{n,m}.beta3);
       Tree{n,m}.alfa=temp*((Tree{n,m}.divert)');
       Tree{n,m}.likehood=(Tree{n,m}.alfa*Tree{n,m}.beta1');
           %%%% 计算节点p1值
        Tree{n,m}.p1=((Tree{n,m}.alfa).*(Tree{n,m}.beta1))/Tree{n,m}.likehood; 
        TmP=sum(Tree{n,m}.p1);
        Tree{n,m}.p1=Tree{n,m}.p1/TmP;     %%%%%%%  归一化概率值p1 
           %%% 计算节点p2值
        Tree{n,m}.p2=(Tree{n,m}.beta1')*((Tree{n-1,p}.alfa).*(Tree{n,m}.beta3));  
        Tree{n,m}.p2=(Tree{n,m}.p2).*(Tree{n-1,p}.divert);
        Tree{n,m}.p2=Tree{n,m}.p2/Tree{n,m}.likehood;
        TmP=sum(sum(Tree{n,m}.p2));
        Tree{n,m}.p2=Tree{n,m}.p2/TmP;   %%%%%%%%   归一化概率值p2
     end
  end
end %%% 向上向下算法结束
  %%%%%%%%%%%%%%%  对参数更新 
     p1=[0 0];
     p2=[0 0;0 0];
     u=[0,0];
     ro=[0,0];
   for k=1:K     %%% K个根节点
      Tree{1,k}.likehood=(Tree{1,k}.alfa)*(Tree{1,k}.beta1');
      Tree{1,k}.p1=((Tree{1,k}.alfa).*(Tree{1,k}.beta1))/Tree{1,k}.likehood;  %%% 根节点p1
      TmP=sum(Tree{1,k}.p1);
      Tree{1,k}.p1=Tree{1,k}.p1/TmP;  
          %%% 
      p1=Tree{1,k}.p1+p1; 
      u=Tree{1,k}.value*Tree{1,k}.p1+u; 
   end
   for k=1:K
      temp=Tree{1,k}.p1;  %% 保存原p1，以用于更新方差 
      Tree{1,k}.p1=p1/K;  %% 更新p1
      Tree{1,k}.u=u./(K*Tree{1,k}.p1);  %%% 更新均值
      ro=((Tree{1,k}.value-Tree{1,k}.u).^2).*temp+ro;  %%% 更新均值后才算方差  Tree{1,k}.p1
   end
   for k=1:K
      Tree{1,k}.ro=ro./(K*Tree{1,k}.p1);  %%% 更新方差 
   end
  %%%%%%%%%%%%  其余节点
   for n=2:N
      for m=1:2^(n-1)
            p1=[0 0];
            p2=[0 0;0 0];
            u=[0,0];
            ro=[0,0];
          for k=1:K
             r=2^(n-1)*(k-1)+m;
             p1=Tree{n,r}.p1+p1;
             p2=Tree{n,r}.p2+p2;
             u=Tree{n,r}.value*Tree{n,r}.p1+u;
          end
         for k=1:K   %%% 更新
             r=2^(n-1)*(k-1)+m;
             p=floor((r+1)/2);
             temp=Tree{n,r}.p1;     %%% 保存原p1，以用于更新方差
             Tree{n,r}.p1=p1/K;                %%% 更新p1
             Tree{n,r}.divert=p2./(K*[Tree{n-1,p}.p1;Tree{n-1,p}.p1]); %%% 更新转移概率
             Tree{n,r}.u=u./(K*Tree{n,r}.p1);
             ro=((Tree{n,r}.value-Tree{n,r}.u).^2).*temp+ro;  %%% 更新均值后才算方差 Tree{n,r}.p1
         end
         for k=1:K
             r=2^(n-1)*(k-1)+m;
             Tree{n,r}.ro=ro./(K*Tree{n,r}.p1);  %%% 更新方差 
         end
     end
   end
 maxlike=1;
for k=1:K
    maxlike=maxlike*((Tree{1,k}.alfa.*(Tree{1,k}.beta1))*(Tree{1,k}.p1'));     %%%% 小波系数 W 的似然函数
end   %%% 结果
  %if mod(testn,200)==1
     maxlike
     %end
end  %%%%%%%%%%%%%%%%% end_while寻优结束
%%%%%%%%%%%%%%%%%%%%%%%%%  Bayes 估计真实小波系数
for n=1:N
    for m=1:2^(n-1)*K
        ro=noise(1);     %%% 噪声方差
        temp=(abs(Tree{n,m}.ro-ro))./Tree{n,m}.ro;
        value=(Tree{n,m}.p1)*(temp')*Tree{n,m}.value;
        D2(n,m)=value;
    end
end
ED=D-D2;
  %%%%%%%%%%%%%%%%%%%%%%%%%%% 重构
 C=c;
L1=L(1)+1;
for n=1:N
    L2=L(n)+L1-1;
    C(L1:L2)=D2(n,1:L(n));
    L1=L2+1;
end     
Reconstruct= waverec(C,l,wname);  % 用去噪后的小波系数重构信号
e=s-Reconstruct;
SNR=10*log10((s*s')/(e*e'))
figure;
subplot(3,1,1),plot(s),title('原信号');
subplot(3,1,2),plot(noisyx),title('加噪后信号');
subplot(3,1,3),plot(Reconstruct),title('经过处理后的信号');
toc;