scilab/modules/signal_processing/macros/unwrap.sci

   1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   2 // Copyright (C) - 2013 - Samuel GOUGEON
   3 //
   4 // Copyright (C) 2012 - 2016 - Scilab Enterprises
   5 //
   6 // This file is hereby licensed under the terms of the GNU GPL v2.0,
   7 // pursuant to article 5.3.4 of the CeCILL v.2.1.
   8 // This file was originally licensed under the terms of the CeCILL v2.1,
   9 // and continues to be available under such terms.
  10 // For more information, see the COPYING file which you should have received
  11 // along with this program.
  12
  13 function [retval, K] = unwrap(a, varargin)
  14     // [unwr, cuspPoints] = unwrap(a)        // row or column vector or matrix
  15     // [unwr, cuspPoints] = unwrap(a, jump)  // default jump = 2*%pi
  16     // jump=0 => = average local slope
  17     // [unwr, cuspPoints] = unwrap(a, jump, direc)
  18     //          direc="c": unwrap/unfold along columns
  19     //          direc="r": along rows
  20     //          else: along both directions (default)
  21     // [unwr, cuspPoints] = unwrap(a, "unfold")
  22
  23     // Initializations
  24     // --------------
  25     [lhs,rhs] = argn(0);
  26     transposed = %f
  27     unfold = %f     // unfold (1D) instead of unwrap (1D|2D)
  28     jump = 2*%pi
  29     direc  = ""
  30     retval = []
  31     K = []
  32
  33     // EXAMPLES
  34     // --------
  35     if rhs==0 then
  36         %_unwrap()      // display 1D examples
  37         halt(_("Press return to display 2D examples"))
  38         %_unwrap("2D")  // display 2D examples
  39         return
  40     end
  41
  42     // =================== CHECKING PARAMETERS ====================
  43     if type(a)~=1 then
  44         msg = _("%s: Wrong type for input argument #%d: Real expected.\n")
  45         error(msprintf(msg, "unwrap",1))
  46     end
  47     if rhs>1
  48         if typeof(varargin(1))=="string" then
  49             if varargin(1)=="unfold" then
  50                 unfold = %t
  51             else
  52                 msg = _("%s: Wrong value for input argument #%d: Must be in the set  {%s}.\n")
  53                 error(msprintf(msg, "unwrap",2,"""unfold"""))
  54             end
  55         else
  56             jump = varargin(1)
  57             if type(jump)~=1 then
  58                 msg = _("%s: Wrong type for input argument #%d: Real expected.\n")
  59                 error(msprintf(msg, "unwrap",2))
  60             else
  61                 jump = abs(jump(1))
  62             end
  63         end
  64         if ~unfold & rhs>=3 then
  65             direc = varargin(2)
  66             if typeof(direc)~="string" then
  67                 msg = _("%s: Wrong type for input argument #%d: String expected.\n")
  68                 error(msprintf(msg, "unwrap", 3))
  69             else
  70                 direc = direc(1)
  71                 if and(direc~=["r" "c"]) then
  72                     msg = _("%s: Wrong value for input argument #%d: Must be in the set  {%s}.\n")
  73                     error(msprintf(msg, "unwrap", 3, """r"",""c"""))
  74                 end
  75             end
  76         end
  77     end
  78
  79     if direc=="c" | size(a,2)==1
  80         a = a.'
  81         transposed = %t
  82     end
  83
  84     [m,n] = size(a)
  85     if (n < 4)
  86         if transposed then
  87             msg = _("%s: Wrong size for input argument #%d: Must have at least 4 rows.\n")
  88         else
  89             msg = _("%s: Wrong size for input argument #%d: Must have at least 4 columns.\n")
  90         end
  91         error(msprintf(msg, "unwrap",1))
  92     end
  93     if unfold & m>1 then
  94         msg = _("%s: Wrong size for input argument #%d: Vector expected.\n")
  95         error(msprintf(msg, "unwrap",1))
  96     end
  97
  98     // ============================ PROCESSING ==============================
  99
 100     // Columns #1 and #$ are duplicated:
 101     rae = [a(:,1)  a  a(:,$)];       // => n+2 columns
 102     // Derivative along rows (assuming equally spaced x values)
 103     d = rae(:,2:$) - rae(:, 1:$-1);     // n+1 columns
 104     clear rae
 105
 106     // ---------------------------  U N W R A P  ---------------------------
 107     if ~unfold then
 108         // jump (may be)  (n-1 columns)
 109         ju = d(:, 2:$-1)
 110         // average Local derivative (before/after jump) (n-1 columns) :
 111         avL = (d(:, 1:$-2)+d(:, 3:$))/2
 112         // where
 113         wh = abs(ju)>4*abs(avL) & d(:, 1:$-2).*d(:, 3:$)>0
 114         wh = abs(ju)>5*abs(avL)
 115         K = find(wh)
 116         d = d(:, 2:$-1)                  // n-1 columns
 117         if jump~=0 then
 118             d(K) = ju(K) - sign(ju(K)).*jump
 119             // Cleaning wrongly inserted jumps :
 120             k = find(abs(d(K)-mean(d))>5*stdev(d))
 121             if ~isempty(k)
 122                 d(K(k)) = d(K(k)) - sign(d(K(k))).*jump
 123             end
 124         else
 125             d(K) = avL(K)
 126         end
 127
 128         if and(size(a)~=1) & direc==""
 129             // 2D case: levels all rows according to the unwrapped 1st column
 130             p = unwrap(a(:,1),jump)
 131         else
 132             p  = a(:,1)
 133         end
 134         r = [p d]
 135
 136         // ---------------------------  U N F O L D  ---------------------------
 137     else
 138         // Local radius of curvature: calculation
 139         Rc = ((1+d(:,2:$).^2).^1.5) ./ abs(d(:,2:$)-d(:,1:$-1))
 140         // criterium to detect cusp points:
 141         K = find(Rc < min(mean(Rc)/30, min(Rc)/(5e3*%eps)))
 142         [nR,nC] = size(Rc)
 143
 144         // trimming edges (due to duplication of 1st & last columns):
 145         tmp = find(K==1 | K==nC)
 146         K(tmp) = []
 147
 148         // trimming duplicates (if any)
 149         tmp = find(K(2:$)-K(1:$-1)==1)
 150         K(tmp+1) = []
 151
 152         // fine targetting of cusp points:
 153         k = find(d(K).*d(K+1)>0)
 154         K(k) = K(k)+1
 155
 156         // We unfold:
 157         // The level and the slope at the left edge are taken as references.
 158         // The slope may be not of convenient sign. In this case, the
 159         //  unfolded profile will be upside down
 160         // One section over two must be kept as is, others must be flipped
 161
 162         nCusp = length(K)   // number of detected cusp points
 163         if nCusp>0 then
 164             flip = %t
 165             // Main loop over cusp points:
 166             for p = 1:nCusp
 167                 if flip then
 168                     iLeft = K(p)+1
 169                     if p<nCusp then
 170                         iRight = K(p+1)
 171                     else
 172                         iRight = nC
 173                     end
 174                     tmp = iLeft:iRight
 175                     d(tmp) = -d(tmp)        // the section is flipped
 176                     d(iLeft) = (d(iLeft-1)+d(iLeft+1))/2 // the reconnexion is smoothed
 177                 end
 178                 flip = ~flip    // one section / 2 is left as is
 179             end
 180         end
 181
 182         // trimming both edges =>  n-1 columns:
 183         d = d(:, 2:$-1)
 184
 185         // The first point is taken as anchor:
 186         r = [a(1)  d]
 187     end
 188
 189     // Building the new profile from its restored derivative
 190     retval = cumsum(r,"c")
 191
 192     // Post-processing
 193     if transposed then
 194         retval = retval.'
 195     end
 196 endfunction