声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 3439|回复: 10

[共享资源] 提升小波变换

[复制链接]
发表于 2006-12-14 21:53 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
简单实现了一下提升小波变换。包括Haar,(1,1),5/3,9/7 提升结构。
欢迎有兴趣的同学一起讨论

  1. function [res]=lifting_transform_forward(im,levels,wavetype)
  2. wavelet = wavetype;
  3. for k =1:levels
  4.     [L,H]=Nonlinear_lifting_forward(im,wavelet);
  5.     [LL,LH{k}] = Nonlinear_lifting_forward(L',wavelet);
  6.     [HL{k},HH{k}] = Nonlinear_lifting_forward(H',wavelet);
  7.     im = LL';
  8. end

  9. for k= levels:-1:1
  10.     temp1 =cat(2,LL,LH{k});
  11.     temp2 =cat(2,HL{k},HH{k});
  12.     lifting_im = cat(1,temp1,temp2);
  13.     LL = lifting_im;
  14. end
  15. res = lifting_im';
复制代码

  1. function [res]=lifting_transform_inverse(im,levels,wavetype)
  2. wavelet = wavetype;
  3. [m,n] = size(im);
  4. min_len = m/2^levels;
  5. LL = im([1:min_len],[1:min_len]);
  6. for k= levels:-1:1
  7.     im_len = m/2^k;
  8.     LH{k} = im([1:im_len],[im_len+1:im_len*2]);
  9.     HL{k} = im([im_len+1:im_len*2],[1:im_len]);
  10.     HH{k} = im([im_len+1:im_len*2],[im_len+1:im_len*2]);
  11. end
  12. for k =levels:-1:1
  13.     up = Nonlinear_lifting_inverse(LL,LH{k},wavelet);
  14.     down = Nonlinear_lifting_inverse(HL{k},HH{k},wavelet);
  15.     im  = Nonlinear_lifting_inverse(up',down',wavelet);
  16.     LL= im';
  17. end
  18. res = im';
复制代码


  1. function [c,d]=Nonlinear_lifting_forward(im,wavelet)
  2. [width,height] = size(im);
  3. %%%%split
  4. for i=1:2:(width-1)
  5.     Xo((i+1)/2,:) = im(i,:);
  6. end
  7. for i=2:2:width
  8.     Xe(i/2,:) = im(i,:);
  9. end
  10. [sub_im_height,sub_im_width]=size(Xe);
  11. %%%%   5/3 lifting %%%%%%%%%%%%%%
  12. if isequal(wavelet, '5.3')
  13.     next_Xe([1:sub_im_height-1],:) = Xe([2:sub_im_height],:);
  14.     next_Xe(sub_im_height,:) = Xe(sub_im_height,:);
  15.     d(:,:) = Xo(:,:) - (0.5*Xe(:,:)+0.5*next_Xe(:,:));
  16.     up_d([2:sub_im_height],:) = d([1:sub_im_height-1],:);
  17.     up_d(1,:) = d(1,:);
  18.     c(:,:) = Xe(:,:) +(0.25*(up_d(:,:)+d(:,:)));
  19. end
  20. %%%%   9/7 lifting %%%%%%%%%%%%%%
  21. if isequal(wavelet, '9.7')
  22.     alfa    = -1.586134342;
  23.     beta    = -0.052980118;
  24.     gamma   =  0.882911075;
  25.     delta   =  0.443506852;
  26.     k       =  1.149604398;
  27.     C0 = Xe;
  28.     D0 = Xo;
  29.     C0_1([1:sub_im_height-1],:) = C0([2:sub_im_height],:);
  30.     C0_1(sub_im_height,:) = C0(sub_im_height,:);
  31.     D1(:,:) = D0(:,:)+alfa.*(C0(:,:)+C0_1(:,:));
  32.     D1_1([2:sub_im_height],:) = D1([1:sub_im_height-1],:);
  33.     D1_1(1,:) = D1(1,:);
  34.     C1(:,:) = C0(:,:)+beta.*(D1(:,:)+D1_1(:,:));
  35.     C1_1([1:sub_im_height-1],:) = C1([2:sub_im_height],:);
  36.     C1_1(sub_im_height,:) = C1(sub_im_height,:);
  37.     D2(:,:) = D1(:,:)+gamma.*(C1(:,:)+C1_1(:,:));
  38.     D2_1([2:sub_im_height],:) = D2([1:sub_im_height-1],:);
  39.     D2_1(1,:) = D2(1,:);
  40.     C2(:,:) = C1(:,:)+delta.*(D2(:,:)+D2_1(:,:));
  41.     c = k.* C2;
  42.     d = D2./k;
  43. end
  44. if  isequal(wavelet, '1.1')
  45.     c(:,:) = 0.5*(Xe(:,:)+Xo(:,:));
  46.     d(:,:) = Xo(:,:)-c(:,:);
  47. end

  48. if isequal(wavelet,'Haar')
  49.     c0 = Xe;
  50.     d0 = Xo;
  51.     d(:,:) = d0(:,:)-c0(:,:);
  52.     c(:,:) = c0(:,:)+0.5.*d(:,:);
  53. end
  54. % figure;
  55. % imshow(uint8(c));
  56. % figure;
  57. % imshow(uint8(d));
  58. %
复制代码


  1. function [res]=Nonlinear_lifting_inverse(L,H,wavelet)
  2. c = L;
  3. d = H;
  4. [sub_im_height,sub_im_width]=size(L);
  5. %%%%   5/3 lifting %%%%%%%%%%%%%%
  6. if isequal(wavelet, '5.3')
  7.     up_d([2:sub_im_height],:) = d([1:sub_im_height-1],:);
  8.     up_d(1,:) = d(1,:);
  9.     Xe(:,:) = c(:,:) - (0.25*(up_d(:,:)+d(:,:)));
  10.    
  11.     next_Xe([1:sub_im_height-1],:) = Xe([2:sub_im_height],:);
  12.     next_Xe(sub_im_height,:) = Xe(sub_im_height,:);
  13.     Xo(:,:) = d(:,:) + (0.5*Xe(:,:)+0.5*next_Xe(:,:));
  14. end
  15. %%%%   9/7 lifting %%%%%%%%%%%%%%
  16. if isequal(wavelet, '9.7')
  17.     alfa    = -1.586134342;
  18.     beta    = -0.052980118;
  19.     gamma   =  0.882911075;
  20.     delta   =  0.443506852;
  21.     k       =  1.149604398;
  22.     D2 = d.*k;
  23.     C2 = c./k;
  24.     D2_1([2:sub_im_height],:) = D2([1:sub_im_height-1],:);
  25.     D2_1(1,:) = D2(1,:);
  26.     C1(:,:) = C2(:,:)-delta.*(D2(:,:)+D2_1(:,:));
  27.     C1_1([1:sub_im_height-1],:) = C1([2:sub_im_height],:);
  28.     C1_1(sub_im_height,:) = C1(sub_im_height,:);
  29.     D1(:,:) = D2(:,:)-gamma.*(C1(:,:)+C1_1(:,:));
  30.     D1_1([2:sub_im_height],:) = D1([1:sub_im_height-1],:);
  31.     D1_1(1,:) = D1(1,:);
  32.     C0(:,:) = C1(:,:)-beta.*(D1(:,:)+D1_1(:,:));
  33.     C0_1([1:sub_im_height-1],:) = C0([2:sub_im_height],:);
  34.     C0_1(sub_im_height,:) = C0(sub_im_height,:);
  35.     D0(:,:) = D1(:,:)-alfa.*(C0(:,:)+C0_1(:,:));
  36.      Xe = C0;
  37.      Xo = D0;
  38. end
  39. if  isequal(wavelet, '1.1')
  40.     d0(:,:) = d(:,:)+c(:,:);
  41.     c0(:,:) = 2.*(c(:,:)-0.5.*d0(:,:));
  42.      Xe = c0;
  43.      Xo = d0;
  44. end
  45. if isequal(wavelet,'Haar')
  46.     c0(:,:) = c(:,:)-0.5.*d(:,:);
  47.     d0(:,:) = d(:,:)+c0(:,:);
  48.     Xe = c0;
  49.     Xo = d0;
  50. end
  51. res([1:2:sub_im_height*2-1],:) = Xo(:,:);
  52. res([2:2:sub_im_height*2],:) = Xe(:,:);
复制代码

  1. function testing_lifting_linjie
  2. close all;
  3. clear all;
  4. clc;

  5. im = double(imread('lena.bmp'));
  6. figure;
  7. subplot(221);
  8. imshow(uint8(im));
  9. title('Original Image');  

  10. levels =3;
  11. [im_h,im_w] = size(im);
  12. lifting_im = lifting_transform_forward(im,levels,'5.3');
  13. subplot(222);
  14. imshow(uint8(lifting_im));
  15. title('lifting coefficients');  

  16. result_im = lifting_transform_inverse(lifting_im,levels,'5.3');
  17. high =  sum(sum(abs(lifting_im(:,:))))
  18. MSE=sum(sum((result_im-im).^2))/(im_h*im_w);
  19. PSNR = 20*log((255*255)/MSE);
  20. disp(['     PSNR =  ',num2str(PSNR)]);
  21. subplot(223);
  22. imshow(uint8(result_im));
  23. title('reconstruction image')
复制代码

评分

1

查看全部评分

回复
分享到:

使用道具 举报

发表于 2007-2-2 23:39 | 显示全部楼层
顶啊!!
发表于 2007-2-5 10:38 | 显示全部楼层

很好的东西!

很好的东西!
希望能学习学习!
西安电子科技大学的同胞!
我是西安的,可以联系吗?
发表于 2008-2-21 17:18 | 显示全部楼层
谢谢 学到东西了
发表于 2008-2-25 04:59 | 显示全部楼层
谢谢,赞~
发表于 2008-11-26 16:36 | 显示全部楼层
好东西,谢谢版主,辛苦了
发表于 2008-11-27 11:52 | 显示全部楼层

:handshake ,:victory:
发表于 2009-1-29 05:16 | 显示全部楼层
我下载了代码.
在使用过程中我发现了一个bug.
由于我只需要5.3提升小波,其他的就没去检查了.
在5.3提升小波的实现中,以下代码应该存在bug:

for k= levels:-1:1
    temp1 =cat(2,LL,LH{k});
    temp2 =cat(2,HL{k},HH{k});
    lifting_im = cat(1,temp1,temp2);
    LL = lifting_im;
end
res = im';


for k =levels:-1:1
    up = Nonlinear_lifting_inverse(LL,LH{k},wavelet);
    down = Nonlinear_lifting_inverse(HL{k},HH{k},wavelet);
    im  = Nonlinear_lifting_inverse(up',down',wavelet);
    LL= im';
end

应该改为:
for k= levels:-1:1
    temp1 =cat(1,LL,LH{k});
    temp2 =cat(1,HL{k},HH{k});
    lifting_im = cat(1, temp1' , temp2' );
    LL = lifting_im';
end
res = lifting_im;


for k =levels:-1:1
    up = Nonlinear_lifting_inverse(LL' , LH{k}' , wavelet);
    down = Nonlinear_lifting_inverse(HL{k}' , HH{k}' , wavelet);
    im  = Nonlinear_lifting_inverse(up', down', wavelet);
    LL= im;
end
res = im;

我还修改了函数 Nonlinear_lifting_forward 和 Nonlinear_lifting_inverse 中 5.3提升小波的计算方法,主要是循环扩展,和 Y(2n) = X(2n) + [(Y(2n-1) + Y(2n+1) + 2)/4]; (原文使用的公式为:Y(2n) = X(2n) + [(Y(2n-1) + Y(2n+1))/4])

[ 本帖最后由 zhu1982lin 于 2009-1-29 05:30 编辑 ]
发表于 2009-3-24 20:18 | 显示全部楼层

回复 8楼 zhu1982lin 的帖子

我按照你的改了 可还是不行啊!一直提示有错!
??? Error using ==> cat
CAT arguments dimensions are not consistent.

Error in ==> lifting_transform_forward at 13
    lifting_im = cat(1,temp1',temp2');
发表于 2012-9-12 15:18 | 显示全部楼层
亲有木有提升小波包的代码哇
发表于 2012-9-12 18:58 | 显示全部楼层

  1. function varargout = lwt(x,LS,varargin)
  2. %LWT Lifting wavelet decomposition 1-D.
  3. %   LWT performs a 1-D lifting wavelet decomposition
  4. %   with respect to a particular lifted wavelet that you specify.
  5. %
  6. %   [CA,CD] = LWT(X,W) computes the approximation
  7. %   coefficients vector CA and detail coefficients vector CD,
  8. %   obtained by a lifting wavelet decomposition, of
  9. %   the vector X. W is a lifted wavelet name (see LIFTWAVE).
  10. %
  11. %   X_InPlace = LWT(X,W) computes the approximation and
  12. %   detail coefficients. These coefficients are stored in-place:
  13. %     CA = X_InPlace(1:2:end) and CD = X_InPlace(2:2:end)
  14. %
  15. %   LWT(X,W,LEVEL) computes the lifting wavelet decomposition
  16. %   at level LEVEL.
  17. %
  18. %   X_InPlace = LWT(X,W,LEVEL,'typeDEC',typeDEC) or
  19. %   [CA,CD] = LWT(X,W,LEVEL,'typeDEC',typeDEC) with
  20. %   typeDEC = 'w' or 'wp' computes the wavelet or the
  21. %   wavelet packet decomposition using lifting, at level LEVEL.
  22. %
  23. %   Instead of a lifted wavelet name, you may use the associated
  24. %   lifting scheme LS:
  25. %     LWT(X,LS,...) instead of LWT(X,W,...).
  26. %
  27. %   For more information about lifting schemes type: lsinfo.
  28. %
  29. %   See also ILWT.

  30. %   M. Misiti, Y. Misiti, G. Oppenheim, J.M. Poggi 02-Feb-2000.
  31. %   Last Revision: 10-Jul-2003.
  32. %   Copyright 1995-2004 The MathWorks, Inc.
  33. %   $Revision: 1.1.6.3 $  $Date: 2004/04/13 00:39:57 $

  34. % Check arguments.
  35. nbIn = nargin;
  36. switch nbIn
  37.     case {2,3,5} ,
  38.     case {0,1} , error('Not enough input arguments.');
  39.     case {4}   , error('Invalid number of input arguments.');
  40.     otherwise  , error('Too many input arguments.');
  41. end
  42. if nargout>3
  43.     error('Too many output arguments.');
  44. end

  45. % Default: level and typeDEC.
  46. level = 1;
  47. typeDEC = 'w';
  48. if nargin>2
  49.     level = varargin{1};
  50.     for k = 2:2:length(varargin)-1
  51.         argName = lower( varargin{k});
  52.         switch argName
  53.             case 'typedec' , typeDEC = varargin{k+1};
  54.         end
  55.     end
  56. end
  57. % if level>1
  58. %     msg = nargoutchk(0,1,nargout); error(msg);
  59. % end
  60. if ischar(LS) , LS = liftwave(LS); end

  61. %===================%
  62. % LIFTING ALGORITHM %
  63. %===================%
  64. % Splitting.
  65. lx = length(x);
  66. firstIdxAPP = 1;
  67. firstIdxDET = 1+mod(firstIdxAPP,2);
  68. idxAPP = firstIdxAPP:2:lx;
  69. idxDET = firstIdxDET:2:lx;
  70. lenAPP = length(idxAPP);
  71. lenDET = length(idxDET);

  72. % Lifting.
  73. NBL = size(LS,1);
  74. LStype = LS{NBL,3};
  75. for k = 1:NBL-1
  76.     liftTYPE = LS{k,1};
  77.     liftFILT = LS{k,2};
  78.     DF       = LS{k,3};
  79.     switch liftTYPE
  80.        case 'p' ,
  81.            x(idxAPP) = x(idxAPP) + ...
  82.                lsupdate('v',x(idxDET),liftFILT,DF,lenAPP,LStype);
  83.        case 'd' ,
  84.            x(idxDET) = x(idxDET) + ...
  85.                lsupdate('v',x(idxAPP),liftFILT,DF,lenDET,LStype);
  86.     end
  87. end

  88. % Normalization.
  89. if isempty(LStype)
  90.     x(idxAPP) = LS{NBL,1}*x(idxAPP);
  91.     x(idxDET) = LS{NBL,2}*x(idxDET);
  92. end
  93. %========================================================================%

  94. % Recursion if level > 1.
  95. if level>1
  96.    x(idxAPP) = lwt(x(idxAPP),LS,level-1,'typeDEC',typeDEC);
  97.    if isequal(typeDEC,'wp')
  98.        x(idxDET) = lwt(x(idxDET),LS,level-1,'typeDEC',typeDEC);
  99.    end
  100. end

  101. switch nargout
  102.   case 1 , varargout = {x};
  103.   case 2 , varargout = {x(idxAPP),x(idxDET)};
  104.   case 3 , varargout = {x,x(idxAPP),x(idxDET)};
  105. end
复制代码

评分

1

查看全部评分

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-13 13:03 , Processed in 0.091802 second(s), 19 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表