开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 其他 > 查看源码
LAWMLShrink2.m
资源名称：小波去噪.rar　[点击查看]
上传用户：sla11nk8
上传日期：2013-03-09
资源大小：21k
文件大小：3k
源码类别：

其他

开发平台：

Matlab

LAWMLShrink2.m：源码内容
  1. function out=LAWMLShrink2(wcY,L)
  2. %
  3. %%%%It's My Own Function!!!!!
  4. %
  5. %%%% out=LAWMLShrink2(wcY,L)
  6. %块大小在不同尺度和不同方向变化
  7. [n,J] = dyadlength(wcY);
  8. ws=wcY;
  9. for j=(J-1):-1:L
  10. scale=ScaleNumber(j,n);
  11. for kk=1:3%对3个方向的系数分别处理
  12.      switch kk
  13.         case 1
  14.             %得到当前矩阵
  15.             [t1,t2]=dyad2LH(j);
  16.             wcCorrent=wcY(t1,t2);
  17.             if j~=L
  18.             %得到其父结点矩阵
  19.                [p1,p2]=dyad2LH(j-1);
  20.                wcParent=wcY(p1,p2);
  21.             end;
  22.             if scale>=3
  23.                 linLength=3;
  24.             else
  25.                 linLength=-2*scale+9;
  26.             end;
  27.             
  28.         case 2
  29.             %得到当前矩阵
  30.             [t1,t2]=dyad2HL(j);
  31.             wcCorrent=wcY(t1,t2);
  32.             if j~=L
  33.             %得到其父结点矩阵
  34.                [p1,p2]=dyad2HL(j-1);
  35.                wcParent=wcY(p1,p2);
  36.             end;
  37.             if scale>=3
  38.                 linLength=3;
  39.             else
  40.                 linLength=-2*scale+9;
  41.             end;
  42.             
  43.         case 3
  44.             %得到当前矩阵
  45.             [t1,t2]=dyad2HH(j);
  46.             wcCorrent=wcY(t1,t2);
  47.             if j~=L
  48.             %得到其父结点矩阵
  49.                [p1,p2]=dyad2HH(j-1);
  50.                wcParent=wcY(p1,p2);
  51.             end;
  52.             if scale>3
  53.                 linLength=3;
  54.             else
  55.                 linLength=-2*scale+11;
  56.             end;
  57.             
  58.        end;
  59.          
  60.      nn=2^j;     
  61.      %linLength=5;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  62.      bb=floor(linLength/2);
  63.      %数据延拓
  64.      wcExtent(bb+1:nn+bb,bb+1:nn+bb)=wcCorrent;
  65.      wcExtent(1:bb,bb+1:nn+bb)=wcCorrent(bb:(-1):1,1:nn);%上行
  66.      wcExtent(bb+1:nn+bb,1:bb)=wcCorrent(1:nn,bb:(-1):1);%左列
  67.      wcExtent(nn+bb+1:nn+2*bb,bb+1:nn+bb)=wcCorrent(nn:(-1):nn-bb+1,1:nn);%下行
  68.      wcExtent(bb+1:nn+bb,nn+bb+1:nn+2*bb)=wcCorrent(1:nn,nn:(-1):nn-bb+1);%右列
  69.      wcExtent(1:bb,1:bb)=wcCorrent(bb:(-1):1,bb:(-1):1);%左上角
  70.      wcExtent(1:bb,nn+bb+1:nn+2*bb)=wcCorrent(bb:(-1):1,nn:(-1):nn-bb+1);%右上角
  71.      wcExtent(nn+bb+1:nn+2*bb,1:bb)=wcCorrent(nn:(-1):nn-bb+1,bb:(-1):1);%左下角
  72.      wcExtent(nn+bb+1:nn+2*bb,nn+bb+1:nn+2*bb)=wcCorrent(nn:(-1):nn-bb+1,nn:(-1):nn-bb+1);%右下角
  73.      
  74.      z=zeros(nn+bb*2,nn+bb*2);
  75.      %用领域的点估计Y的方差
  76.      for k1=bb+1:nn+bb
  77.          for k2=bb+1:nn+bb
  78.              u=wcExtent(k1-bb:k1+bb,k2-bb:k2+bb);
  79.              sumY2=sum(sum(u.^2));
  80.              %噪声方差归一化
  81.              sigma=1;
  82.              sitaX2=sumY2/linLength/linLength/1-sigma^2;
  83.              if (sitaX2<0)%噪声太大
  84.                 z(k1,k2)=0;
  85.              else
  86.                 z(k1,k2)=sitaX2*wcExtent(k1,k2)/(sitaX2+sigma^2);
  87.                 
  88.              end;
  89.          end;
  90.      end;
  91.      ws(t1,t2)=z(bb+1:nn+bb,bb+1:nn+bb);
  92. end;
  93. end;
  94. out=ws;