summaryrefslogtreecommitdiffstats
path: root/scilab/modules/signal_processing
diff options
context:
space:
mode:
authorClément David <clement.david@esi-group.com>2021-07-15 15:50:48 +0200
committerClément David <clement.david@esi-group.com>2021-07-15 15:50:48 +0200
commit81a9cc049332de0c712cf56da585fcd25c8e59e3 (patch)
treef365c25a05c1782a16b2985bf7481116c106952d /scilab/modules/signal_processing
parentb0937f19e4b8ddf416ca9a9a433bcbbd3f4ef2c0 (diff)
parent90b4ecb9b57ac6ad4e00d4a99b0d8bd1a6424403 (diff)
downloadscilab-master.zip
scilab-master.tar.gz
Merge remote-tracking branch 'origin/6.1'HEADmaster
Change-Id: I5d26fc380a28efe0bb6d0096fe9684b274b55bfe
Diffstat (limited to 'scilab/modules/signal_processing')
-rw-r--r--scilab/modules/signal_processing/Makefile.in2
-rw-r--r--scilab/modules/signal_processing/help/en_US/filters/filter.xml3
-rwxr-xr-xscilab/modules/signal_processing/help/en_US/filters/sgolay.xml259
-rwxr-xr-xscilab/modules/signal_processing/help/en_US/filters/sgolaydiff.xml188
-rwxr-xr-xscilab/modules/signal_processing/help/en_US/filters/sgolayfilt.xml180
-rw-r--r--scilab/modules/signal_processing/help/en_US/unwrap.xml306
-rw-r--r--scilab/modules/signal_processing/help/images/unwrap_2D.pngbin0 -> 108181 bytes
-rw-r--r--scilab/modules/signal_processing/help/ja_JP/unwrap.xml389
-rw-r--r--scilab/modules/signal_processing/macros/%_unwrap.sci158
-rw-r--r--scilab/modules/signal_processing/macros/sgolay.sci86
-rw-r--r--scilab/modules/signal_processing/macros/sgolaydiff.sci54
-rw-r--r--scilab/modules/signal_processing/macros/sgolayfilt.sci66
-rw-r--r--scilab/modules/signal_processing/macros/unwrap.sci196
-rw-r--r--scilab/modules/signal_processing/tests/nonreg_tests/bug_16305.tst23
-rw-r--r--scilab/modules/signal_processing/tests/unit_tests/unwrap.dia.ref77
-rw-r--r--scilab/modules/signal_processing/tests/unit_tests/unwrap.tst84
16 files changed, 2071 insertions, 0 deletions
diff --git a/scilab/modules/signal_processing/Makefile.in b/scilab/modules/signal_processing/Makefile.in
index d38d068..ed59e26 100644
--- a/scilab/modules/signal_processing/Makefile.in
+++ b/scilab/modules/signal_processing/Makefile.in
@@ -493,9 +493,11 @@ NMEDIT = @NMEDIT@
493OBJDUMP = @OBJDUMP@ 493OBJDUMP = @OBJDUMP@
494OBJEXT = @OBJEXT@ 494OBJEXT = @OBJEXT@
495OCAMLC = @OCAMLC@ 495OCAMLC = @OCAMLC@
496OCAMLCFLAGS = @OCAMLCFLAGS@
496OCAMLDEP = @OCAMLDEP@ 497OCAMLDEP = @OCAMLDEP@
497OCAMLLEX = @OCAMLLEX@ 498OCAMLLEX = @OCAMLLEX@
498OCAMLOPT = @OCAMLOPT@ 499OCAMLOPT = @OCAMLOPT@
500OCAMLOPTFLAGS = @OCAMLOPTFLAGS@
499OCAMLYACC = @OCAMLYACC@ 501OCAMLYACC = @OCAMLYACC@
500OPENMPI_CC = @OPENMPI_CC@ 502OPENMPI_CC = @OPENMPI_CC@
501OPENMPI_CFLAGS = @OPENMPI_CFLAGS@ 503OPENMPI_CFLAGS = @OPENMPI_CFLAGS@
diff --git a/scilab/modules/signal_processing/help/en_US/filters/filter.xml b/scilab/modules/signal_processing/help/en_US/filters/filter.xml
index a35eebe..6c7d479 100644
--- a/scilab/modules/signal_processing/help/en_US/filters/filter.xml
+++ b/scilab/modules/signal_processing/help/en_US/filters/filter.xml
@@ -133,6 +133,9 @@
133 <member> 133 <member>
134 <link linkend="ltitr">ltitr</link> 134 <link linkend="ltitr">ltitr</link>
135 </member> 135 </member>
136 <member>
137 <link linkend="sgolayfilt">sgolayfilt</link>
138 </member>
136 </simplelist> 139 </simplelist>
137 </refsection> 140 </refsection>
138</refentry> 141</refentry>
diff --git a/scilab/modules/signal_processing/help/en_US/filters/sgolay.xml b/scilab/modules/signal_processing/help/en_US/filters/sgolay.xml
new file mode 100755
index 0000000..e6d742a
--- /dev/null
+++ b/scilab/modules/signal_processing/help/en_US/filters/sgolay.xml
@@ -0,0 +1,259 @@
1<?xml version="1.0" encoding="UTF-8"?>
2<!--
3 This file is part of the Cardiovascular Wawes Analysis toolbox
4 Copyright (C) 2014 - INRIA - Serge Steer
5 This file must be used under the terms of the CeCILL.
6 This source file is licensed as described in the file COPYING, which
7 you should have received as part of this distribution. The terms
8 are also available at
9 http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10-->
11<refentry xml:id="sgolay" xml:lang="en" xmlns="http://docbook.org/ns/docbook"
12 xmlns:xlink="http://www.w3.org/1999/xlink"
13 xmlns:svg="http://www.w3.org/2000/svg"
14 xmlns:scilab="http://www.scilab.org"
15 xmlns:ns4="http://www.w3.org/1999/xhtml"
16 xmlns:mml="http://www.w3.org/1998/Math/MathML"
17 xmlns:db="http://docbook.org/ns/docbook">
18 <refnamediv>
19 <refname>sgolay</refname>
20
21 <refpurpose>Savitzky-Golay Filter Design.</refpurpose>
22 </refnamediv>
23
24 <refsynopsisdiv>
25 <title>Calling Sequence</title>
26
27 <synopsis>[B [,C]] = sgolay(k,nf [,w])</synopsis>
28 </refsynopsisdiv>
29
30 <refsection>
31 <title>Arguments</title>
32
33 <variablelist>
34 <varlistentry>
35 <term>k</term>
36
37 <listitem>
38 <para>a positive scalar with integer value: the fitting polynomial
39 degree.</para>
40 </listitem>
41 </varlistentry>
42
43 <varlistentry>
44 <term>nf</term>
45
46 <listitem>
47 <para>a positive scalar with integer value: the filter length, must
48 be odd and greater the k+1.</para>
49 </listitem>
50 </varlistentry>
51
52 <varlistentry>
53 <term>w</term>
54
55 <listitem>
56 <para>a real vector of length nf with positive entries: the weights.
57 If omitted no weights are applied.</para>
58 </listitem>
59 </varlistentry>
60
61 <varlistentry>
62 <term>B</term>
63
64 <listitem>
65 <para>a real n by n array: the set of filter coefficients: the early
66 rows of B smooth based on future values and later rows smooth based
67 on past values, with the middle row using half future and half past.
68 In particular, you can use row i to estimate x(k) based on the i-1
69 preceding values and the n-i following values of x values as y(k) =
70 B(i,:) * x(k-i+1:k+n-i).</para>
71 </listitem>
72 </varlistentry>
73
74 <varlistentry>
75 <term>C</term>
76
77 <listitem>
78 <para>a real n by k+1 array: the matrix of differentiation filters.
79 Each column of C is a differentiation filter for derivatives of
80 order P-1 where P is the column index. Given a signal X with length
81 nf, an estimate of the P-th order derivative of its middle value can
82 be found from: <para/> (P) X ((nf+1)/2) = P!*C(:,P+1)'*X</para>
83 </listitem>
84 </varlistentry>
85 </variablelist>
86 </refsection>
87
88 <refsection>
89 <title>Description</title>
90
91 <para>This function computes the smoothing FIR Savitzky-Golay
92 Filter. This is achieved by fitting successive sub-sets of
93 adjacent data points with a low-degree polynomial by the method of
94 linear least squares.</para>
95
96 <para>This filter can also be used to approximate numerical
97 derivatives.</para>
98 </refsection>
99
100 <refsection>
101 <title>Examples</title>
102
103 <variablelist>
104 <varlistentry>
105 <listitem>
106 <para>The sgolay function can be used to construct FIR filter with
107 no group delay. The following instructions explains what the <link
108 linkend="sgolayfilt">sgolayfilt</link> function does:</para>
109
110 <programlisting role="example"><![CDATA[
111dt=0.01;
112t = (0:0.01:4*%pi)';
113x = sin(t)+0.05*rand(t,'normal');
114
115nf = 61;
116[B,C] = sgolay(3,nf);
117nfs2 = floor(nf/2);
118
119//compute the middle part of the filtered signal
120xfm = filter(B(nfs2+1,:), 1, x);
121xfm = xfm(nf:$);
122//compute the left part of the filtered signal
123xfl = B(1:nfs2,:)*x(1:nf);
124//compute the right part of the filtered signal
125xfr = B(nfs2+2:nf,:)*x($-nf+1:$);
126clf;ax=gca();
127plot(t,x,'b',t,[xfl;xfm;xfr],'r');
128e=gce();set(e.children(1),"thickness",2);
129legend(["Raw","Filtered"]);
130 ]]></programlisting>
131
132 <para>The above script generate this graph:</para>
133
134 <scilab:image><![CDATA[
135dt=0.01;
136t = (0:0.01:4*%pi)';
137x = sin(t)+0.05*rand(t,'normal');
138
139nf = 61;
140[B,C] = sgolay(3,nf);
141nfs2 = floor(nf/2);
142
143//compute the middle part of the filtered signal
144xfm = filter(B(nfs2+1,:), 1, x);
145xfm = xfm(nf:$);
146//compute the left part of the filtered signal
147xfl = B(1:nfs2,:)*x(1:nf);
148//compute the right part of the filtered signal
149xfr = B(nfs2+2:nf,:)*x($-nf+1:$);
150clf;ax=gca();
151plot(t,x,'b',t,[xfl;xfm;xfr],'r');
152e=gce();set(e.children(1),"thickness",2);
153legend(["Raw","Filtered"]);
154 ]]></scilab:image>
155 </listitem>
156 </varlistentry>
157
158 <varlistentry>
159 <listitem>
160 <para>It can also be used to obtain approximate derivatives:</para>
161
162 <programlisting role="example"><![CDATA[
163dt=0.01;
164t = (0:0.01:4*%pi)';
165x = sin(t)+0.03*rand(t,'normal');
166
167nf = 61;
168[B,C] = sgolay(3,nf);
169nfs2 = floor(nf/2);
170
171dx=-conv(x,C(:,2)',"valid");
172d2x=2*conv(x,C(:,3)',"valid");
173
174clf();
175subplot(211)
176plot(t(2:$),diff(x)/dt,'g',t,cos(t),'b',t(nfs2+1:$-nfs2),dx/dt,'r');
177legend(["diff","theoretical","sgolay"]);
178title("First order differentiation")
179subplot(212)
180plot(t,-sin(t),'b',t(nfs2+1:$-nfs2),d2x/(dt^2),'r');
181legend(["theoretical","sgolay"]);
182title("Second order differentiation")
183 ]]></programlisting>
184
185 <scilab:image><![CDATA[
186dt=0.01;
187t = (0:0.01:4*%pi)';
188x = sin(t)+0.03*rand(t,'normal');
189
190nf = 61;
191[B,C] = sgolay(3,nf);
192nfs2 = floor(nf/2);
193
194dx=-conv(x,C(:,2)',"valid");
195d2x=2*conv(x,C(:,3)',"valid");
196
197clf();
198subplot(211)
199plot(t(2:$),diff(x)/dt,'g',t,cos(t),'b',t(nfs2+1:$-nfs2),dx/dt,'r');
200legend(["diff","theoretical","sgolay"]);
201title("First order differentiation")
202subplot(212)
203plot(t,-sin(t),'b',t(nfs2+1:$-nfs2),d2x/(dt^2),'r');
204legend(["theoretical","sgolay"]);
205title("Second order differentiation")
206 ]]></scilab:image>
207 </listitem>
208 </varlistentry>
209 </variablelist>
210 </refsection>
211
212 <refsection>
213 <title>See Also</title>
214
215 <simplelist type="inline">
216 <member><link linkend="sgolayfilt">sgolayfilt</link></member>
217
218 <member><link linkend="conv">conv</link></member>
219 </simplelist>
220 </refsection>
221
222 <refsection>
223 <title>Authors</title>
224
225 <simplelist type="vert">
226 <member>Serge Steer, INRIA</member>
227 </simplelist>
228 </refsection>
229
230 <refsection>
231 <title>Bibliography</title>
232
233 <para>Abraham Savitzky et Marcel J. E. Golay, « Smoothing and
234 Differentiation of Data by Simplified Least Squares Procedures »,
235 Analytical Chemistry, vol. 8, no 36,‎ 1964, p. 1627–1639 (DOI
236 10.1021/ac60214a047)</para>
237
238 <para><ulink
239 url="http://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter">http://en.wikipedia.org/wiki/Savitzky-Golay_filter</ulink>.</para>
240 </refsection>
241
242 <refsection>
243 <title>History</title>
244
245 <revhistory>
246 <revision>
247 <revnumber>0.0</revnumber>
248
249 <revdescription>Function added</revdescription>
250 </revision>
251 </revhistory>
252 </refsection>
253
254 <refsection>
255 <title>Used Functions</title>
256
257 <para><link linkend="pinv">pinv</link></para>
258 </refsection>
259</refentry>
diff --git a/scilab/modules/signal_processing/help/en_US/filters/sgolaydiff.xml b/scilab/modules/signal_processing/help/en_US/filters/sgolaydiff.xml
new file mode 100755
index 0000000..c50bcb0
--- /dev/null
+++ b/scilab/modules/signal_processing/help/en_US/filters/sgolaydiff.xml
@@ -0,0 +1,188 @@
1<?xml version="1.0" encoding="UTF-8"?>
2<!--
3 Copyright (C) 2014 - INRIA - Serge Steer
4 This file must be used under the terms of the CeCILL.
5 This source file is licensed as described in the file COPYING, which
6 you should have received as part of this distribution. The terms
7 are also available at
8 http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
9
10-->
11<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:lang="en" xml:id="sgolaydiff">
12 <refnamediv>
13 <refname>sgolaydiff</refname>
14 <refpurpose>Numerical derivatives computation using Savitzky-Golay filter.</refpurpose>
15 </refnamediv>
16 <refsynopsisdiv>
17 <title>Calling Sequence</title>
18 <synopsis>D = sgolaydiff(X,order,k,nf [,w])</synopsis>
19 <synopsis>D = sgolaydiff(X,order,C)</synopsis>
20 </refsynopsisdiv>
21 <refsection>
22 <title>Arguments</title>
23 <variablelist>
24 <varlistentry>
25 <term>X</term>
26 <listitem>
27 <para>
28 It can be either a real vector or a general real array. If
29 it is an array the derivation is applied along the first
30 non singleton dimension.
31 </para>
32 </listitem>
33 </varlistentry>
34 <varlistentry>
35 <term>order</term>
36 <listitem>
37 <para>
38 a positive scalar with integer value, the order of derivation.
39 </para>
40 </listitem>
41 </varlistentry>
42 <varlistentry>
43 <term>C</term>
44 <listitem>
45 <para>a real nf by k+1 array: the matrix of differentiation
46 filters computed by the <link linkend="sgolay">sgolay</link> function.</para>
47 </listitem>
48 </varlistentry>
49 <varlistentry>
50 <term>k</term>
51 <listitem>
52 <para>a positive scalar with integer value: the fitting polynomial
53 degree.</para>
54 </listitem>
55 </varlistentry>
56 <varlistentry>
57 <term>nf</term>
58 <listitem>
59 <para>a positive scalar with integer value: the filter length, must
60 be odd and greater the k+1.</para>
61 </listitem>
62 </varlistentry>
63 <varlistentry>
64 <term>w</term>
65
66 <listitem>
67 <para>a real vector of length nf with positive entries: the weights.
68 If omitted no weights are applied.</para>
69 </listitem>
70 </varlistentry>
71 <varlistentry>
72 <term>D</term>
73 <listitem>
74 <para>A vector with identical shape as X, the estimated derivative.</para>
75 </listitem>
76 </varlistentry>
77 </variablelist>
78 </refsection>
79 <refsection>
80 <title>Description</title>
81 <para>
82 <para>This function computes numerical derivatives using the FIR
83 Savitzky-Golay smoothing Filter. This is achieved by fitting
84 successive sub-sets of adjacent data points with a low-degree
85 polynomial by the method of linear least squares.</para>
86 </para>
87 </refsection>
88 <refsection>
89 <title>More information</title>
90 <caution><para>order must be less than k</para></caution>
91 <note><para>The (nf−1)/2 first and last derivative values are
92 estimated by adding to the input signal, in reverse order and vertically symmetrized, copies of the first
93 (nf−1)/2 points at the beginning and copies of the last (nf−1)/2
94 points at the end .</para></note>
95 </refsection>
96 <refsection>
97 <title>Examples</title>
98 <programlisting role="example"><![CDATA[
99//generate a noisy signal
100dt=0.01;
101t = (0:0.01:4*%pi)';
102x = sin(t)+0.03*rand(t,"normal");
103
104clf();
105//first order derivative
106subplot(211)
107dx=sgolaydiff(x,1,3,61);
108plot(t(2:$),diff(x)/dt,'g',t,cos(t),'b',t,dx/dt,'r');
109legend(["diff","theoretical","sgolaydiff"]);
110//second order derivative
111subplot(212)
112d2x=sgolaydiff(x,2,3,101);
113plot(t,-sin(t),'b',t,d2x/(dt^2),'r');
114legend(["theoretical","sgolaydiff"]);
115title("Second order differentiation")
116 ]]></programlisting>
117 <scilab:image><![CDATA[
118//generate a noisy signal
119dt=0.01;
120t = (0:0.01:4*%pi)';
121x = sin(t)+0.03*rand(t,"normal");
122
123clf();
124//first order derivative
125subplot(211)
126dx=sgolaydiff(x,1,3,61);
127plot(t(2:$),diff(x)/dt,'g',t,cos(t),'b',t,dx/dt,'r');
128legend(["diff","theoretical","sgolaydiff"]);
129//second order derivative
130subplot(212)
131d2x=sgolaydiff(x,2,3,101);
132plot(t,-sin(t),'b',t,d2x/(dt^2),'r');
133legend(["theoretical","sgolaydiff"]);
134title("Second order differentiation")
135 ]]></scilab:image>
136 </refsection>
137 <refsection>
138 <title>See Also</title>
139 <simplelist type="inline">
140 <member>
141 <link linkend="sgolay">sgolay</link>
142 </member>
143 <member>
144 <link linkend="sgolayfilt">sgolayfilt</link>
145 </member>
146 <member>
147 <link linkend="diff">diff</link>
148 </member>
149 <member>
150 <link linkend="numdiff">numdiff</link>
151 </member>
152 </simplelist>
153 </refsection>
154 <refsection>
155 <title>Authors</title>
156 <simplelist type="vert">
157 <member>Serge Steer, INRIA</member>
158 </simplelist>
159 </refsection>
160 <refsection>
161 <title>Bibliography</title>
162 <para>Abraham Savitzky et Marcel J. E. Golay, « Smoothing and
163 Differentiation of Data by Simplified Least Squares Procedures »,
164 Analytical Chemistry, vol. 8, no 36,‎ 1964, p. 1627–1639 (DOI
165 10.1021/ac60214a047)</para>
166
167 <para><ulink
168 url="http://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter">http://en.wikipedia.org/wiki/Savitzky-Golay_filter</ulink>.</para>
169 </refsection>
170
171 <refsection>
172 <title>History</title>
173 <revhistory>
174 <revision>
175 <revnumber>0.0</revnumber>
176 <revdescription>Function added</revdescription>
177 </revision>
178 </revhistory>
179 </refsection>
180 <refsection>
181 <title>Used Functions</title>
182 <para> This function is based on the
183 <link linkend="sgolay">sgolay</link> function which computes the matrix of differentiation
184 filters and on the
185 <link linkend="conv">conv</link> function to compute the smoothed derivative/
186 </para>
187 </refsection>
188</refentry>
diff --git a/scilab/modules/signal_processing/help/en_US/filters/sgolayfilt.xml b/scilab/modules/signal_processing/help/en_US/filters/sgolayfilt.xml
new file mode 100755
index 0000000..15ea486
--- /dev/null
+++ b/scilab/modules/signal_processing/help/en_US/filters/sgolayfilt.xml
@@ -0,0 +1,180 @@
1<?xml version="1.0" encoding="UTF-8"?>
2<!--
3 This file is part of the Cardiovascular Wawes Analysis toolbox
4 Copyright (C) 2014 - INRIA - Serge Steer
5 This file must be used under the terms of the CeCILL.
6 This source file is licensed as described in the file COPYING, which
7 you should have received as part of this distribution. The terms
8 are also available at
9 http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10-->
11<refentry xml:id="sgolayfilt" xml:lang="en"
12 xmlns="http://docbook.org/ns/docbook"
13 xmlns:xlink="http://www.w3.org/1999/xlink"
14 xmlns:svg="http://www.w3.org/2000/svg"
15 xmlns:scilab="http://www.scilab.org"
16 xmlns:ns4="http://www.w3.org/1999/xhtml"
17 xmlns:mml="http://www.w3.org/1998/Math/MathML"
18 xmlns:db="http://docbook.org/ns/docbook">
19 <refnamediv>
20 <refname>sgolayfilt</refname>
21
22 <refpurpose>Filter signal using Savitzky-Golay Filter.</refpurpose>
23 </refnamediv>
24
25 <refsynopsisdiv>
26 <title>Calling Sequence</title>
27
28 <synopsis>xf=sgolayfilt(X,B)</synopsis>
29
30 <synopsis>xf=sgolayfilt(X,k,nf [,w]) </synopsis>
31 </refsynopsisdiv>
32
33 <refsection>
34 <title>Arguments</title>
35
36 <variablelist>
37 <varlistentry>
38 <term>X</term>
39
40 <listitem>
41 <para>a real array. If it is a vector it contains the signal to be
42 filtered else each column of X is filtered.</para>
43 </listitem>
44 </varlistentry>
45
46 <varlistentry>
47 <term>B</term>
48
49 <listitem>
50 <para>a real n by n array: the set of filter coefficients produced
51 by the <link linkend="sgolay">sgolay</link> function.</para>
52 </listitem>
53 </varlistentry>
54
55 <varlistentry>
56 <term>k</term>
57
58 <listitem>
59 <para>a positive scalar with integer value: the fitting polynomial
60 degree.</para>
61 </listitem>
62 </varlistentry>
63
64 <varlistentry>
65 <term>k</term>
66
67 <listitem>
68 <para>a positive scalar with integer value: the filter length, must
69 be odd and greater the k+1.</para>
70 </listitem>
71 </varlistentry>
72
73 <varlistentry>
74 <term>xf</term>
75
76 <listitem>
77 <para>a real array: the filtered signal(s).</para>
78 </listitem>
79 </varlistentry>
80 </variablelist>
81 </refsection>
82
83 <refsection>
84 <title>Description</title>
85
86 <para>This function applies a sgolay filter which can be specified either
87 by its set of filter coefficients as produced by <link
88 linkend="sgolay">sgolay</link> either by the degree of the polynomial and
89 the filter length.</para>
90 <para>According to the <ulink
91 url="http://www.hpl.hp.com/techreports/2010/HPL-2010-109.pdf">article</ulink> the -3dB cutoff frequency can be approximated by
92 (k+1)/(3.2*nf-4.6).</para>
93 </refsection>
94
95 <refsection>
96 <title>More information</title>
97
98 <warning>
99 <para>Signal size should be greater than the filter length</para>
100 </warning>
101 </refsection>
102
103 <refsection>
104 <title>Examples</title>
105
106 <programlisting role="example"><![CDATA[
107dt=0.01;
108t = (0:0.01:4*%pi)';
109x = sin(t)+0.05*rand(t,'normal');
110clf;plot(t,x,'b',t,sgolayfilt(x,3,51),'r');
111e=gce();set(e.children(1),"thickness",2);
112legend(["Raw","Filtered"]);
113 ]]></programlisting>
114
115 <scilab:image><![CDATA[
116dt=0.01;
117t = (0:0.01:4*%pi)';
118x = sin(t)+0.05*rand(t,'normal');
119clf;plot(t,x,'b',t,sgolayfilt(x,3,51),'r');
120e=gce();set(e.children(1),"thickness",2);
121legend(["Raw","Filtered"]);
122 ]]></scilab:image>
123 </refsection>
124
125 <refsection>
126 <title>See Also</title>
127
128 <simplelist type="inline">
129 <member><link linkend="sgolay">sgolay</link></member>
130
131 <member><link linkend="filter">filter</link></member>
132
133 <member><link linkend="wfir">wfir</link></member>
134 </simplelist>
135 </refsection>
136
137 <refsection>
138 <title>Authors</title>
139
140 <simplelist type="vert">
141 <member>Serge Steer, INRIA</member>
142 </simplelist>
143 </refsection>
144
145 <refsection>
146 <title>Bibliography</title>
147
148 <para>Abraham Savitzky et Marcel J. E. Golay, « Smoothing and
149 Differentiation of Data by Simplified Least Squares Procedures »,
150 Analytical Chemistry, vol. 8, no 36,‎ 1964, p. 1627–1639 (DOI
151 10.1021/ac60214a047)</para>
152
153 <para><ulink
154 url="http://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter">http://en.wikipedia.org/wiki/Savitzky-Golay_filter</ulink>.</para>
155
156 <para><ulink
157 url="http://www.hpl.hp.com/techreports/2010/HPL-2010-109.pdf">On
158 the Frequency-Domain Properties of Savitzky-Golay
159 Filters</ulink>.</para>
160 </refsection>
161
162 <refsection>
163 <title>History</title>
164
165 <revhistory>
166 <revision>
167 <revnumber>0.0</revnumber>
168
169 <revdescription>Function sgolayfilt added.</revdescription>
170 </revision>
171 </revhistory>
172 </refsection>
173
174 <refsection>
175 <title>Used Functions</title>
176
177 <para><link linkend="filter">filter</link>,<link
178 linkend="sgolay">sgolay</link></para>
179 </refsection>
180</refentry>
diff --git a/scilab/modules/signal_processing/help/en_US/unwrap.xml b/scilab/modules/signal_processing/help/en_US/unwrap.xml
new file mode 100644
index 0000000..58ef427
--- /dev/null
+++ b/scilab/modules/signal_processing/help/en_US/unwrap.xml
@@ -0,0 +1,306 @@
1<?xml version="1.0" encoding="UTF-8"?>
2<!--
3 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
4 * Copyright (C) 2013 - Samuel GOUGEON
5 *
6 * Copyright (C) 2012 - 2016 - Scilab Enterprises
7 *
8 * This file is hereby licensed under the terms of the GNU GPL v2.0,
9 * pursuant to article 5.3.4 of the CeCILL v.2.1.
10 * This file was originally licensed under the terms of the CeCILL v2.1,
11 * and continues to be available under such terms.
12 * For more information, see the COPYING file which you should have received
13 * along with this program.
14 *
15 -->
16<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns5="http://www.w3.org/1999/xhtml" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:id="unwrap" xml:lang="en">
17 <refnamediv>
18 <refname>unwrap</refname>
19 <refpurpose>unwrap a Y(x) profile or a Z(x,y) surface. Unfold a Y(x) profile</refpurpose>
20 </refnamediv>
21 <refsynopsisdiv>
22 <title>Syntaxes</title>
23 <synopsis>unwrap() // runs some examples
24 [U, breakPoints] = unwrap(Y)
25 [U, breakPoints] = unwrap(Y, z_jump)
26 [U, cuspPoints] = unwrap(Y, "unfold")
27 U = unwrap(Z)
28 U = unwrap(Z, z_jump)
29 U = unwrap(Z, z_jump, dir)
30 </synopsis>
31 </refsynopsisdiv>
32 <refsection>
33 <title>Arguments</title>
34 <variablelist>
35 <varlistentry>
36 <term>Y</term>
37 <listitem>
38 <para>Vector of real numbers: the profile to unwrap or unfold. Implicit abscissae X are assumed to be equispaced.</para>
39 </listitem>
40 </varlistentry>
41 <varlistentry>
42 <term>Z</term>
43 <listitem>
44 <para>Matrix of real numbers: the surface to unwrap. Implicit abscissae (X,Y) are assumed to be cartesian and equispaced (constant steps may be different along X versus along Y).</para>
45 </listitem>
46 </varlistentry>
47 <varlistentry>
48 <term>z_jump</term>
49 <listitem>
50 <para>
51 Scalar real positive number used in unwrapping mode: the jump's height applied at breakpoints, performing the unwrapping. Only its absolute value is considered. The jump actually applied has the sign of the slope on both sides of each breakpoint. The default value is <literal>z_jump = 2*%pi</literal>. The special value <literal>z_jump = 0</literal> applies jumps equal to the average slope around each breakpoint, restoring a continuous slope over the whole profile or surface.
52 </para>
53 </listitem>
54 </varlistentry>
55 <varlistentry>
56 <term>dir</term>
57 <listitem>
58 <para>"c" | "r" | "" (default): direction along which unwrapping is performed. "c" unwraps along columns, "r" unwraps along rows, "" unwraps in both directions.</para>
59 </listitem>
60 </varlistentry>
61 <varlistentry>
62 <term>"unfold"</term>
63 <listitem>
64 <para>Provide this switch to unfold the given curve if it is folded, instead of unwrapping it.</para>
65 </listitem>
66 </varlistentry>
67 <varlistentry>
68 <term>U</term>
69 <listitem>
70 <para>
71 Unwrapped profile or surface, or unfolded profile. <varname>U</varname> has the same sizes as <varname>Y</varname> or <varname>Z</varname>.
72 </para>
73 </listitem>
74 </varlistentry>
75 <varlistentry>
76 <term>breakPoints, cuspPoints</term>
77 <listitem>
78 <para>
79 Vector of indices of points in <varname>Y</varname> where wrapping or folding has been detected and processed.
80 </para>
81 </listitem>
82 </varlistentry>
83 </variablelist>
84 </refsection>
85 <refsection>
86 <title>Description</title>
87 <para>UNWRAPPING</para>
88 <para>
89 <function>unwrap()</function> will be useful to process profiles or even surfaces wrapped for instance by a periodic and monotonic function such as
90 <literal>Y = modulo(X,w)</literal> or <literal>Y = atan(tan(X))</literal>. It aims to somewhat invert these functions, recovering the input <literal>X</literal> over it full range instead of the limited <literal>w</literal> or <literal>[-%pi/2, %pi/2]</literal> one.
91 </para>
92 <para>A breakpoint of a wrapped profile is detected as a point where slopes on both neighbouring sides of the point are almost equal but much smaller (in absolute value) from and opposite to the slope at the considered point: at the point, there is a jump breaking and opposed to the neighbouring slope.
93 </para>
94 <para>This detection strategy avoids considering any particular level as a wrapping one. It allows to process wrapped profiles to which a constant (or even a trend) has been added afterwards.
95 </para>
96 <para>
97 Unwrapping consists in reducing every detected jump and somewhat restoring a continuous slope (initially assumed to be so). At each point, it is performed by applying a Y-shift on a whole side of the profile, with respect to the other. The Y-shift may be the same for all breakpoints, as provided by the user. If the user specifies a null Y-shift, <function>unwrap()</function> applies a jump equal to the average neighbouring slope, depending on each breakpoint.
98 </para>
99 <warning>An unwrapped profile is always defined but for a constant.</warning>
100 <para>
101 Unless <varname>dir</varname> is used, <function>unwrap()</function> applied to a surface unwraps its first column. Each point of this one is then taken as reference level for unwrapping each line starting with it.
102 </para>
103 <para> </para>
104 <para>UNFOLDING</para>
105 <para>
106 If the <varname>"unfold"</varname> keyword is used and a profile -- not a surface -- is provided, the profile is assumed to be folded instead of being wrapped.
107 </para>
108 <para>At a folding point -- or "cusp point" --, the profile is continuous, but its slope is broken: the slope has almost the same absolute value on both sides of the point, but is somewhat opposed from one side to the other.
109 </para>
110 <para>
111 Folding may occur for instance when a non-monotonic periodic function and its inverse are applied to a profile X, like with <varname>Y= acos(cos(X))</varname>. Recovering X from Y is quite more difficult than if it was wrapped. Indeed, some ambiguous situations may exist, like when the profile is tangentially folded on one of its quasi-horizontal sections (if any).
112 </para>
113 <para>When a cusp point is detected, a) one side of the profile starting from it is opposed (upside down), and b) the continuity of the profile and of its slope are preserved and retrieved at the considered point (this may need adding a small jump by the local slope).</para>
114 <warning>The slope on the left edge of the input profile is used as starting reference. The unfolded profile may be upside down with respect to the original true one. In addition, as for unwrapping, it is defined but for a constant.
115 </warning>
116 <para>Known limitations: presently, folded surfaces can't be processed.</para>
117 </refsection>
118 <refsection>
119 <title>Examples</title>
120 <para>Unwrapping or unfolding 1D profiles:</para>
121 <programlisting role="example"><![CDATA[
122// 1D EXAMPLES
123// -----------
124f = scf();
125f.figure_size = [800 1000];
126f.figure_position(2) = 0;
127f.figure_name = "unwrap() & ""unfold""" + _(": 1-D examples ");
128ax = gda();
129ax.y_label.font_size=2;
130drawlater()
131
132// Original 1D profile
133t = linspace(-4,4.2,800);
134alpha = t.^2 + t -1;
135subplot(5,2,1)
136titlepage("unwrap(): unwrap | unfold")
137subplot(5,2,2)
138plot(t,alpha)
139t2 = "$\text{Original profile: } \alpha=t^2+t-1$";
140ax = gca();
141ax.tight_limits = "on";
142yT = max(alpha) - strange(alpha)*0.1;
143xstring(0,yT,t2)
144e = gce();
145e.text_box_mode = "centered";
146e.font_size = 2;
147
148// Loop over cases
149for i=1:4
150 subplot(5,2,2*i+1)
151 if i==1 then
152 // Wrapping by atan(tan())
153 ralpha = atan(tan(alpha)); // Raw recovered alpha [pi]
154 ylabel("$atan(tan(\alpha))$")
155 [u, K] = unwrap(ralpha, %pi); // arctan
156 t2 = "$\text{unwrap(Y, \%pi)}$";
157 elseif i==2
158 // Wrapping by modulo() + Y-shift
159 c = (rand(1,1)-0.5)*4;
160 ralpha = pmodulo(alpha, 5) + c;
161 ylabel("$modulo(\alpha,\ 5)"+msprintf("%+5.2f",c)+"$")
162 [u, K] = unwrap(ralpha, 0);
163 t2 = "$\text{unwrap(Y, 0)}$";
164 elseif i==3
165 // Folding by asin(sin()) + Y-shift
166 ralpha = 1+asin(sin(alpha)); // Raw recovered alpha [2.pi]
167 ylabel("$1+asin(sin(\alpha))$")
168 [u, K] = unwrap(ralpha, "unfold");
169 t2 = "$\text{unwrap(Y,""unfold"")}$";
170 else
171 // Folding by acos(cos()) + a trend
172 ralpha = 1+alpha/10+acos(cos(alpha)); // Raw recovered alpha [2.pi]
173 ylabel("$1+\frac{\alpha}{10}+acos(cos(\alpha))$")
174 [u, K] = unwrap(ralpha, "unfold");
175 t2 = "$\text{unwrap(Y,""unfold"")}$";
176 end
177 // Plotting the profile to be processed
178 plot(t, ralpha)
179 // Staring the breakpoints or the cusp points on the curve:
180 if K~=[] then
181 plot(t(K), ralpha(K),"*")
182 end
183 // Plotting the processed (unwrapped/unfolded) profile:
184 subplot(5,2,2*i+2)
185 plot(t,u)
186 ax = gca();
187 ax.tight_limits = "on";
188 // Adding a legend:
189 yT = max(u) - strange(u)*0.2;
190 xstring(0,yT,t2)
191 e = gce();
192 e.text_box_mode = "centered";
193 e.font_size = 2;
194end
195sda();
196drawnow()
197 ]]></programlisting>
198 <scilab:image>
199 %_unwrap()
200 </scilab:image>
201 <para>Unwrapping 2-D surfaces:</para>
202 <programlisting role="example"><![CDATA[
203// 2-D EXAMPLES
204// -----------
205ax = gda();
206ax.title.font_size = 2;
207f = scf();
208f.color_map = hotcolormap(100);
209f.axes_size = [900 450];
210f.figure_position(2) = 0;
211f.figure_name = "unwrap()" + _(": 2-D examples");
212drawlater()
213
214nx = 300;
215ny = 400;
216rmax = 8.8;
217x = linspace(-rmax/2, rmax/2, nx)-1;
218y = linspace(-rmax/2, rmax/2, ny)+1;
219[X, Y] = meshgrid(x,y);
220for ex=0:1 // examples
221 // Original surface
222 // Generating the surface
223 if ex==0 then
224 z = X.^2 + Y.^2;
225 else
226 z = sqrt(0.3+sinc(sqrt(z)*3))*17-7;
227 end
228 // PLotting it in 3D
229 subplot(2,4,1+4*ex)
230 surf(x, y, z)
231 title("Original profile Z")
232 e = gce();
233 e.thickness = 0; // removes the mesh
234 e.parent.tight_limits = "on";
235
236 // Wrapped surface (flat)
237 m = 2.8;
238 zw = pmodulo(z, m); // wraps it
239 subplot(2,4,2+4*ex)
240 grayplot(x, y, zw')
241 title(msprintf("Zw = pmodulo(Z, %g) (flat)",m))
242 ax0 = gca();
243 ax0.tight_limits = "on";
244
245 // Unwrapped surfaces (flat):
246 // in both directions:
247 u = unwrap(zw, 0);
248 subplot(2,4,3+4*ex)
249 grayplot(x, y, u')
250 title(msprintf("unwrap(Zw, %g) (flat)", 0))
251 ax3 = gca();
252 ax3.tight_limits = "on";
253
254 if ex==0 then
255 direc = "r";
256 else
257 direc = "c";
258 end
259 // Along a single direction:
260 u = unwrap(zw, m, direc);
261 subplot(2,4,4+4*ex)
262 grayplot(x, y, u')
263 title(msprintf("unwrap(Zw, %g, ""%s"") (flat)",m,direc))
264 ax1 = gca();
265 ax1.tight_limits = "on";
266end
267sda();
268drawnow()
269 ]]></programlisting>
270 <para/>
271 <inlinemediaobject>
272 <imageobject>
273 <imagedata fileref="../images/unwrap_2D.png"/>
274 </imageobject>
275 </inlinemediaobject>
276 </refsection>
277 <refsection role="see also">
278 <title>See also</title>
279 <simplelist type="inline">
280 <member>
281 <link linkend="atan">atan</link>
282 </member>
283 <member>
284 <link linkend="modulo">modulo</link>
285 </member>
286 </simplelist>
287 </refsection>
288 <refsection role="bibliography">
289 <title>Bibliography</title>
290 <para>
291 <ulink url="http://siptoolbox.sourceforge.net/doc/sip-0.7.0-reference/unwrapl.html">SIP toolbox: unwrap linear</ulink>
292 </para>
293 <para>
294 <ulink url="http://siptoolbox.sourceforge.net/doc/sip-0.7.0-reference/unwrapp.html">SIP toolbox: unwrap along path</ulink>
295 </para>
296 </refsection>
297 <refsection role="history tag">
298 <title>History</title>
299 <revhistory>
300 <revision>
301 <revnumber>5.5.0</revnumber>
302 <revdescription>unwrap() function introduced</revdescription>
303 </revision>
304 </revhistory>
305 </refsection>
306</refentry>
diff --git a/scilab/modules/signal_processing/help/images/unwrap_2D.png b/scilab/modules/signal_processing/help/images/unwrap_2D.png
new file mode 100644
index 0000000..355acde
--- /dev/null
+++ b/scilab/modules/signal_processing/help/images/unwrap_2D.png
Binary files differ
diff --git a/scilab/modules/signal_processing/help/ja_JP/unwrap.xml b/scilab/modules/signal_processing/help/ja_JP/unwrap.xml
new file mode 100644
index 0000000..0524920
--- /dev/null
+++ b/scilab/modules/signal_processing/help/ja_JP/unwrap.xml
@@ -0,0 +1,389 @@
1<?xml version="1.0" encoding="UTF-8"?>
2<!--
3 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
4 * Copyright (C) 2013 - Samuel GOUGEON
5 *
6 * Copyright (C) 2012 - 2016 - Scilab Enterprises
7 *
8 * This file is hereby licensed under the terms of the GNU GPL v2.0,
9 * pursuant to article 5.3.4 of the CeCILL v.2.1.
10 * This file was originally licensed under the terms of the CeCILL v2.1,
11 * and continues to be available under such terms.
12 * For more information, see the COPYING file which you should have received
13 * along with this program.
14 *
15 -->
16<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns5="http://www.w3.org/1999/xhtml" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:id="unwrap" xml:lang="ja">
17 <refnamediv>
18 <refname>unwrap</refname>
19 <refpurpose>
20 Y(x)輪郭またはZ(x,y)面をアンラップする. Y(x)輪郭を展開する
21 </refpurpose>
22 </refnamediv>
23 <refsynopsisdiv>
24 <title>呼び出し手順</title>
25 <synopsis>unwrap() // 例を実行
26 [U, breakPoints] = unwrap(Y)
27 [U, breakPoints] = unwrap(Y, z_jump)
28 [U, cuspPoints] = unwrap(Y, "unfold")
29 U = unwrap(Z)
30 U = unwrap(Z, z_jump)
31 U = unwrap(Z, z_jump, dir)
32 </synopsis>
33 </refsynopsisdiv>
34 <refsection>
35 <title>引数</title>
36 <variablelist>
37 <varlistentry>
38 <term>Y</term>
39 <listitem>
40 <para>
41 実数ベクトル: アンラップまたは展開する外形.
42 暗黙の座標軸Xが等間隔であることを仮定します.
43 </para>
44 </listitem>
45 </varlistentry>
46 <varlistentry>
47 <term>Z</term>
48 <listitem>
49 <para>
50 実数の行列: アンラップする面. 暗黙の座標軸 (X,Y)
51 は直交軸で等間隔(間隔はXとYで異なることが可能)と仮定されます.
52 </para>
53 </listitem>
54 </varlistentry>
55 <varlistentry>
56 <term>z_jump</term>
57 <listitem>
58 <para>
59 アンラップモードで使用される正のスカラー実数:
60 アンラップを実行する.不連続点で適用されるジャンプ高さ.
61 絶対値のみが考慮されます.
62 実際に適用されるジャンプは各不連続点の
63 両端をつなぐ傾きの符号を有します.
64 デフォルト値は, <literal>z_jump = 2*%pi</literal>です.
65 特別な値 <literal>z_jump = 0</literal> は,
66 各不連続点の周辺の平均的な傾きに
67 等しいジャンプを適用し,
68 外形または面全体に連続的な傾きを回復します.
69 </para>
70 </listitem>
71 </varlistentry>
72 <varlistentry>
73 <term>dir</term>
74 <listitem>
75 <para>
76 "c" | "r" | "" (デフォルト):
77 アンラップを行う方向.
78 "c" は列方向にアンラップ, "r" は行方向にアンラップ,
79 "" は両方向にアンラップします.
80 </para>
81 </listitem>
82 </varlistentry>
83 <varlistentry>
84 <term>"unfold"</term>
85 <listitem>
86 <para>
87 指定した曲線が折り畳まれている場合,アンラップではなく,展開します.
88 </para>
89 </listitem>
90 </varlistentry>
91 <varlistentry>
92 <term>U</term>
93 <listitem>
94 <para>
95 外形または面をアンラップ, または外形を展開します.
96 <varname>U</varname> は<varname>Y</varname> または
97 <varname>Z</varname>と同じ大きさとなります.
98 </para>
99 </listitem>
100 </varlistentry>
101 <varlistentry>
102 <term>breakPoints, cuspPoints</term>
103 <listitem>
104 <para>
105 <varname>Y</varname>の点の添字のベクトル.
106 ただし,ラップ及び折りたたみが検出され,処理されています.
107 </para>
108 </listitem>
109 </varlistentry>
110 </variablelist>
111 </refsection>
112 <refsection>
113 <title>説明</title>
114 <para>アンラップ処理</para>
115 <para>
116 <function>unwrap()</function> は,例えば
117 <literal>Y = modulo(X,w)</literal> または
118 <literal>Y = atan(tan(X))</literal>のような
119 周期的で単調な関数により外形またはラップされた面を
120 処理する際に有用です.
121 この関数は,これらの関数の逆関数のように,
122 制限された
123 <literal>w</literal> または <literal>[-%pi/2, %pi/2]</literal>
124 の範囲ではなく,全領域で入力 <literal>X</literal> を復元する
125 目的で使用されます.
126 </para>
127 <para>
128 ラップされた外形のブレークポイントは,
129 点の近傍の両端の傾きがほぼ同じだが(絶対値が)非常に小さく,
130 考慮される点における傾きとは逆となる点として検出されます:
131 この点において, 近傍の傾きを分断し, 逆となるジャンプがあります.
132 </para>
133 <para>
134 この検出方法はラップするものを特定のレベルで考慮する
135 ことを避けます.
136 これにより,後で定数(またはトレンドすらも)を追加された
137 ラップされた外形を処理できるようになります.
138 </para>
139 <para>
140 アンラップ処理は,全ての検出されたジャンプを減らし,
141 (元が連続だったとすると)一部を連続した傾きに回復します.
142 この処理は,各点における外形のその他の面に関する全側面への
143 Y-シフト適用により行われます.
144 Yシフトは,ユーザにより指定されるため,
145 全てのブレークポイントで同じとなる可能性があります.
146 ユーザがYシフトにヌルを指定した場合,
147 <function>unwrap()</function> は,
148 各ブレークポイントに依存して近傍の平均的傾きに
149 等しいジャンプを適用します.
150 </para>
151 <warning>
152 アンラップされた外形は常に定義されますが,
153 一定ではありません.
154 </warning>
155 <para>
156 <varname>dir</varname> が使用されない限り,
157 <function>unwrap()</function> は面の最初の列を
158 アンラップします.
159 これらのデータの各点がその点から始まる各線をアンラップする
160 基準レベルとして使用されます.
161 </para>
162 <para> </para>
163 <para>展開処理</para>
164 <para>
165 キーワード<varname>"unfold"</varname> が使用され,
166 面ではなく外形が指定された場合,
167 外形はラップされているのではなく折り畳まれていると
168 仮定されます.
169 </para>
170 <para>
171 折りたたみ点 -- または "尖点" -- において,
172 外形は連続ですが,その傾きは不連続になります:
173 傾きは絶対値がその両側の点でほぼ等しくなりますが,
174 両側で逆方向になります.
175 </para>
176 <para>
177 折りたたみは例えば<varname>Y= acos(cos(X))</varname>
178 のように非単調周期関数およびその逆関数が
179 外形Xに適用された場合に発生します.
180 YからXを回復するのはラップされた場合と比べてはるかに困難です.
181 実際,
182 外形が準水平切片の一つの上に接線方向に折り畳まれた場合のように,
183 あいまいな状況が存在する可能性があります.
184 </para>
185 <para>
186 尖点を検出した場合,
187 a) そこから始まる外形の片側が逆方向(反対)となり,
188 かつ,
189 b) 外形と傾きの連続性が保持され,考慮する点
190 (局所的な傾きに基づく小さなジャンプを付加する必要が
191 あるかもしれません)で回復されます.
192 </para>
193 <warning>
194 入力される外形の左端の傾きが最初の基準として
195 使用されます.展開された外形は元の本当の外形に
196 対して上下逆となる可能性があります.
197 加えて,アンラップする場合,定義はされますが
198 一定となります.
199 </warning>
200 <para>既知の制限: 現在, 折り畳まれている面は処理できません.</para>
201 </refsection>
202 <refsection>
203 <title>例</title>
204 <para>1次元の外形をアンラップまたは展開:</para>
205 <programlisting role="example"><![CDATA[
206// 1次元の例
207// -----------
208f = scf();
209f.figure_size = [800 1000];
210f.figure_position(2) = 0;
211f.figure_name = "unwrap() & ""unfold""" + _(": 1-D examples ");
212ax = gda();
213ax.y_label.font_size=2;
214drawlater()
215// 元の1次元の外形
216t = linspace(-4,4.2,800);
217alpha = t.^2 + t -1;
218subplot(5,2,1)
219titlepage("unwrap(): unwrap | unfold")
220subplot(5,2,2)
221plot(t,alpha)
222t2 = "$\text{Original profile: } \alpha=t^2+t-1$";
223ax = gca();
224ax.tight_limits = "on";
225yT = max(alpha) - strange(alpha)*0.1;
226xstring(0,yT,t2)
227e = gce();
228e.text_box_mode = "centered";
229e.font_size = 2;
230// 複数の例をループ
231for i=1:4
232 subplot(5,2,2*i+1)
233 if i==1 then
234 // atan(tan())でアンラップ
235 ralpha = atan(tan(alpha)); // 回復したalpha [pi]
236 ylabel("$atan(tan(\alpha))$")
237 [u, K] = unwrap(ralpha, %pi); // arctan
238 t2 = "$\text{unwrap(Y, \%pi)}$";
239 elseif i==2
240 // modulo() + Y-シフトでアンラップ
241 c = (rand(1,1)-0.5)*4;
242 ralpha = pmodulo(alpha, 5) + c;
243 ylabel("$modulo(\alpha,\ 5)"+msprintf("%+5.2f",c)+"$")
244 [u, K] = unwrap(ralpha, 0);
245 t2 = "$\text{unwrap(Y, 0)}$";
246 elseif i==3
247 // asin(sin()) + Y-シフトで展開
248 ralpha = 1+asin(sin(alpha)); // 回復したalpha [2.pi]
249 ylabel("$1+asin(sin(\alpha))$")
250 [u, K] = unwrap(ralpha, "unfold");
251 t2 = "$\text{unwrap(Y,""unfold"")}$";
252 else
253 // acos(cos()) + トレンド で折りたたみ
254 ralpha = 1+alpha/10+acos(cos(alpha)); // 回復したalpha [2.pi]
255 ylabel("$1+\frac{\alpha}{10}+acos(cos(\alpha))$")
256 [u, K] = unwrap(ralpha, "unfold");
257 t2 = "$\text{unwrap(Y,""unfold"")}$";
258 end
259 // 処理する外形をプロット
260 plot(t, ralpha)
261 // 曲線の不連続点または尖点に星印を付けます:
262 if K~=[] then
263 plot(t(K), ralpha(K),"*")
264 end
265 // 処理(展開またはアンラップ)した外形をプロットします:
266 subplot(5,2,2*i+2)
267 plot(t,u)
268 ax = gca();
269 ax.tight_limits = "on";
270 // Adding a legend:
271 yT = max(u) - strange(u)*0.2;
272 xstring(0,yT,t2)
273 e = gce();
274 e.text_box_mode = "centered";
275 e.font_size = 2;
276end
277sda();
278drawnow()
279 ]]></programlisting>
280 <scilab:image>
281 %_unwrap()
282 </scilab:image>
283 <para>2次元面をアンラップ:</para>
284 <programlisting role="example"><![CDATA[
285// 2次元の例
286// -------
287ax = gda();
288ax.title.font_size = 2;
289f = scf();
290f.color_map = hotcolormap(100);
291f.axes_size = [900 450];
292f.figure_position(2) = 0;
293f.figure_name = "unwrap()" + _(": 2-D examples");
294drawlater()
295
296nx = 300;
297ny = 400;
298rmax = 8.8;
299x = linspace(-rmax/2, rmax/2, nx)-1;
300y = linspace(-rmax/2, rmax/2, ny)+1;
301[X, Y] = meshgrid(x,y);
302for ex=0:1 // 例
303 // 元の面
304 // 面を生成
305 if ex==0 then
306 z = X.^2 + Y.^2;
307 else
308 z = sqrt(0.3+sinc(sqrt(z)*3))*17-7;
309 end
310
311 // 3次元でプロットします
312 subplot(2,4,1+4*ex)
313 surf(x, y, z)
314 title("Original profile Z")
315 e = gce();
316 e.thickness = 0; // メッシュを削除
317 e.parent.tight_limits = "on";
318
319 // ラップされた面 (平面)
320 m = 2.8;
321 zw = pmodulo(z, m); // ラップする
322 subplot(2,4,2+4*ex)
323 grayplot(x, y, zw')
324 title(msprintf("Zw = pmodulo(Z, %g) (flat)",m))
325 ax0 = gca();
326 ax0.tight_limits = "on";
327
328 // アンラップされた面 (平面):
329 // in both directions:
330 u = unwrap(zw, 0);
331 subplot(2,4,3+4*ex)
332 grayplot(x, y, u')
333 title(msprintf("unwrap(Zw, %g) (flat)", 0))
334 ax3 = gca();
335 ax3.tight_limits = "on";
336
337 if ex==0 then
338 direc = "r";
339 else
340 direc = "c";
341 end
342 // 一方向にアンラップ:
343 u = unwrap(zw, m, direc);
344 subplot(2,4,4+4*ex)
345 grayplot(x, y, u')
346 title(msprintf("unwrap(Zw, %g, ""%s"") (flat)",m,direc))
347 ax1 = gca();
348 ax1.tight_limits = "on";
349end
350sda();
351drawnow()
352 ]]></programlisting>
353 <para/>
354 <inlinemediaobject>
355 <imageobject>
356 <imagedata fileref="../images/unwrap_2D.png"/>
357 </imageobject>
358 </inlinemediaobject>
359 </refsection>
360 <refsection role="see also">
361 <title>参照</title>
362 <simplelist type="inline">
363 <member>
364 <link linkend="atan">atan</link>
365 </member>
366 <member>
367 <link linkend="modulo">modulo</link>
368 </member>
369 </simplelist>
370 </refsection>
371 <refsection role="bibliography">
372 <title>参考文献</title>
373 <para>
374 <ulink url="http://siptoolbox.sourceforge.net/doc/sip-0.7.0-reference/unwrapl.html">SIP toolbox: unwrap linear</ulink>
375 </para>
376 <para>
377 <ulink url="http://siptoolbox.sourceforge.net/doc/sip-0.7.0-reference/unwrapp.html">SIP toolbox: unwrap along path</ulink>
378 </para>
379 </refsection>
380 <refsection role="history tag">
381 <title>履歴</title>
382 <revhistory>
383 <revision>
384 <revnumber>5.5.0</revnumber>
385 <revdescription>unwrap() 関数が導入されました</revdescription>
386 </revision>
387 </revhistory>
388 </refsection>
389</refentry>
diff --git a/scilab/modules/signal_processing/macros/%_unwrap.sci b/scilab/modules/signal_processing/macros/%_unwrap.sci
new file mode 100644
index 0000000..a77d71b
--- /dev/null
+++ b/scilab/modules/signal_processing/macros/%_unwrap.sci
@@ -0,0 +1,158 @@
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
13function %_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
158endfunction
diff --git a/scilab/modules/signal_processing/macros/sgolay.sci b/scilab/modules/signal_processing/macros/sgolay.sci
new file mode 100644
index 0000000..b53a029
--- /dev/null
+++ b/scilab/modules/signal_processing/macros/sgolay.sci
@@ -0,0 +1,86 @@
1//This file is part of the Cardiovascular Wawes Analysis toolbox
2//Copyright (C) 2014 - INRIA - Serge Steer
3//This file must be used under the terms of the CeCILL.
4//This source file is licensed as described in the file COPYING, which
5//you should have received as part of this distribution. The terms
6//are also available at
7//http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
8function [B,C] = sgolay(k,nf,w)
9// Savitzky-Golay Filter Design.
10// k : a positive scalar with integer value: the polynomial order
11// nf : a positive scalar with integer value: the FIR filter length,
12// must be odd and greater the k+1
13// w : a real vector of length nf with positive entries: the
14// weights. If omitted no weights are applied.
15
16// B : a real n by n array: the set of filter coefficients: the early rows
17// of B smooth based on future values and later rows smooth based on
18// past values, with the middle row using half future and half past.
19// In particular, you can use row i to estimate x(k) based on the i-1
20// preceding values and the n-i following values of x values as
21// y(k) = B(i,:) * x(k-i+1:k+n-i).
22
23// C : a real n by k+1 array: the matrix of differentiation filters. Each
24// column of C is a differentiation filter for derivatives of order
25// P-1 where P is the column index. Given a signal X with length nf,
26// an estimate of the P-th order derivative of its middle value can be
27// found from:
28// (P)
29// X ((nf+1)/2) = P!*C(:,P+1)'*X
30// References:
31
32// - Abraham Savitzky et Marcel J. E. Golay, « Smoothing and
33// Differentiation of Data by Simplified Least Squares Procedures »,
34// Analytical Chemistry, vol. 8, no 36,‎ 1964, p. 1627–1639 (DOI 10.1021/ac60214a047)
35
36// - http://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter
37
38 if type(k)<>1|int(k)<>k|size(k,'*')<>1 then
39 error(msprintf(_("%s: Wrong value for input argument #%d: An integer value expected.\n"),"sgolay",2))
40 end
41
42 if k > nf-1 then
43 error('The degree must be less than the frame length.'),
44 end
45 if type(k)<>1|int(k)<>k|size(k,'*')<>1 then
46 error(msprintf(_("%s: Wrong type for input argument #%d: a scalar with integer value expected.\n"),...
47 "sgolay",1))
48 end
49
50 if type(nf)<>1|int(nf)<>nf|size(nf,'*')<>1 then
51 error(msprintf(_("%s: Wrong type for input argument #%d: a scalar with integer value expected.\n"),...
52 "sgolay",2))
53 end
54
55 if modulo(nf,2)<>1 then
56 error(msprintf(_("%s: Wrong value for input argument #%d: Must be odd.\n"),"sgolay",2))
57 end
58 if nf<=k then
59 error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Argument #%d expected to be less than argument #%d.\n"),"sgolay",1,2,2,1))
60 end
61
62 //compute the Vandermonde matrix
63 J=(ones(k+1,1)*(-(nf-1)./2:(nf-1)./2)).^((0:k)'*ones(1,nf));
64 if argn(2)==3 then
65 if type(w)<>1|~isreal(w) then
66 error(msprintf(_("%s: Wrong type for argument %d: Real vector expected.\n"),...
67 "sgolay",3))
68 end
69 if size(w,"*")<>nf then
70 error(msprintf(_("%s: Incompatible input arguments #%d and #%d\n"),"sgolay",2,3))
71 end
72 // Check to see if all elements are positive
73 if min(w) <= 0 then
74 error(msprintf(_("%s: Wrong values for input argument #%d: Elements must be positive.\n"),...
75 "sgolay",3))
76
77 end
78 // Diagonalize the vector to form the weighting matrix
79 J = J*diag(sqrt(w))
80 end
81 //Compute the matrix of differentiators C=J'/(J*W*J')
82 C = pinv(J);
83 // Compute the projection matrix B
84 B = C*J;
85
86endfunction
diff --git a/scilab/modules/signal_processing/macros/sgolaydiff.sci b/scilab/modules/signal_processing/macros/sgolaydiff.sci
new file mode 100644
index 0000000..c7225a3
--- /dev/null
+++ b/scilab/modules/signal_processing/macros/sgolaydiff.sci
@@ -0,0 +1,54 @@
1//This file is part of the Cardiovascular Wawes Analysis toolbox
2//Copyright (C) 2014 - INRIA - Serge Steer
3//Copyright (C) 2020 - Stéphane Mottelet (improvement of computation at boudaries)
4//This file must be used under the terms of the CeCILL.
5//This source file is licensed as described in the file COPYING, which
6//you should have received as part of this distribution. The terms
7//are also available at
8//http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
9function D=sgolaydiff(X,order,varargin)
10 if type(X)<>1|~isreal(X) then
11 error(msprintf(_("%s: Wrong type for argument %d: Real matrix expected.\n"),"sgolaydiff",1))
12 end
13 sz=size(X)
14 if size(find(sz==1),'*')==1 then X = X(:);end
15 i1=find(sz>1,1)
16 if i1<>[] then X=matrix(X,sz(i1),-1);end
17
18 if or(size(varargin)==[2 3]) then
19 k=varargin(1);
20 nf=varargin(2);
21 if type(k)<>1|int(k)<>k|size(k,'*')<>1 then
22 error(msprintf(_("%s: Wrong value for input argument #%d: An integer value expected.\n"),"sgolay",2))
23 end
24 if type(nf)<>1|int(nf)<>nf|size(nf,'*')<>1 then
25 error(msprintf(_("%s: Wrong type for input argument #%d: a scalar with integer value expected.\n"),...
26 "sgolay",2))
27 end
28 if modulo(nf,2)==0 then nf=nf+1;varargin(2)=nf;end
29 [B,C] = sgolay(varargin(:))
30 else
31 C=varargin(1)
32 if type(C)<>1|~isreal(C) then
33 error(msprintf(_("%s: Wrong type for argument %d: Real matrix expected.\n"),"sgolaydiff",2))
34 end
35 [nf,n]=size(C);
36 k=n-1;
37 end
38 if size(X,1) < nf then
39 error(msprintf(_("%s: Wrong size for argument %d: At least %d expected.\n"),"sgolaydiff",1,nf))
40 end
41 if order>k then
42 error(msprintf(_("%s: Wrong value for argument %d: At most %d expected.\n"),"sgolaydiff",2,n-1))
43 end
44 nf2=floor(nf/2);
45 //extend X on its boundaries by adding, in reverse order and vertically
46 //symmetrized, copies of the first (nf−1)/2 points at the beginning and
47 //copies of the last (nf−1)/2 points at the end.
48 X = [X(1)-X(nf2+1:-1:2,:);X;2*X($)-X($-1:-1:$-nf2,:)];
49 for k=1:size(X,2)
50 D(:,k)=((-1)^order*order)*conv(X(:,k),C(:,order+1)',"valid");
51 end
52
53 D=matrix(D,sz)
54endfunction
diff --git a/scilab/modules/signal_processing/macros/sgolayfilt.sci b/scilab/modules/signal_processing/macros/sgolayfilt.sci
new file mode 100644
index 0000000..343ea80
--- /dev/null
+++ b/scilab/modules/signal_processing/macros/sgolayfilt.sci
@@ -0,0 +1,66 @@
1//This file is part of the Cardiovascular Wawes Analysis toolbox
2//Copyright (C) 2014 - INRIA - Serge Steer
3//This file must be used under the terms of the CeCILL.
4//This source file is licensed as described in the file COPYING, which
5//you should have received as part of this distribution. The terms
6//are also available at
7//http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
8function y=sgolayfilt(X,varargin)
9//filter signal using Savitzky-Golay Filter.
10// y=sgolayfilt(X,B)
11// y=sgolayfilt(X,k,nf [,w])
12//
13// X :
14// B : a real n by n array: the set of filter coefficients: the early rows
15// of B smooth based on future values and later rows smooth based on
16// past values, with the middle row using half future and half past.
17// In particular, you can use row i to estimate x(k) based on the i-1
18// preceding values and the n-i following values of x values as
19// y(k) = B(i,:) * x(k-i+1:k+n-i).
20
21// k : a positive scalar with integer value: the polynomial order, must
22// be odd
23// nf : a positive scalar with integer value: the FIR filter length,
24// must be odd and greater the k+1
25// w : a real vector of length nf with positive entries: the weights.
26// If omitted no weights are applied.
27//
28// References:
29
30// - Abraham Savitzky et Marcel J. E. Golay, « Smoothing and
31// Differentiation of Data by Simplified Least Squares Procedures »,
32// Analytical Chemistry, vol. 8, no 36,‎ 1964, p. 1627–1639 (DOI 10.1021/ac60214a047)
33
34// - http://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter
35 if type(X)<>1|~isreal(X) then
36 error(msprintf(_("%s: Wrong type for argument %d: Real matrix expected.\n"),...
37 "sgolayfilt",1))
38 end
39 if size(X,1) == 1 then X = X(:);end
40
41 if or(size(varargin)==[2 3]) then
42 B = sgolay(varargin(:))
43 else
44 B=varargin(1)
45 if type(B)<>1|~isreal(B) then
46 error(msprintf(_("%s: Wrong type for argument %d: Real matrix expected.\n"),...
47 "sgolayfilt",2))
48 end
49 end
50 n=size(B,1);
51 if size(X,1) < n,
52 error(msprintf(_("%s: Wrong size for argument %d: At least %d expected.\n"),"sgolayfilt",1,size(B,1)))
53 end
54
55 // The first k rows of B are used to filter the first k points
56 // of the data set based on the first n points of the data set.
57 // The last k rows of B are used to filter the last k points
58 // of the data set based on the last n points of the dataset.
59 // The remaining data is filtered using the central row of B.
60 k = floor(n/2);
61 z=[]
62 for kc=1:size(X,2)
63 z(:,kc) = filter(B(k+1,:), 1, X(:,kc));
64 end
65 y = [B(1:k,:)*X(1:n,:) ; z(n:$,:) ; B(k+2:n,:)*X($-n+1:$,:) ];
66endfunction
diff --git a/scilab/modules/signal_processing/macros/unwrap.sci b/scilab/modules/signal_processing/macros/unwrap.sci
new file mode 100644
index 0000000..9faa852
--- /dev/null
+++ b/scilab/modules/signal_processing/macros/unwrap.sci
@@ -0,0 +1,196 @@
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
13function [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
196endfunction
diff --git a/scilab/modules/signal_processing/tests/nonreg_tests/bug_16305.tst b/scilab/modules/signal_processing/tests/nonreg_tests/bug_16305.tst
new file mode 100644
index 0000000..fd8fdbe
--- /dev/null
+++ b/scilab/modules/signal_processing/tests/nonreg_tests/bug_16305.tst
@@ -0,0 +1,23 @@
1// =============================================================================
2// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3// Copyright (C) 2021 - Samuel GOUGEON
4//
5// This file is distributed under the same license as the Scilab package.
6// =============================================================================
7//
8// <-- CLI SHELL MODE -->
9// <-- NO CHECK REF -->
10// <-- Non-regression test for bug 16305 -->
11//
12// <-- Bugzilla URL -->
13// http://bugzilla.scilab.org/16305
14//
15// <-- Short Description -->
16// unwrap(,"unfold") failed unfolding straight segments
17
18x = linspace(0,10,1000);
19y = acos(cos(x));
20r = unwrap(y,'unfold');
21
22Ref = [0. 1.2412412 2.4824825 3.7190723 4.9603135 6.2015547 7.436594 8.6778353 9.9051221];
23assert_checkalmostequal(r(1:124:1000), Ref, 1e-7);
diff --git a/scilab/modules/signal_processing/tests/unit_tests/unwrap.dia.ref b/scilab/modules/signal_processing/tests/unit_tests/unwrap.dia.ref
new file mode 100644
index 0000000..72cc289
--- /dev/null
+++ b/scilab/modules/signal_processing/tests/unit_tests/unwrap.dia.ref
@@ -0,0 +1,77 @@
1// =============================================================================
2// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3// Copyright (C) 2013 - Samuel Gougeon
4//
5// This file is distributed under the same license as the Scilab package.
6// =============================================================================
7//
8// <-- CLI SHELL MODE -->
9// Example #1 (1D wrapped)
10t = linspace(-4, 4.2, 800);
11Y = t.^2 + t -1;
12Yw = atan(tan(Y)); // Raw wrapped recovered Y [pi]
13Yunw = unwrap(Yw, %pi);
14assert_checkequal(size(Yunw), size(Y));
15assert_checkalmostequal(max(Yunw), 8.273629385640850486539, [], 1e-12);
16assert_checkalmostequal(min(Yunw), -13.816370473381940797708, [], 1e-12);
17// Example #2 (1D wrapped)
18rand("seed", 0);
19c = (rand(1, 1)-0.5)*4;
20Yw = pmodulo(Y, 5) + c; // unwrap() is unsensitive to a constant and even to a trend
21[Yunw, K] = unwrap(Yw, 0);
22assert_checkequal(size(Yunw), size(Y));
23assert_checkalmostequal(max(Yunw), 9.6852994618564913764658, [], 1e-12);
24assert_checkalmostequal(min(Yunw), -12.404700397166287473283, [], 1e-12);
25assert_checkequal(K, [15 98 233 450 585 668 734 791]);
26// Example #3 (1D folded)
27Yf = 1 + asin(sin(Y)); // Raw folded recovered Y [2.pi], + constant Y shift
28Yunf = unwrap(Yf, "unfold");
29assert_checkequal(size(Yunf), size(Y));
30assert_checkalmostequal(max(Yunf), 9.4279967527199310950436, [], 1e-12);
31assert_checkalmostequal(min(Yunf), -12.436144421661456505035, [], 1e-12);
32// Example #4 (1D folded)
33Yf = 1 + Y/10 + acos(cos(Y)); // Raw recovered Y [2.pi], + linear trend
34[Yunf, K] = unwrap(Yf, "unfold");
35assert_checkequal(size(Yunf), size(Y));
36assert_checkalmostequal(max(Yunf), 15.801371238872524926933, [], 1e-12);
37assert_checkalmostequal(min(Yunf), -6.1759877188984493301405, [], 1e-12);
38assert_checkequal(K, [24 75 138 233 451 546 609 660 704 743 779]);
39//Unwrapping 2D surfaces:
40nx = 300;
41ny = 400;
42rmax = 8.8;
43x = linspace(-rmax/2, rmax/2, nx)-1;
44y = linspace(-rmax/2, rmax/2, ny)+1;
45[X, Y] = meshgrid(x, y);
46// Example #5 : 2D wrapped
47Z = X.^2 + Y.^2; // Paraboloïd
48m = 2.8;
49Zw = pmodulo(Z, m); // >raps it
50Zunw = unwrap(Zw, 0); // Unwraps it in both directions, with free jumps
51assert_checkequal(size(Zunw), size(Z));
52assert_checkalmostequal(max(Zunw), 19.120000000000032969183, [], 1e-12);
53assert_checkalmostequal(min(Zunw), -39.199790375290433531, [], 1e-12);
54ZunwR = unwrap(Zw, 0, "r"); // Unwraps it along rows
55assert_checkequal(size(ZunwR), size(Z));
56assert_checkalmostequal(max(ZunwR), 2.78738299382542109583, [], 1e-12);
57assert_checkalmostequal(min(ZunwR), -28.995103166214036605197, [], 1e-12);
58ZunwC = unwrap(Zw, 0, "c"); // Unwraps it along columns
59assert_checkequal(size(ZunwC), size(Z));
60assert_checkalmostequal(max(ZunwC), 20.398076084160148724322, [], 1e-12);
61assert_checkalmostequal(min(ZunwC), -11.462196183827872530969, [], 1e-12);
62// Example #6 : 2D wrapped
63Z = sqrt(0.3+sinc(sqrt(Z)*3))*17-7; // Defines the surface
64m = 2.8; // Wrapping step (jump's amplitude)
65Zw = pmodulo(Z, m); // Wraps it
66Zunw = unwrap(Zw, m); // Unwraps it in both directions, with given jump
67assert_checkequal(size(Zunw), size(Z));
68assert_checkalmostequal(max(Zunw), 12.380638179897072603808, [], 1e-12);
69assert_checkalmostequal(min(Zunw), -2.109245308884592162713, [], 1e-12);
70ZunwR = unwrap(Zw, m, "r"); // Unwraps it along rows
71assert_checkequal(size(ZunwR), size(Z));
72assert_checkalmostequal(max(ZunwR), 12.380638179897072603808, [], 1e-12);
73assert_checkalmostequal(min(ZunwR), -2.109245308884592162713, [], 1e-12);
74ZunwC = unwrap(Zw, m, "c"); // Unwraps it along columns
75assert_checkequal(size(ZunwC), size(Z));
76assert_checkalmostequal(max(ZunwC), 12.380638179897072603808, [], 1e-12);
77assert_checkalmostequal(min(ZunwC), -2.1092453088845903863557, [], 1e-12);
diff --git a/scilab/modules/signal_processing/tests/unit_tests/unwrap.tst b/scilab/modules/signal_processing/tests/unit_tests/unwrap.tst
new file mode 100644
index 0000000..7391fe9
--- /dev/null
+++ b/scilab/modules/signal_processing/tests/unit_tests/unwrap.tst
@@ -0,0 +1,84 @@
1// =============================================================================
2// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3// Copyright (C) 2013 - Samuel Gougeon
4//
5// This file is distributed under the same license as the Scilab package.
6// =============================================================================
7//
8// <-- CLI SHELL MODE -->
9
10// Example #1 (1D wrapped)
11t = linspace(-4, 4.2, 800);
12Y = t.^2 + t -1;
13Yw = atan(tan(Y)); // Raw wrapped recovered Y [pi]
14Yunw = unwrap(Yw, %pi);
15assert_checkequal(size(Yunw), size(Y));
16assert_checkalmostequal(max(Yunw), 8.273629385640850486539, [], 1e-12);
17assert_checkalmostequal(min(Yunw), -13.816370473381940797708, [], 1e-12);
18
19// Example #2 (1D wrapped)
20rand("seed", 0);
21c = (rand(1, 1)-0.5)*4;
22Yw = pmodulo(Y, 5) + c; // unwrap() is unsensitive to a constant and even to a trend
23[Yunw, K] = unwrap(Yw, 0);
24assert_checkequal(size(Yunw), size(Y));
25assert_checkalmostequal(max(Yunw), 9.6852994618564913764658, [], 1e-12);
26assert_checkalmostequal(min(Yunw), -12.404700397166287473283, [], 1e-12);
27assert_checkequal(K, [15 98 233 450 585 668 734 791]);
28
29// Example #3 (1D folded)
30Yf = 1 + asin(sin(Y)); // Raw folded recovered Y [2.pi], + constant Y shift
31Yunf = unwrap(Yf, "unfold");
32assert_checkequal(size(Yunf), size(Y));
33assert_checkalmostequal(max(Yunf), 9.4279967527199310950436, [], 1e-12);
34assert_checkalmostequal(min(Yunf), -12.436144421661456505035, [], 1e-12);
35
36// Example #4 (1D folded)
37Yf = 1 + Y/10 + acos(cos(Y)); // Raw recovered Y [2.pi], + linear trend
38[Yunf, K] = unwrap(Yf, "unfold");
39assert_checkequal(size(Yunf), size(Y));
40assert_checkalmostequal(max(Yunf), 15.801371238872524926933, [], 1e-12);
41assert_checkalmostequal(min(Yunf), -6.1759877188984493301405, [], 1e-12);
42assert_checkequal(K, [24 75 138 233 451 546 609 660 704 743 779]);
43
44//Unwrapping 2D surfaces:
45nx = 300;
46ny = 400;
47rmax = 8.8;
48x = linspace(-rmax/2, rmax/2, nx)-1;
49y = linspace(-rmax/2, rmax/2, ny)+1;
50[X, Y] = meshgrid(x, y);
51
52// Example #5 : 2D wrapped
53Z = X.^2 + Y.^2; // Paraboloïd
54m = 2.8;
55Zw = pmodulo(Z, m); // >raps it
56Zunw = unwrap(Zw, 0); // Unwraps it in both directions, with free jumps
57assert_checkequal(size(Zunw), size(Z));
58assert_checkalmostequal(max(Zunw), 19.120000000000032969183, [], 1e-12);
59assert_checkalmostequal(min(Zunw), -39.199790375290433531, [], 1e-12);
60ZunwR = unwrap(Zw, 0, "r"); // Unwraps it along rows
61assert_checkequal(size(ZunwR), size(Z));
62assert_checkalmostequal(max(ZunwR), 2.78738299382542109583, [], 1e-12);
63assert_checkalmostequal(min(ZunwR), -28.995103166214036605197, [], 1e-12);
64ZunwC = unwrap(Zw, 0, "c"); // Unwraps it along columns
65assert_checkequal(size(ZunwC), size(Z));
66assert_checkalmostequal(max(ZunwC), 20.398076084160148724322, [], 1e-12);
67assert_checkalmostequal(min(ZunwC), -11.462196183827872530969, [], 1e-12);
68
69// Example #6 : 2D wrapped
70Z = sqrt(0.3+sinc(sqrt(Z)*3))*17-7; // Defines the surface
71m = 2.8; // Wrapping step (jump's amplitude)
72Zw = pmodulo(Z, m); // Wraps it
73Zunw = unwrap(Zw, m); // Unwraps it in both directions, with given jump
74assert_checkequal(size(Zunw), size(Z));
75assert_checkalmostequal(max(Zunw), 12.380638179897072603808, [], 1e-12);
76assert_checkalmostequal(min(Zunw), -2.109245308884592162713, [], 1e-12);
77ZunwR = unwrap(Zw, m, "r"); // Unwraps it along rows
78assert_checkequal(size(ZunwR), size(Z));
79assert_checkalmostequal(max(ZunwR), 12.380638179897072603808, [], 1e-12);
80assert_checkalmostequal(min(ZunwR), -2.109245308884592162713, [], 1e-12);
81ZunwC = unwrap(Zw, m, "c"); // Unwraps it along columns
82assert_checkequal(size(ZunwC), size(Z));
83assert_checkalmostequal(max(ZunwC), 12.380638179897072603808, [], 1e-12);
84assert_checkalmostequal(min(ZunwC), -2.1092453088845903863557, [], 1e-12);