adaptshrink2.m

上課所教所有matlab 小波去燥的教程！~非常適合各位的參考！非常好喔！

function out=AdaptShrink2(wcY,L)
%
%%%%It's My Own Function!!!!!
%
%%%% out=AdaptShrink(wcY,L)
[n,J] = dyadlength(wcY);
ws=wcY;
for j=(J-1):-1:L
 for kk=1:3%对3个方向的系数分别处理
     switch kk
        case 1
            %得到当前矩阵
            [t1,t2]=dyad2LH(j);
            wcCorrent=wcY(t1,t2);
            if j~=L
            %得到其父结点矩阵
               [p1,p2]=dyad2LH(j-1);
               wcParent=wcY(p1,p2);
            end;
            
        case 2
            %得到当前矩阵
            [t1,t2]=dyad2HL(j);
            wcCorrent=wcY(t1,t2);
            if j~=L
            %得到其父结点矩阵
               [p1,p2]=dyad2HL(j-1);
               wcParent=wcY(p1,p2);
            end;
        case 3
            %得到当前矩阵
            [t1,t2]=dyad2HH(j);
            wcCorrent=wcY(t1,t2);
            if j~=L
            %得到其父结点矩阵
               [p1,p2]=dyad2HH(j-1);
               wcParent=wcY(p1,p2);
            end;
            
       end;
         
    
     %求出z(k1,k2)
     z=abs(wcCorrent);
     nn=2^j;
     for k1=2:(nn-1)
         for k2=2:(nn-1)
             %求向量u
             u(1)=wcCorrent(k1-1,k2-1);
             u(2)=wcCorrent(k1,k2-1);
             u(3)=wcCorrent(k1+1,k2-1);
             u(4)=wcCorrent(k1+1,k2);
             u(5)=wcCorrent(k1+1,k2+1);
             u(6)=wcCorrent(k1,k2+1);
             u(7)=wcCorrent(k1-1,k2+1);
             u(8)=wcCorrent(k1-1,k2);
             u(9)=wcCorrent(k1,k2);
             if (j==L)
                 u(10)=u(9);
             else
                 u(10)=wcParent(round(k1/2),round(k2/2));
             end;
             %求向量w
             %w=ones(1,10)/10;
             w=[1 1 1 1 1 1 1 1 0 1]./9;
             z(k1,k2)=sum(w.*abs(u));%求出z(k1,k2)
         end;
     end;
        
     %对z矩阵的元素排序
     z_1d=[];
     for k1=1:nn
         z_1d=[z_1d z(k1,:)];
     end;
     [zOrder,index]=sort(z_1d);
     %计算窗口大小LL
     LL=max(50,floor(0.02*nn*nn));
     
     [k1,k2]=MapTo2d(index,nn);%1维向量的标号映射到2维标号
     
     z=wcCorrent;
     for p=1:nn*nn
         %计算实际窗口大小
         begin=max(p-LL,1);
         eend=min(p+LL,nn*nn);
         total=eend-begin+1;
              
      if (p==1)
             sumY2=0;
             for t=1:(LL+1)
                 sumY2=sumY2+wcCorrent(k1(t),k2(t))^2;
             end;
         elseif (p-LL)<2%窗口左侧数据不够
                 sumY2=sumY2+wcCorrent(k1(p+LL),k2(p+LL))^2;
             elseif (p+LL)>nn*nn%窗口右侧数据不够
                     tempp=wcCorrent(k1(p-LL-1),k2(p-LL-1))^2;
                     sumY2=sumY2-tempp;
                 else%窗口数据足够
                     tempp=wcCorrent(k1(p-LL-1),k2(p-LL-1))^2;
                     sumY2=sumY2-tempp;
                     tempp=wcCorrent(k1(p+LL),k2(p+LL))^2;
                     sumY2=sumY2+tempp;
         end;
      
        %噪声方差归一化
         sigma=1;
         sitaX2=sumY2/total-sigma^2;
         if (sitaX2<0)%噪声太大
             z(k1(p),k2(p))=0;
         else
             thresh=sigma^2/sqrt(sitaX2);
             z(k1(p),k2(p))=SoftThresh(wcCorrent(k1(p),k2(p)),thresh);
         end;
     end;
     ws(t1,t2)=z;
 end;
end;
out=ws;      
     
     function [k1,k2]=MapTo2d(k,n)
     %1维向量的标号映射到2维标号
     k1=ceil(k./n);
     k2=k-(k1-1).*n;