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 %_unwrap(typexample)
  14     if argn(2)==0 then
  15         typexample = ""
  16     end
  17     if convstr(typexample)~="2d" then
  18         // 1-D EXAMPLES
  19         // -----------
  20         f = scf();
  21         f.figure_size = [800 1000];
  22         f.figure_position(2) = 0;
  23         f.figure_name = "unwrap() & ""unfold""" + _(": 1-D examples ");
  24         ax = gda();
  25         ax.y_label.font_size=2;
  26         drawlater()
  27
  28         // Original 1-D profile
  29         t = linspace(-4,4.2,800);
  30         alpha = t.^2 + t -1;
  31         subplot(5,2,1)
  32         titlepage("unwrap(): unwrap | unfold")
  33         subplot(5,2,2)
  34         plot(t,alpha)
  35         t2 = "$\text{Original profile: } \alpha=t^2+t-1$";
  36         ax = gca();
  37         ax.tight_limits = "on";
  38         yT = max(alpha) - strange(alpha)*0.1;
  39         xstring(0,yT,t2)
  40         e = gce();
  41         e.text_box_mode = "centered";
  42         e.font_size = 2;
  43
  44         // Loop over cases
  45         for i=1:4
  46             subplot(5,2,2*i+1)
  47             if i==1 then
  48                 // Wrapping by atan(tan())
  49                 ralpha = atan(tan(alpha));            // raw recovered alpha [pi]
  50                 ylabel("$\mbox{Y= }atan(tan(\alpha))$")
  51                 [u, K] = unwrap(ralpha, %pi);         // arctan
  52                 t2 = "$\text{unwrap(Y, \%pi)}$";
  53             elseif i==2
  54                 // Wrapping by modulo() + Y-shift
  55                 c = (rand(1,1)-0.5)*4;
  56                 ralpha = pmodulo(alpha, 5) + c;
  57                 ylabel("$\mbox{Y= }modulo(\alpha,\ 5)"+sprintf("%+5.2f",c)+"$")
  58                 [u, K] = unwrap(ralpha, 0);
  59                 t2 = "$\text{unwrap(Y, 0)}$";
  60             elseif i==3
  61                 // Folding by asin(sin()) + Y-shift
  62                 ralpha = 1+asin(sin(alpha));          // raw recovered alpha [2.pi]
  63                 ylabel("$\mbox{Y= }1+asin(sin(\alpha))$")
  64                 [u, K] = unwrap(ralpha, "unfold");
  65                 t2 = "$\text{unwrap(Y,""unfold"")}$";
  66             else
  67                 // Folding by acos(cos()) + trend
  68                 ralpha = 1+alpha/10+acos(cos(alpha)); // raw recovered alpha [2.pi]
  69                 ylabel("$\mbox{Y= }1+\frac{\alpha}{10}+acos(cos(\alpha))$")
  70                 [u, K] = unwrap(ralpha, "unfold");
  71                 t2 = "$\text{unwrap(Y,""unfold"")}$";
  72             end
  73             // Plotting the profile to be processed
  74             plot(t, ralpha)
  75             // Staring the breakpoints or the cusp points on the curve:
  76             if K~=[]
  77                 plot(t(K), ralpha(K),"*")
  78             end
  79             // Plotting the processed (unwrapped/unfolded) profile:
  80             subplot(5,2,2*i+2)
  81             plot(t,u)
  82             ax = gca();
  83             ax.tight_limits = "on";
  84             // Adding a legend:
  85             yT = max(u) - strange(u)*0.2;
  86             xstring(0,yT,t2)
  87             e = gce();
  88             e.text_box_mode = "centered";
  89             e.font_size = 2;
  90         end
  91         sda();
  92         drawnow()
  93
  94     else
  95         // 2D EXAMPLES
  96         // -----------
  97         ax = gda();
  98         ax.title.font_size = 2;
  99         f = scf();
 100         f.color_map = hotcolormap(100);
 101         f.figure_size = [475 1050];
 102         f.figure_position(2) = 0;
 103         f.figure_name = "unwrap()" + _(": 2-D examples");
 104         drawlater()
 105
 106         nx = 300;
 107         ny = 400;
 108         rmax = 8.8;
 109         x = linspace(-rmax/2, rmax/2, nx)-1;
 110         y = linspace(-rmax/2, rmax/2, ny)+1;
 111         [X, Y] = meshgrid(x,y);
 112         for ex=0:1   // examples
 113             // Original surface
 114             // Generating the surface
 115             if ex==0 then
 116                 z = X.^2 + Y.^2;
 117             else
 118                 z = sqrt(0.3+sinc(sqrt(z)*3))*17-7;
 119             end
 120             // PLotting it in 3D
 121             subplot(4,2,1+ex)
 122             surf(x, y, z)
 123             title("Original profile Z")
 124             e = gce();
 125             e.thickness = 0; // removes the mesh
 126             e.parent.tight_limits = "on";
 127
 128             // WRAPPED surface (flat)
 129             m = 2.8;
 130             zw = pmodulo(z, m); // wraps it
 131             subplot(4,2,3+ex)
 132             grayplot(x, y, zw')
 133             title(msprintf("Zw = pmodulo(Z, %g)  (flat)",m))
 134             ax0 = gca();
 135             ax0.tight_limits = "on";
 136
 137             // UNWRAPPED surfaces (flat):
 138             // in both directions:
 139             u = unwrap(zw, 0);
 140             subplot(4,2,5+ex)
 141             grayplot(x, y, u')
 142             title(msprintf("unwrap(Zw, %g)  (flat)", 0))
 143             ax3 = gca();
 144             ax3.tight_limits = "on";
 145
 146             if ex==0, direc="r", else direc="c", end
 147             // along a single direction:
 148             u = unwrap(zw, m, direc);
 149             subplot(4,2,7+ex)
 150             grayplot(x, y, u')
 151             title(msprintf("unwrap(Zw, %g, ""%s"")  (flat)",m,direc))
 152             ax1 = gca();
 153             ax1.tight_limits = "on";
 154         end
 155         sda();
 156         drawnow()
 157     end
 158 endfunction