- 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:
-----------
<member>
<link linkend="numderivative">numderivative</link>
</member>
+ <member>
+ <link linkend="sgolaydiff">sgolaydiff</link>
+ </member>
</simplelist>
</refsection>
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
<member>
<link linkend="ltitr">ltitr</link>
</member>
+ <member>
+ <link linkend="sgolayfilt">sgolayfilt</link>
+ </member>
</simplelist>
</refsection>
</refentry>
--- /dev/null
+<?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>
--- /dev/null
+<?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>
--- /dev/null
+<?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>
--- /dev/null
+//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
--- /dev/null
+//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
--- /dev/null
+//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
projects / scilab.git / commitdiff
