[signal_processing] add Savitsky-Golay filtering macros 99/21499/10
Serge Steer [Fri, 12 Jun 2020 13:31:43 +0000 (15:31 +0200)]
Courtesy of S. Steer. These filters were too usefull to be kept in
Cardiovascular Wave Analysis toolbox.

sgolaydiff was improved at boundaries (odd mirroring is better than
even mirroring for approximating derivative).

Change-Id: I0c4d46d2c18fb3e969db0d8bc924a9f0d9fe92af

14 files changed:
scilab/CHANGES.md
scilab/modules/differential_equations/help/en_US/diff.xml
scilab/modules/helptools/etc/images_md5.txt
scilab/modules/helptools/images/sgolay_1.png [new file with mode: 0644]
scilab/modules/helptools/images/sgolay_2.png [new file with mode: 0644]
scilab/modules/helptools/images/sgolaydiff_1.png [new file with mode: 0644]
scilab/modules/helptools/images/sgolayfilt_1.png [new file with mode: 0644]
scilab/modules/signal_processing/help/en_US/filters/filter.xml
scilab/modules/signal_processing/help/en_US/filters/sgolay.xml [new file with mode: 0755]
scilab/modules/signal_processing/help/en_US/filters/sgolaydiff.xml [new file with mode: 0755]
scilab/modules/signal_processing/help/en_US/filters/sgolayfilt.xml [new file with mode: 0755]
scilab/modules/signal_processing/macros/sgolay.sci [new file with mode: 0644]
scilab/modules/signal_processing/macros/sgolaydiff.sci [new file with mode: 0644]
scilab/modules/signal_processing/macros/sgolayfilt.sci [new file with mode: 0644]

index f40be82..053255d 100644 (file)
@@ -213,6 +213,7 @@ Feature changes and additions on 6.1.1
   - Decimal positive integers > 2^52 up to 2^1024 = number_properties("huge") can now be processed, with bitnum up to 1024 instead of 52.
   - bitnum can actually be an array. It is now optional as well for input decimal integers.
 * The `Arnoldi` module is now internal.
+* `sgolay` and the companion `sgolayfilter` and `sgolaydiff` functions have been added to implement Savitsky-Golay filters.
 
 Help pages:
 -----------
index 1002f90..03566f0 100644 (file)
@@ -119,6 +119,9 @@ norm(dy-cos(t(1:$-1)),%inf)
             <member>
                 <link linkend="numderivative">numderivative</link>
             </member>
+            <member>
+                <link linkend="sgolaydiff">sgolaydiff</link>
+            </member>
         </simplelist>
     </refsection>
 
index ecfe3fd..e7a0587 100644 (file)
@@ -1022,6 +1022,10 @@ scatter_9.png=9f987520c26d42c41f2ee6d3ec930529
 scifunc_block_m_1.png=a26368cd2a69c472967105b9d8deb5ed
 sec_1.png=4eff2516ae7932bb126b6e2fa133f046
 sech_1.png=62b827923764aaaf8ce355b0256bef3b
+sgolay_1.png=a68d1b4aa9bbf7ccfcd38dfd6158cf01
+sgolay_2.png=9ec7d521df8c827f7ef37e2997613bd6
+sgolaydiff_1.png=68aacc1223db752ac39fbede1935ac75
+sgolayfilt_1.png=e6c2afd9293c4a59cddb04e2a4dd8ce5
 sgrid_1.png=a5e9d0fb3e290982e4c01a50b800b8ff
 sgrid_2.png=260f87be847f70d84df4d1013bd9113b
 sgrid_3.png=b26fd988424a754a724bd483a4f26e0f
diff --git a/scilab/modules/helptools/images/sgolay_1.png b/scilab/modules/helptools/images/sgolay_1.png
new file mode 100644 (file)
index 0000000..1fb6b91
Binary files /dev/null and b/scilab/modules/helptools/images/sgolay_1.png differ
diff --git a/scilab/modules/helptools/images/sgolay_2.png b/scilab/modules/helptools/images/sgolay_2.png
new file mode 100644 (file)
index 0000000..3b7e5f5
Binary files /dev/null and b/scilab/modules/helptools/images/sgolay_2.png differ
diff --git a/scilab/modules/helptools/images/sgolaydiff_1.png b/scilab/modules/helptools/images/sgolaydiff_1.png
new file mode 100644 (file)
index 0000000..23e02e3
Binary files /dev/null and b/scilab/modules/helptools/images/sgolaydiff_1.png differ
diff --git a/scilab/modules/helptools/images/sgolayfilt_1.png b/scilab/modules/helptools/images/sgolayfilt_1.png
new file mode 100644 (file)
index 0000000..012a785
Binary files /dev/null and b/scilab/modules/helptools/images/sgolayfilt_1.png differ
index a35eebe..6c7d479 100644 (file)
             <member>
                 <link linkend="ltitr">ltitr</link>
             </member>
