lawmlshrink2.m

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

function out=LAWMLShrink2(wcY,L)
%
%%%%It's My Own Function!!!!!
%
%%%% out=LAWMLShrink2(wcY,L)
%块大小在不同尺度和不同方向变化
[n,J] = dyadlength(wcY);
ws=wcY;
for j=(J-1):-1:L
scale=ScaleNumber(j,n);
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;
            if scale>=3
                linLength=3;
            else
                linLength=-2*scale+9;
            end;
            
        case 2
            %得到当前矩阵
            [t1,t2]=dyad2HL(j);
            wcCorrent=wcY(t1,t2);
            if j~=L
            %得到其父结点矩阵
               [p1,p2]=dyad2HL(j-1);
               wcParent=wcY(p1,p2);
            end;
            if scale>=3
                linLength=3;
            else
                linLength=-2*scale+9;
            end;
            
        case 3
            %得到当前矩阵
            [t1,t2]=dyad2HH(j);
            wcCorrent=wcY(t1,t2);
            if j~=L
            %得到其父结点矩阵
               [p1,p2]=dyad2HH(j-1);
               wcParent=wcY(p1,p2);
            end;
            if scale>3
                linLength=3;
            else
                linLength=-2*scale+11;
            end;
            
       end;
         
     nn=2^j;     
     %linLength=5;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     bb=floor(linLength/2);
     %数据延拓
     wcExtent(bb+1:nn+bb,bb+1:nn+bb)=wcCorrent;
     wcExtent(1:bb,bb+1:nn+bb)=wcCorrent(bb:(-1):1,1:nn);%上行
     wcExtent(bb+1:nn+bb,1:bb)=wcCorrent(1:nn,bb:(-1):1);%左列
     wcExtent(nn+bb+1:nn+2*bb,bb+1:nn+bb)=wcCorrent(nn:(-1):nn-bb+1,1:nn);%下行
     wcExtent(bb+1:nn+bb,nn+bb+1:nn+2*bb)=wcCorrent(1:nn,nn:(-1):nn-bb+1);%右列
     wcExtent(1:bb,1:bb)=wcCorrent(bb:(-1):1,bb:(-1):1);%左上角
     wcExtent(1:bb,nn+bb+1:nn+2*bb)=wcCorrent(bb:(-1):1,nn:(-1):nn-bb+1);%右上角
     wcExtent(nn+bb+1:nn+2*bb,1:bb)=wcCorrent(nn:(-1):nn-bb+1,bb:(-1):1);%左下角
     wcExtent(nn+bb+1:nn+2*bb,nn+bb+1:nn+2*bb)=wcCorrent(nn:(-1):nn-bb+1,nn:(-1):nn-bb+1);%右下角
     
     z=zeros(nn+bb*2,nn+bb*2);
     %用领域的点估计Y的方差
     for k1=bb+1:nn+bb
         for k2=bb+1:nn+bb
             u=wcExtent(k1-bb:k1+bb,k2-bb:k2+bb);
             sumY2=sum(sum(u.^2));
             %噪声方差归一化
             sigma=1;
             sitaX2=sumY2/linLength/linLength/1-sigma^2;
             if (sitaX2<0)%噪声太大
                z(k1,k2)=0;
             else
                z(k1,k2)=sitaX2*wcExtent(k1,k2)/(sitaX2+sigma^2);
                
             end;
         end;
     end;
     ws(t1,t2)=z(bb+1:nn+bb,bb+1:nn+bb);
end;
end;
out=ws;