wextend.m
上传用户:haiyisale
上传日期:2013-01-09
资源大小:3246k
文件大小:13k
源码类别:

波变换

开发平台:

Matlab

  1. function x = wextend(type,mode,x,lf,location)
  2. %WEXTEND Extend a Vector or a Matrix.
  3. %
  4. %   Y = WEXTEND(TYPE,MODE,X,L,LOC) or
  5. %   Y = WEXTEND(TYPE,MODE,X,L)
  6. %
  7. %   The valid extension types (TYPE) are:
  8. %     1,'1','1d' or '1D'    : 1-D extension
  9. %     2,'2','2d' or '2D'    : 2-D extension 
  10. %     'ar' or 'addrow'      : add rows  
  11. %     'ac' or 'addcol       : add columns   
  12. %
  13. %   The valid extension modes (MODE) are:
  14. %     'zpd' zero extension.
  15. %     'sp0' smooth extension of order 0.
  16. %     'spd' (or 'sp1') smooth extension of order 1.
  17. %     'sym' (or 'symh') symmetric extension (half-point).
  18. %     'symw' symmetric extension (whole-point).
  19. %     'asym' (or 'asymh') antisymmetric extension (half-point).
  20. %     'asymw' antisymmetric extension (whole-point).
  21. %     'ppd' periodized extension (1).
  22. %     'per' periodized extension (2):
  23. %        If the signal length is odd, WEXTEND adds an extra-sample
  24. %        equal to the last value on the right and performs extension
  25. %        using the 'ppd' mode. Otherwise, 'per' reduces to 'ppd'.
  26. %        The same kind of rule stands for images.
  27. %
  28. %   With TYPE = {1,'1','1d' or '1D'}: 
  29. %     LOC = 'l' (or 'u') for left (up) extension.
  30. %     LOC = 'r' (or 'd') for right (down) extension.
  31. %     LOC = 'b' for extension on both sides.
  32. %     LOC = 'n' nul extension
  33. %     The default is LOC = 'b'.
  34. %     L is the length of the extension.
  35. %
  36. %   With TYPE = {'ar','addrow'}
  37. %     LOC is a 1D extension location.
  38. %     The default is LOC = 'b'.
  39. %     L is the number of rows to add.
  40. %
  41. %   With TYPE = {'ac','addcol'}
  42. %     LOC is a 1D extension location.
  43. %     The default is LOC = 'b'.
  44. %     L is the number of columns to add.
  45. %
  46. %   With TYPE = {2,'2','2d' or '2D'}:
  47. %     LOC = [locrow,loccol] where locrow and loccol are 1D
  48. %     extension locations or 'n' (none).
  49. %     The default is LOC = 'bb'.
  50. %     L = [lrow,lcol] where lrow is the number of rows 
  51. %     to add and lcol is the number of columns to add.
  52. %   M. Misiti, Y. Misiti, G. Oppenheim, J.M. Poggi 15-Nov-97.
  53. %   Last Revision: 14-Jul-2003.
  54. %   Copyright 1995-2004 The MathWorks, Inc.
  55. %   $Revision: 1.10.4.2 $  $Date: 2004/03/15 22:42:42 $
  56. type = lower(type);
  57. switch type
  58.   case {1,'1','1d'}
  59.       if nargin<5 , loc = 'b'; else , loc = testLoc(location); end
  60.       if isequal(loc,'n') , return; end
  61.       isROW = (size(x,1)<2);
  62.       x  = x(:)';
  63.       sx = length(x);
  64.       switch mode
  65.         case 'zpd'            % Zero Padding.
  66.           [ext_L,ext_R] = getZPD_ext(1,lf,loc);
  67.           x = [ext_L x ext_R];
  68.         case {'sym','symh'}   % Half-point Symmetrization .
  69.           I = getSymIndices(sx,lf,loc);
  70.           x  = x(I);
  71.         case {'symw'}         % Whole-point Symmetrization.
  72.           x = WP_SymExt(x,lf,loc);
  73.           
  74.         case {'asym','asymh'} % Half-point Anti-Symmetrization.
  75.           x = HP_AntiSymExt(x,lf,loc);  
  76.         case {'asymw'}        % Whole-point Anti-Symmetrization.
  77.           x = WP_AntiSymExt(x,lf,loc);  
  78.         case 'sp0'            % Smooth padding of order 0.
  79.           [ext_L,ext_R] = getSP0_ext('row',x(1),x(sx),lf,loc);
  80.           x = [ext_L x ext_R];
  81.         case {'spd','sp1'}    % Smooth padding of order 1.
  82.           if sx<2 , d = 0; else , d = 1; end
  83.           d_L = x(1)-x(1+d);
  84.           d_R = x(sx)-x(sx-d);
  85.           [ext_L,ext_R] = getSP1_ext('row',x(1),x(sx),d_L,d_R,lf,loc);                  
  86.           x = [ext_L x ext_R];
  87.         case {'ppd'}          % Periodization.
  88.           I = getPerIndices(sx,lf,loc);
  89.           x = x(I);
  90.         case {'per'}          % Periodization.
  91.           if rem(sx,2) , x(sx+1) = x(sx); sx = sx+1; end
  92.           I = getPerIndices(sx,lf,loc);
  93.           x = x(I);
  94.         otherwise
  95.           errargt(mfilename,'Invalid Extension Mode for DWT!','msg');
  96.           error('*')
  97.       end
  98.       if ~isROW , x = x'; end
  99.   case {2,'2','2d'}
  100.       if nargin<5
  101.           locRow = 'b'; locCol = 'b';
  102.       else
  103.           if length(location)<2 , location(2) = location(1); end
  104.           locRow = testLoc(location(1)); 
  105.           locCol = testLoc(location(2));
  106.       end
  107.       if length(lf)<2 , lf = [lf , lf]; end
  108.       [rx,cx] = size(x); 
  109.       switch mode
  110.         case 'zpd'            % Zero Padding.
  111.             if ~isequal(locCol,'n') ,
  112.                 [ext_L,ext_R] = getZPD_ext(rx,lf(2),locCol);
  113.                 x  = [ext_L x ext_R];
  114.             end
  115.             if ~isequal(locRow,'n') ,
  116.                 cx = size(x,2);
  117.                 [ext_L,ext_R] = getZPD_ext(lf(1),cx,locRow);
  118.                 x = [ext_L ; x ; ext_R];
  119.             end
  120.         case {'sym','symh'}   % Symmetrization half-point.
  121.             if ~isequal(locCol,'n') ,
  122.                 I = getSymIndices(cx,lf(2),locCol); x = x(:,I);
  123.             end
  124.             if ~isequal(locRow,'n') ,
  125.                 I = getSymIndices(rx,lf(1),locRow); x = x(I,:);
  126.             end
  127.         case {'symw'}         % Symmetrization whole-point.
  128.             if ~isequal(locCol,'n') ,
  129.                 x = WP_SymExt(x,lf(2),locCol);
  130.             end
  131.             if ~isequal(locRow,'n') ,
  132.                 x = WP_SymExt(x',lf(1),locRow); x = x';
  133.             end
  134.           
  135.         case {'asym','asymh'} % Half-point Anti-Symmetrization.
  136.           if ~isequal(locCol,'n') , 
  137.               x = HP_AntiSymExt(x,lf(2),locCol);
  138.           end
  139.           if ~isequal(locRow,'n') , 
  140.               x = HP_AntiSymExt(x',lf(1),locRow); x = x';
  141.           end
  142.         case {'asymw'}        % Whole-point Anti-Symmetrization.
  143.             if ~isequal(locCol,'n') , 
  144.                 x = WP_AntiSymExt(x,lf(2),locCol);
  145.             end
  146.             if ~isequal(locRow,'n') , 
  147.                 x = WP_AntiSymExt(x',lf(1),locRow); x = x';
  148.             end
  149.         case 'sp0'            % Smooth padding of order 0.
  150.             if ~isequal(locCol,'n') ,   
  151.                 [ext_L,ext_R] = getSP0_ext('row',x(:,1),x(:,cx),lf(2),locCol);
  152.                 x = [ext_L x ext_R];
  153.             end
  154.             if ~isequal(locRow,'n') , 
  155.                 [ext_L,ext_R] = getSP0_ext('col',x(1,:),x(rx,:),lf(1),locRow);
  156.                 x = [ext_L ; x ; ext_R];
  157.             end
  158.         case {'spd','sp1'}    % Smooth padding of order 1.
  159.             if ~isequal(locCol,'n') ,   
  160.                 if cx<2 , d = 0; else , d = 1; end
  161.                 d_L = x(:,1) - x(:,1+d);
  162.                 d_R = x(:,cx)- x(:,cx-d);
  163.                 [ext_L,ext_R] = getSP1_ext('row',x(:,1),x(:,cx),d_L,d_R,lf(2),locCol);
  164.                 x = [ext_L x ext_R];
  165.             end
  166.             if ~isequal(locRow,'n') , 
  167.                 if rx<2 , d = 0; else , d = 1; end          
  168.                 d_L = x(1,:) - x(1+d,:);
  169.                 d_R = x(rx,:)- x(rx-d,:);        
  170.                 [ext_L,ext_R] = getSP1_ext('col',x(1,:),x(rx,:),d_L,d_R,lf(1),locRow);
  171.                 x = [ext_L ; x ; ext_R];
  172.             end
  173.         case 'ppd'            % Periodization.
  174.             if ~isequal(locCol,'n') ,   
  175.                 I = getPerIndices(cx,lf(2),locCol); x = x(:,I);
  176.             end
  177.             if ~isequal(locRow,'n') , 
  178.                 I = getPerIndices(rx,lf(1),locRow); x = x(I,:);
  179.             end
  180.         case 'per'            % Periodization.
  181.             if ~isequal(locCol,'n') ,   
  182.                 if rem(cx,2) , x(:,cx+1) = x(:,cx); cx = cx+1; end
  183.                 I = getPerIndices(cx,lf(2),locCol); x = x(:,I);
  184.             end
  185.             if ~isequal(locRow,'n') , 
  186.                 if rem(rx,2) , x(rx+1,:) = x(rx,:); rx = rx+1; end
  187.                 I = getPerIndices(rx,lf(1),locRow); x = x(I,:);
  188.             end
  189.         otherwise
  190.             errargt(mfilename,'Invalid Extension Mode!','msg');
  191.             error('*')
  192.       end
  193.  
  194.    case {'ar','addrow'}
  195.       if nargin<5, loc = 'b'; else , loc = testLoc(location(1)); end
  196.       location = [loc , 'n'];
  197.       x = wextend('2D',mode,x,lf,location);
  198.    case {'ac','addcol'}      
  199.       if nargin<5, loc = 'b'; else , loc = testLoc(location(1)); end
  200.       location = ['n' , loc];
  201.       x = wextend('2D',mode,x,lf,location);
  202. end
  203. %-------------------------------------------------------------------------%
  204. % Internal Function(s)
  205. %-------------------------------------------------------------------------%
  206. function location = testLoc(location)
  207. if ~ischar(location) , location = 'b'; return; end
  208. switch location
  209.    case {'n','l','u','b','r','d'}
  210.    otherwise , location = 'b';
  211. end
  212. %-------------------------------------------------------------------------%
  213. function [ext_L,ext_R] = getZPD_ext(nbr,nbc,location)
  214. switch location
  215.   case {'n'}     , ext_L = [];             ext_R = [];
  216.   case {'l','u'} , ext_L = zeros(nbr,nbc); ext_R = [];
  217.   case {'b'}     , ext_L = zeros(nbr,nbc); ext_R = zeros(nbr,nbc);
  218.   case {'r','d'} , ext_L = [];             ext_R = zeros(nbr,nbc);
  219. end
  220. %-------------------------------------------------------------------------%
  221. function [ext_L,ext_R] = getSP0_ext(type,x_L,x_R,lf,location)
  222. switch type(1)
  223.   case 'r' , ext_V = ones(1,lf);
  224.   case 'c' , ext_V = ones(lf,1);
  225. end
  226. switch location
  227.   case {'n'}     , ext_L = [];              ext_R = [];
  228.   case {'l','u'} , ext_L = kron(ext_V,x_L); ext_R = [];
  229.   case {'b'}     , ext_L = kron(ext_V,x_L); ext_R = kron(ext_V,x_R);
  230.   case {'r','d'} , ext_L = [];              ext_R = kron(ext_V,x_R);
  231. end
  232. %-------------------------------------------------------------------------%
  233. function [ext_L,ext_R] = getSP1_ext(type,x_L,x_R,d_L,d_R,lf,location)
  234. switch type(1)
  235.   case 'r' , ext_V0 = ones(1,lf); ext_VL = [lf:-1:1];  ext_VR = [1:lf];
  236.   case 'c' , ext_V0 = ones(lf,1); ext_VL = [lf:-1:1]'; ext_VR = [1:lf]';
  237. end
  238. switch location
  239.   case {'n'}
  240.     ext_L = [];
  241.     ext_R = [];
  242.   case {'l','u'}
  243.     ext_L = kron(ext_V0,x_L) + kron(ext_VL,d_L);
  244.     ext_R = [];
  245.   case {'b'}
  246.     ext_L = kron(ext_V0,x_L) + kron(ext_VL,d_L);
  247.     ext_R = kron(ext_V0,x_R) + kron(ext_VR,d_R);
  248.   case {'r','d'}
  249.     ext_L = [];
  250.     ext_R = kron(ext_V0,x_R) + kron(ext_VR,d_R);
  251. end
  252. %-------------------------------------------------------------------------%
  253. function I = getPerIndices(lx,lf,location)
  254. switch location
  255.     case {'n'}     , I = [1:lx];
  256.     case {'l','u'} , I = [lx-lf+1:lx , 1:lx];
  257.     case {'b'}     , I = [lx-lf+1:lx , 1:lx , 1:lf];
  258.     case {'r','d'} , I = [1:lx , 1:lf];
  259. end
  260. if lx<lf
  261.     I = mod(I,lx);
  262.     I(I==0) = lx;
  263. end
  264. %-------------------------------------------------------------------------%
  265. function I = getSymIndices(lx,lf,location)
  266. switch location
  267.     case {'n'}     , I = [1:lx];
  268.     case {'l','u'} , I = [lf:-1:1 , 1:lx];
  269.     case {'b'}     , I = [lf:-1:1 , 1:lx , lx:-1:lx-lf+1];
  270.     case {'r','d'} , I = [1:lx , lx:-1:lx-lf+1];
  271. end
  272. if lx<lf
  273.     K = (I<1);
  274.     I(K) = 1-I(K);
  275.     J = (I>lx);
  276.     while any(J)
  277.         I(J) = 2*lx+1-I(J);
  278.         K = (I<1);
  279.         I(K) = 1-I(K);
  280.         J = (I>lx);
  281.     end
  282. end
  283. %-------------------------------------------------------------------------%
  284. function Y = HP_SymExt(X,N,LOC) % Half-point Symmetrization.
  285. C = size(X,2);
  286. switch LOC
  287.     case {'n'}     , I = [1:C];
  288.     case {'l','u'} , I = [N:-1:1 , 1:C];
  289.     case {'b'}     , I = [N:-1:1 , 1:C , C:-1:C-N+1];
  290.     case {'r','d'} , I = [1:C , C:-1:C-N+1];
  291. end
  292. if C<N
  293.     K = (I<1);
  294.     I(K) = 1-I(K);
  295.     J = (I>C);
  296.     while any(J)
  297.         I(J) = 2*C+1-I(J);
  298.         K = (I<1);
  299.         I(K) = 1-I(K);
  300.         J = (I>C);
  301.     end
  302. end
  303. Y = X(:,I);
  304. %-------------------------------------------------------------------------%
  305. function Y = WP_SymExt(X,N,LOC) % Whole-point Symmetrization.
  306. [r,c] = size(X);
  307. while (N+1)>c
  308.     N = N-(c-1);
  309.     X = WP_SymExt(X,c-1,LOC);
  310.     [r,c] = size(X);
  311. end
  312. switch LOC
  313.   case {'l','u'} , I = [N+1:-1:2 , 1:c];
  314.   case {'b'}     , I = [N+1:-1:2 , 1:c , c-1:-1:c-N];
  315.   case {'r','d'} , I = [1:c , c-1:-1:c-N];
  316. end
  317. Y = X(:,I);
  318. %-------------------------------------------------------------------------%
  319. function Y = HP_AntiSymExt(X,N,LOC) % Half-point Anti-Symmetrization.
  320. [r,c] = size(X);
  321. while N>c
  322.     N = N-c;
  323.     X = HP_AntiSymExt(X,c,LOC);
  324.     [r,c] = size(X);
  325. end
  326. switch LOC
  327.     case {'l','u'}
  328.         Y = [fliplr(-X(:,1:N)), X];
  329.     case {'r','d'}
  330.         Y = [X , fliplr(-X(:,end-N+1:end))];
  331.     case 'b'
  332.         Y = [fliplr(-X(:,1:N)), X];
  333.         Y = [Y , fliplr(-X(:,end-N+1:end))];
  334. end
  335. %-------------------------------------------------------------------------%
  336. function Y = WP_AntiSymExt(X,N,LOC) % Whole-point Anti-Symmetrization.
  337. [r,c] = size(X);
  338. while (N+1)>c
  339.     N = N-(c-1);
  340.     X = WP_AntiSymExt(X,c-1,LOC);
  341.     [r,c] = size(X);
  342. end
  343. N = N+1;
  344. switch LOC
  345.     case {'l','u'}
  346.         Y = [fliplr(-X(:,2:N)+ 2*X(:,ones(1,N-1))) , X];
  347.     case {'r','d'}
  348.         Y = [X , fliplr(-X(:,end-N+1:end-1)+ 2*X(:,c*ones(1,N-1)))];
  349.     case 'b'
  350.         Y = [fliplr(-X(:,2:N)+ 2*X(:,ones(1,N-1))) , X];
  351.         Y = [Y , fliplr(-X(:,end-N+1:end-1)+ 2*X(:,c*ones(1,N-1)))];
  352. end
  353. %-------------------------------------------------------------------------%