+            <member>
+                <link linkend="sgolayfilt">sgolayfilt</link>
+            </member>
         </simplelist>
     </refsection>
 </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 (executable)
index 0000000..e6d742a
--- /dev/null
@@ -0,0 +1,259 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!--
+ This file is part of the Cardiovascular Wawes Analysis toolbox
+ Copyright (C) 2014 - INRIA - Serge Steer
+ This file must be used under the terms of the CeCILL.
+ This source file is licensed as described in the file COPYING, which
+ you should have received as part of this distribution.  The terms
+ are also available at
+ http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+-->
+<refentry xml:id="sgolay" xml:lang="en" xmlns="http://docbook.org/ns/docbook"
+          xmlns:xlink="http://www.w3.org/1999/xlink"
+          xmlns:svg="http://www.w3.org/2000/svg"
+          xmlns:scilab="http://www.scilab.org"
+          xmlns:ns4="http://www.w3.org/1999/xhtml"
+          xmlns:mml="http://www.w3.org/1998/Math/MathML"
+          xmlns:db="http://docbook.org/ns/docbook">
+  <refnamediv>
+    <refname>sgolay</refname>
+
+    <refpurpose>Savitzky-Golay Filter Design.</refpurpose>
+  </refnamediv>
+
+  <refsynopsisdiv>
+    <title>Calling Sequence</title>
+
+    <synopsis>[B [,C]] = sgolay(k,nf [,w])</synopsis>
+  </refsynopsisdiv>
+
+  <refsection>
+    <title>Arguments</title>
+
+    <variablelist>
+      <varlistentry>
+        <term>k</term>
+
+        <listitem>
+          <para>a positive scalar with integer value: the fitting polynomial
+          degree.</para>
+        </listitem>
+      </varlistentry>
+
+      <varlistentry>
+        <term>nf</term>
+
+        <listitem>
+          <para>a positive scalar with integer value: the filter length, must
+          be odd and greater the k+1.</para>
+        </listitem>
+      </varlistentry>
+
+      <varlistentry>
+        <term>w</term>
+
+        <listitem>
+          <para>a real vector of length nf with positive entries: the weights.
+          If omitted no weights are applied.</para>
+        </listitem>
+      </varlistentry>
+
+      <varlistentry>
+        <term>B</term>
+
+        <listitem>
+          <para>a real n by n array: the set of filter coefficients: the early
+          rows of B smooth based on future values and later rows smooth based
+          on past values, with the middle row using half future and half past.
+          In particular, you can use row i to estimate x(k) based on the i-1
+          preceding values and the n-i following values of x values as y(k) =
+          B(i,:) * x(k-i+1:k+n-i).</para>
+        </listitem>
+      </varlistentry>
+
+      <varlistentry>
+        <term>C</term>
+
+        <listitem>
+          <para>a real n by k+1 array: the matrix of differentiation filters.
+          Each column of C is a differentiation filter for derivatives of
+          order P-1 where P is the column index. Given a signal X with length
+          nf, an estimate of the P-th order derivative of its middle value can
+          be found from: <para/> (P) X ((nf+1)/2) = P!*C(:,P+1)'*X</para>
+        </listitem>
+      </varlistentry>
+    </variablelist>
+  </refsection>
+
+  <refsection>
+    <title>Description</title>
+
+    <para>This function computes the smoothing FIR Savitzky-Golay
+    Filter. This is achieved by fitting successive sub-sets of
+    adjacent data points with a low-degree polynomial by the method of
+    linear least squares.</para>
+
+    <para>This filter can also be used to approximate numerical
+    derivatives.</para>
+  </refsection>
+
+  <refsection>
+    <title>Examples</title>
+
+    <variablelist>
+      <varlistentry>
+        <listitem>
+          <para>The sgolay function can be used to construct FIR filter with
+          no group delay. The following instructions explains what the <link
+          linkend="sgolayfilt">sgolayfilt</link> function does:</para>
+
+          <programlisting role="example"><![CDATA[
+dt=0.01;
+t = (0:0.01:4*%pi)';
+x = sin(t)+0.05*rand(t,'normal');
+
+nf = 61;
+[B,C] = sgolay(3,nf); 
+nfs2 = floor(nf/2);
+
+//compute the middle part of the filtered signal
+xfm = filter(B(nfs2+1,:), 1, x);
+xfm = xfm(nf:$);
+//compute the left part of the filtered signal
+xfl = B(1:nfs2,:)*x(1:nf);
+//compute the right part of the filtered signal
+xfr = B(nfs2+2:nf,:)*x($-nf+1:$);
+clf;ax=gca();
+plot(t,x,'b',t,[xfl;xfm;xfr],'r');
+e=gce();set(e.children(1),"thickness",2);
+legend(["Raw","Filtered"]);
+    ]]></programlisting>
+
+          <para>The above script generate this graph:</para>
+
+          <scilab:image><![CDATA[
+dt=0.01;
+t = (0:0.01:4*%pi)';
+x = sin(t)+0.05*rand(t,'normal');
+
+nf = 61;
+[B,C] = sgolay(3,nf); 
+nfs2 = floor(nf/2);
+
+//compute the middle part of the filtered signal
+xfm = filter(B(nfs2+1,:), 1, x);
+xfm = xfm(nf:$);
+//compute the left part of the filtered signal
+xfl = B(1:nfs2,:)*x(1:nf);
+//compute the right part of the filtered signal
+xfr = B(nfs2+2:nf,:)*x($-nf+1:$);
+clf;ax=gca();
+plot(t,x,'b',t,[xfl;xfm;xfr],'r');
+e=gce();set(e.children(1),"thickness",2);
+legend(["Raw","Filtered"]);
+    ]]></scilab:image>
+        </listitem>
+      </varlistentry>
+
+      <varlistentry>
+        <listitem>
+          <para>It can also be used to obtain approximate derivatives:</para>
+
+          <programlisting role="example"><![CDATA[
+dt=0.01;
+t = (0:0.01:4*%pi)';
+x = sin(t)+0.03*rand(t,'normal');
+
+nf = 61;
+[B,C] = sgolay(3,nf); 
+nfs2 = floor(nf/2);
+
+dx=-conv(x,C(:,2)',"valid");
+d2x=2*conv(x,C(:,3)',"valid");
+
+clf();
+subplot(211)
+plot(t(2:$),diff(x)/dt,'g',t,cos(t),'b',t(nfs2+1:$-nfs2),dx/dt,'r');
+legend(["diff","theoretical","sgolay"]);
+title("First order differentiation")
+subplot(212)
+plot(t,-sin(t),'b',t(nfs2+1:$-nfs2),d2x/(dt^2),'r');
+legend(["theoretical","sgolay"]);
+title("Second order differentiation")
+    ]]></programlisting>
+
+          <scilab:image><![CDATA[
+dt=0.01;
+t = (0:0.01:4*%pi)';
+x = sin(t)+0.03*rand(t,'normal');
+
+nf = 61;
+[B,C] = sgolay(3,nf); 
+nfs2 = floor(nf/2);
+
+dx=-conv(x,C(:,2)',"valid");
+d2x=2*conv(x,C(:,3)',"valid");
+
+clf();
+subplot(211)
+plot(t(2:$),diff(x)/dt,'g',t,cos(t),'b',t(nfs2+1:$-nfs2),dx/dt,'r');
+legend(["diff","theoretical","sgolay"]);
+title("First order differentiation")
+subplot(212)
+plot(t,-sin(t),'b',t(nfs2+1:$-nfs2),d2x/(dt^2),'r');
+legend(["theoretical","sgolay"]);
+title("Second order differentiation")
+    ]]></scilab:image>
+        </listitem>
+      </varlistentry>
+    </variablelist>
+  </refsection>
+
+  <refsection>
+    <title>See Also</title>
+
+    <simplelist type="inline">
+      <member><link linkend="sgolayfilt">sgolayfilt</link></member>
+
+      <member><link linkend="conv">conv</link></member>
+    </simplelist>
+  </refsection>
+
+  <refsection>
+    <title>Authors</title>
+
+    <simplelist type="vert">
+      <member>Serge Steer, INRIA</member>
+    </simplelist>
+  </refsection>
+
+  <refsection>
+    <title>Bibliography</title>
+
+    <para>Abraham Savitzky et Marcel J. E. Golay, « Smoothing and
+    Differentiation of Data by Simplified Least Squares Procedures »,
+    Analytical Chemistry, vol. 8, no 36,‎ 1964, p. 1627–1639 (DOI
+    10.1021/ac60214a047)</para>
+
+    <para><ulink
+    url="http://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter">http://en.wikipedia.org/wiki/Savitzky-Golay_filter</ulink>.</para>
+  </refsection>
+
+  <refsection>
+    <title>History</title>
+
+    <revhistory>
+      <revision>
+        <revnumber>0.0</revnumber>
+
+        <revdescription>Function added</revdescription>
+      </revision>
+    </revhistory>
+  </refsection>
+
+  <refsection>
+    <title>Used Functions</title>
+
+    <para><link linkend="pinv">pinv</link></para>
+  </refsection>
+</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 (executable)
index 0000000..c50bcb0
--- /dev/null
@@ -0,0 +1,188 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!--
+ Copyright (C) 2014 - INRIA - Serge Steer
+ This file must be used under the terms of the CeCILL.
+ This source file is licensed as described in the file COPYING, which
+ you should have received as part of this distribution.  The terms
+ are also available at
+ http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+
+-->
+<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">
+  <refnamediv>
+    <refname>sgolaydiff</refname>
+    <refpurpose>Numerical derivatives computation using Savitzky-Golay filter.</refpurpose>
+  </refnamediv>
+  <refsynopsisdiv>
+    <title>Calling Sequence</title>
+    <synopsis>D = sgolaydiff(X,order,k,nf [,w])</synopsis>
+    <synopsis>D = sgolaydiff(X,order,C)</synopsis>
+  </refsynopsisdiv>
+  <refsection>
+    <title>Arguments</title>
+    <variablelist>
+      <varlistentry>
+        <term>X</term>
+        <listitem>
+          <para>
+            It can be either a real vector or a general real array. If
+            it is an array the derivation is applied along the first
+            non singleton dimension.
+          </para>
+        </listitem>
+      </varlistentry>
+      <varlistentry>
+        <term>order</term>
+        <listitem>
+          <para>
+            a positive scalar with integer value, the order of derivation.
+          </para>
+        </listitem>
+      </varlistentry>
+      <varlistentry>
+        <term>C</term>
+        <listitem>
+          <para>a real nf by k+1 array: the matrix of differentiation
+          filters computed by the <link linkend="sgolay">sgolay</link> function.</para>
+        </listitem>
+      </varlistentry>
+     <varlistentry>
+        <term>k</term>
+        <listitem>
+          <para>a positive scalar with integer value: the fitting polynomial
+          degree.</para>
+        </listitem>
+      </varlistentry>
+      <varlistentry>
+        <term>nf</term>
+        <listitem>
+          <para>a positive scalar with integer value: the filter length, must
+          be odd and greater the k+1.</para>
+        </listitem>
+      </varlistentry>
+      <varlistentry>
+        <term>w</term>
+
+        <listitem>
+          <para>a real vector of length nf with positive entries: the weights.
+          If omitted no weights are applied.</para>
+        </listitem>
+      </varlistentry>
+      <varlistentry>
+        <term>D</term>
+        <listitem>
+          <para>A vector with identical shape as X, the estimated derivative.</para>
+        </listitem>
+      </varlistentry>
+    </variablelist>
+  </refsection>
+  <refsection>
+    <title>Description</title>
+    <para>
+      <para>This function computes numerical derivatives using the FIR
+      Savitzky-Golay smoothing Filter. This is achieved by fitting
+      successive sub-sets of adjacent data points with a low-degree
+      polynomial by the method of linear least squares.</para>
+    </para>
+  </refsection>
+  <refsection>
+    <title>More information</title>
+    <caution><para>order must be less than k</para></caution>
+    <note><para>The (nf−1)/2 first and last derivative values are
+    estimated by adding to the input signal, in reverse order and vertically symmetrized, copies of the first
+    (nf−1)/2 points at the beginning and copies of the last (nf−1)/2
+    points at the end .</para></note>
+  </refsection>
+  <refsection>
+    <title>Examples</title>
+    <programlisting role="example"><![CDATA[
+//generate a noisy signal
+dt=0.01;
+t = (0:0.01:4*%pi)';
+x = sin(t)+0.03*rand(t,"normal");
+
+clf();
+//first order derivative
+subplot(211)
+dx=sgolaydiff(x,1,3,61);
+plot(t(2:$),diff(x)/dt,'g',t,cos(t),'b',t,dx/dt,'r');
+legend(["diff","theoretical","sgolaydiff"]);
+//second order derivative
+subplot(212)
+d2x=sgolaydiff(x,2,3,101);
+plot(t,-sin(t),'b',t,d2x/(dt^2),'r');
+legend(["theoretical","sgolaydiff"]);
+title("Second order differentiation")
+    ]]></programlisting>
+    <scilab:image><![CDATA[
+//generate a noisy signal
+dt=0.01;
+t = (0:0.01:4*%pi)';
+x = sin(t)+0.03*rand(t,"normal");
+
+clf();
+//first order derivative
+subplot(211)
+dx=sgolaydiff(x,1,3,61);
+plot(t(2:$),diff(x)/dt,'g',t,cos(t),'b',t,dx/dt,'r');
+legend(["diff","theoretical","sgolaydiff"]);
+//second order derivative
+subplot(212)
+d2x=sgolaydiff(x,2,3,101);
+plot(t,-sin(t),'b',t,d2x/(dt^2),'r');
+legend(["theoretical","sgolaydiff"]);
+title("Second order differentiation")
+    ]]></scilab:image>
+  </refsection>
+  <refsection>
+    <title>See Also</title>
+    <simplelist type="inline">
+      <member>
+        <link linkend="sgolay">sgolay</link>
+      </member>
+      <member>
+        <link linkend="sgolayfilt">sgolayfilt</link>
+      </member>
+      <member>
+        <link linkend="diff">diff</link> 
+      </member>
+      <member>
+        <link linkend="numdiff">numdiff</link> 
+      </member>
+    </simplelist>
+  </refsection>
+  <refsection>
+    <title>Authors</title>
+    <simplelist type="vert">
+      <member>Serge Steer, INRIA</member>
+    </simplelist>
+  </refsection>
+  <refsection>
+     <title>Bibliography</title>
+    <para>Abraham Savitzky et Marcel J. E. Golay, « Smoothing and
+    Differentiation of Data by Simplified Least Squares Procedures »,
+    Analytical Chemistry, vol. 8, no 36,‎ 1964, p. 1627–1639 (DOI
+    10.1021/ac60214a047)</para>
+
+    <para><ulink
+    url="http://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter">http://en.wikipedia.org/wiki/Savitzky-Golay_filter</ulink>.</para>
+  </refsection>
+
+    <refsection>
+       <title>History</title>
+      <revhistory>
+        <revision>
+          <revnumber>0.0</revnumber>
+          <revdescription>Function added</revdescription>
+        </revision>
+      </revhistory>
+    </refsection>
+  <refsection>
+     <title>Used Functions</title>
+       <para> This function is based on the 
+         <link linkend="sgolay">sgolay</link> function which computes the   matrix of differentiation
+          filters and on the
+         <link linkend="conv">conv</link> function to compute the smoothed derivative/
+       </para>
+     </refsection>
+</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 (executable)
index 0000000..15ea486
--- /dev/null
@@ -0,0 +1,180 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<!--
+ This file is part of the Cardiovascular Wawes Analysis toolbox
+ Copyright (C) 2014 - INRIA - Serge Steer
+ This file must be used under the terms of the CeCILL.
+ This source file is licensed as described in the file COPYING, which
+ you should have received as part of this distribution.  The terms
+ are also available at
+ http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+-->
+<refentry xml:id="sgolayfilt" xml:lang="en"
+          xmlns="http://docbook.org/ns/docbook"
+          xmlns:xlink="http://www.w3.org/1999/xlink"
+          xmlns:svg="http://www.w3.org/2000/svg"
+          xmlns:scilab="http://www.scilab.org"
+          xmlns:ns4="http://www.w3.org/1999/xhtml"
+          xmlns:mml="http://www.w3.org/1998/Math/MathML"
+          xmlns:db="http://docbook.org/ns/docbook">
+  <refnamediv>
+    <refname>sgolayfilt</refname>
+
+    <refpurpose>Filter signal using Savitzky-Golay Filter.</refpurpose>
+  </refnamediv>
+
+  <refsynopsisdiv>
+    <title>Calling Sequence</title>
+
+    <synopsis>xf=sgolayfilt(X,B)</synopsis>
+
+    <synopsis>xf=sgolayfilt(X,k,nf [,w]) </synopsis>
+  </refsynopsisdiv>
+
+  <refsection>
+    <title>Arguments</title>
+
+    <variablelist>
+      <varlistentry>
+        <term>X</term>
+
+        <listitem>
+          <para>a real array. If it is a vector it contains the signal to be
+          filtered else each column of X is filtered.</para>
+        </listitem>
+      </varlistentry>
+
+      <varlistentry>
+        <term>B</term>
+
+        <listitem>
+          <para>a real n by n array: the set of filter coefficients produced
+          by the <link linkend="sgolay">sgolay</link> function.</para>
+        </listitem>
+      </varlistentry>
+
+      <varlistentry>
+        <term>k</term>
+
+        <listitem>
+          <para>a positive scalar with integer value: the fitting polynomial
+          degree.</para>
+        </listitem>
+      </varlistentry>
+
+      <varlistentry>
+        <term>k</term>
+
+        <listitem>
+          <para>a positive scalar with integer value: the filter length, must
+          be odd and greater the k+1.</para>
+        </listitem>
+      </varlistentry>
+
+      <varlistentry>
+        <term>xf</term>
+
+        <listitem>
+          <para>a real array: the filtered signal(s).</para>
+        </listitem>
+      </varlistentry>
+    </variablelist>
+  </refsection>
+
+  <refsection>
+    <title>Description</title>
+
+    <para>This function applies a sgolay filter which can be specified either
+    by its set of filter coefficients as produced by <link
+    linkend="sgolay">sgolay</link> either by the degree of the polynomial and
+    the filter length.</para>
+   <para>According to the <ulink
+    url="http://www.hpl.hp.com/techreports/2010/HPL-2010-109.pdf">article</ulink> the -3dB cutoff frequency can be approximated by
+   (k+1)/(3.2*nf-4.6).</para>
+  </refsection>
+
+  <refsection>
+    <title>More information</title>
+
+    <warning>
+      <para>Signal size should be greater than the filter length</para>
+    </warning>
+  </refsection>
+
+  <refsection>
+    <title>Examples</title>
+
+    <programlisting role="example"><![CDATA[
+dt=0.01;
+t = (0:0.01:4*%pi)';
+x = sin(t)+0.05*rand(t,'normal');
+clf;plot(t,x,'b',t,sgolayfilt(x,3,51),'r');
+e=gce();set(e.children(1),"thickness",2);
+legend(["Raw","Filtered"]);
+    ]]></programlisting>
+
+    <scilab:image><![CDATA[
+dt=0.01;
+t = (0:0.01:4*%pi)';
+x = sin(t)+0.05*rand(t,'normal');
+clf;plot(t,x,'b',t,sgolayfilt(x,3,51),'r');
+e=gce();set(e.children(1),"thickness",2);
+legend(["Raw","Filtered"]);
+    ]]></scilab:image>
+  </refsection>
+
+  <refsection>
+    <title>See Also</title>
+
+    <simplelist type="inline">
+      <member><link linkend="sgolay">sgolay</link></member>
+
+      <member><link linkend="filter">filter</link></member>
+
+      <member><link linkend="wfir">wfir</link></member>
+    </simplelist>
+  </refsection>
+
+  <refsection>
+    <title>Authors</title>
+
+    <simplelist type="vert">
+      <member>Serge Steer, INRIA</member>
+    </simplelist>
+  </refsection>
+
+  <refsection>
+    <title>Bibliography</title>
+
+    <para>Abraham Savitzky et Marcel J. E. Golay, « Smoothing and
+    Differentiation of Data by Simplified Least Squares Procedures »,
+    Analytical Chemistry, vol. 8, no 36,‎ 1964, p. 1627–1639 (DOI
+    10.1021/ac60214a047)</para>
+
+    <para><ulink
+    url="http://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter">http://en.wikipedia.org/wiki/Savitzky-Golay_filter</ulink>.</para>
+
+    <para><ulink
+    url="http://www.hpl.hp.com/techreports/2010/HPL-2010-109.pdf">On
+    the Frequency-Domain Properties of Savitzky-Golay
+    Filters</ulink>.</para>
+  </refsection>
+
+  <refsection>
+    <title>History</title>
+
+    <revhistory>
+      <revision>
+        <revnumber>0.0</revnumber>
+
+        <revdescription>Function sgolayfilt added.</revdescription>
+      </revision>
+    </revhistory>
+  </refsection>
+
+  <refsection>
+    <title>Used Functions</title>
+
+    <para><link linkend="filter">filter</link>,<link
+    linkend="sgolay">sgolay</link></para>
+  </refsection>
+</refentry>
diff --git a/scilab/modules/signal_processing/macros/sgolay.sci b/scilab/modules/signal_processing/macros/sgolay.sci
new file mode 100644 (file)
index 0000000..b53a029
--- /dev/null
@@ -0,0 +1,86 @@
+//This file is part of the Cardiovascular Wawes Analysis toolbox
+//Copyright (C) 2014 - INRIA - Serge Steer
+//This file must be used under the terms of the CeCILL.
+//This source file is licensed as described in the file COPYING, which
+//you should have received as part of this distribution.  The terms
+//are also available at
+//http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+function [B,C] = sgolay(k,nf,w)
+// Savitzky-Golay Filter Design.
+//   k : a positive scalar with integer value: the polynomial order
+//   nf : a positive scalar with integer value: the FIR filter length,
+//        must be odd and greater the k+1 
+//   w : a real vector of length nf with positive entries: the
+//       weights. If omitted no weights are applied.
+
+//   B : a real n by n array: the set of filter coefficients: the early rows
+//       of B smooth based on future values and later rows smooth based on
+//       past values, with the middle row using half future and half past.
+//       In particular, you can use row i to estimate x(k) based on the i-1
+//       preceding values and the n-i following values of x values as 
+//       y(k) = B(i,:) * x(k-i+1:k+n-i).
+
+//   C : a real n by k+1 array: the matrix of differentiation filters.  Each
+//       column of C is a differentiation filter for derivatives of order
+//       P-1 where P is the column index.  Given a signal X with length nf,
+//       an estimate of the P-th order derivative of its middle value can be
+//       found from:
+//                      (P)
+//                     X   ((nf+1)/2) = P!*C(:,P+1)'*X
+//   References:
+  
+//   - Abraham Savitzky et Marcel J. E. Golay, « Smoothing and
+//     Differentiation of Data by Simplified Least Squares Procedures »,
+//     Analytical Chemistry, vol. 8, no 36,‎ 1964, p. 1627–1639 (DOI 10.1021/ac60214a047)
+
+//   - http://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter
+
+  if type(k)<>1|int(k)<>k|size(k,'*')<>1 then
+    error(msprintf(_("%s: Wrong value for input argument #%d: An integer value expected.\n"),"sgolay",2))
+  end
+  
+  if k > nf-1 then
+    error('The degree must be less than the frame length.'),
+  end
+  if type(k)<>1|int(k)<>k|size(k,'*')<>1 then
+    error(msprintf(_("%s: Wrong type for input argument #%d: a scalar with integer value expected.\n"),...
+                   "sgolay",1))
+  end
+  if type(nf)<>1|int(nf)<>nf|size(nf,'*')<>1 then
+    error(msprintf(_("%s: Wrong type for input argument #%d: a scalar with integer value expected.\n"),...
+                   "sgolay",2))
+  end
+  
+  if modulo(nf,2)<>1 then
+    error(msprintf(_("%s: Wrong value for input argument #%d: Must be odd.\n"),"sgolay",2))
+  end
+  if nf<=k then
+    error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Argument #%d expected to be less than argument #%d.\n"),"sgolay",1,2,2,1))
+  end
+  
+  //compute the Vandermonde matrix
+  J=(ones(k+1,1)*(-(nf-1)./2:(nf-1)./2)).^((0:k)'*ones(1,nf));
+  if argn(2)==3 then
+    if type(w)<>1|~isreal(w) then
+      error(msprintf(_("%s: Wrong type for argument %d: Real vector expected.\n"),...
+                     "sgolay",3))
+    end
+    if size(w,"*")<>nf then
+      error(msprintf(_("%s: Incompatible input arguments #%d and #%d\n"),"sgolay",2,3))
+    end
+    // Check to see if all elements are positive
+    if min(w) <= 0 then
+      error(msprintf(_("%s: Wrong values for input argument #%d: Elements must be positive.\n"),...
+                     "sgolay",3))
+
+    end
+    // Diagonalize the vector to form the weighting matrix
+    J = J*diag(sqrt(w))
+  end
+  //Compute the matrix of differentiators C=J'/(J*W*J')
+  C = pinv(J); 
+  // Compute the projection matrix B
+  B = C*J;
+
+endfunction
diff --git a/scilab/modules/signal_processing/macros/sgolaydiff.sci b/scilab/modules/signal_processing/macros/sgolaydiff.sci
new file mode 100644 (file)
index 0000000..c7225a3
--- /dev/null
@@ -0,0 +1,54 @@
+//This file is part of the Cardiovascular Wawes Analysis toolbox
+//Copyright (C) 2014 - INRIA - Serge Steer
+//Copyright (C) 2020 - Stéphane Mottelet (improvement of computation at boudaries)
+//This file must be used under the terms of the CeCILL.
+//This source file is licensed as described in the file COPYING, which
+//you should have received as part of this distribution.  The terms
+//are also available at
+//http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+function D=sgolaydiff(X,order,varargin)
+  if type(X)<>1|~isreal(X) then
+    error(msprintf(_("%s: Wrong type for argument %d: Real matrix expected.\n"),"sgolaydiff",1))
+  end
+  sz=size(X)
+  if size(find(sz==1),'*')==1 then X = X(:);end
+  i1=find(sz>1,1)
+  if i1<>[] then X=matrix(X,sz(i1),-1);end
+
+  if or(size(varargin)==[2 3]) then
+    k=varargin(1);
+    nf=varargin(2);
+    if type(k)<>1|int(k)<>k|size(k,'*')<>1 then
+      error(msprintf(_("%s: Wrong value for input argument #%d: An integer value expected.\n"),"sgolay",2))
+    end
+    if type(nf)<>1|int(nf)<>nf|size(nf,'*')<>1 then
+      error(msprintf(_("%s: Wrong type for input argument #%d: a scalar with integer value expected.\n"),...
+                     "sgolay",2))
+    end
+    if modulo(nf,2)==0 then nf=nf+1;varargin(2)=nf;end
+    [B,C] = sgolay(varargin(:))
+  else
+    C=varargin(1)
+    if type(C)<>1|~isreal(C) then
+      error(msprintf(_("%s: Wrong type for argument %d: Real matrix expected.\n"),"sgolaydiff",2))
+    end
+    [nf,n]=size(C);
+    k=n-1;
+  end
+   if size(X,1) < nf then
+    error(msprintf(_("%s: Wrong size for argument %d: At least %d expected.\n"),"sgolaydiff",1,nf))
+  end
+  if order>k then
+    error(msprintf(_("%s: Wrong value for argument %d: At most %d expected.\n"),"sgolaydiff",2,n-1))
+  end
+  nf2=floor(nf/2);
+  //extend X on its boundaries by adding, in reverse order and vertically
+  //symmetrized, copies of the first (nf−1)/2 points at the beginning and 
+  //copies of the last (nf−1)/2 points at the end.
+  X = [X(1)-X(nf2+1:-1:2,:);X;2*X($)-X($-1:-1:$-nf2,:)];
+  for k=1:size(X,2)
+    D(:,k)=((-1)^order*order)*conv(X(:,k),C(:,order+1)',"valid");
+  end
+
+  D=matrix(D,sz)
+endfunction
diff --git a/scilab/modules/signal_processing/macros/sgolayfilt.sci b/scilab/modules/signal_processing/macros/sgolayfilt.sci
new file mode 100644 (file)
index 0000000..343ea80
--- /dev/null
@@ -0,0 +1,66 @@
+//This file is part of the Cardiovascular Wawes Analysis toolbox
+//Copyright (C) 2014 - INRIA - Serge Steer
+//This file must be used under the terms of the CeCILL.
+//This source file is licensed as described in the file COPYING, which
+//you should have received as part of this distribution.  The terms
+//are also available at
+//http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+function y=sgolayfilt(X,varargin)
+//filter signal using  Savitzky-Golay Filter.
+// y=sgolayfilt(X,B)
+// y=sgolayfilt(X,k,nf [,w]) 
+//
+//   X :
+//   B : a real n by n array: the set of filter coefficients: the early rows
+//       of B smooth based on future values and later rows smooth based on
+//       past values, with the middle row using half future and half past.
+//       In particular, you can use row i to estimate x(k) based on the i-1
+//       preceding values and the n-i following values of x values as 
+//       y(k) = B(i,:) * x(k-i+1:k+n-i).
+
+//   k : a positive scalar with integer value: the polynomial order, must
+//        be odd
+//   nf : a positive scalar with integer value: the FIR filter length,
+//        must be odd and greater the k+1 
+//   w : a real vector of length nf with positive entries: the weights.
+//         If omitted no weights are applied.   
+//   
+//   References:
+  
+//   - Abraham Savitzky et Marcel J. E. Golay, « Smoothing and
+//     Differentiation of Data by Simplified Least Squares Procedures »,
+//     Analytical Chemistry, vol. 8, no 36,‎ 1964, p. 1627–1639 (DOI 10.1021/ac60214a047)
+
+//   - http://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter
+  if type(X)<>1|~isreal(X) then
+    error(msprintf(_("%s: Wrong type for argument %d: Real matrix expected.\n"),...
+                                      "sgolayfilt",1))
+  end
+  if size(X,1) == 1 then X = X(:);end
+
+  if or(size(varargin)==[2 3]) then
+    B = sgolay(varargin(:))
+  else
+    B=varargin(1)
+    if type(B)<>1|~isreal(B) then
+      error(msprintf(_("%s: Wrong type for argument %d: Real matrix expected.\n"),...
+                                      "sgolayfilt",2))
+    end
+  end
+  n=size(B,1);
+  if size(X,1) < n, 
+    error(msprintf(_("%s: Wrong size for argument %d: At least %d expected.\n"),"sgolayfilt",1,size(B,1)))
+  end
+
+  // The first k rows of B are used to filter the first k points
+  // of the data set based on the first n points of the data set.
+  // The last k rows of B are used to filter the last k points
+  // of the data set based on the last n points of the dataset.
+  // The remaining data is filtered using the central row of B.
+  k = floor(n/2);
+  z=[]
+  for kc=1:size(X,2)
+    z(:,kc) = filter(B(k+1,:), 1, X(:,kc));
+  end
+  y = [B(1:k,:)*X(1:n,:) ; z(n:$,:) ; B(k+2:n,:)*X($-n+1:$,:) ];
+endfunction