move unwrap() => [signal_processing]
[scilab.git] / 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