Return to previous page

Contents of file 'sgfilter/sgfilter.w':




Deprecated: Function split() is deprecated in /storage/content/45/2011745/jonsson.eu/public_html/php/htmlicise.php on line 46
    1   % File:        sgfilter.w [CWEB source code]
    2   % Residence:   http://jonsson.eu/programs/cweb/sgfilter/
    3   % Created:     January 23, 2006 [v.1.0]
    4   % Last change: December 4, 2011 [v.1.6]
    5   % Author:      Fredrik Jonsson <http://jonsson.eu>
    6   % Description: The CWEB source code for the SGFILTER stand-alone implementation
    7   %              of the Savitzky-Golay smoothing filter. For information on the
    8   %              CWEB programming language, see
    9   %              http://www.literateprogramming.com.
   10   % Compilation: Compile this program by using the enclosed Makefile, or use
   11   %              the blocks of the Makefile as listed in the documentation
   12   %              (file sgfilter.ps or sgfilter.pdf). The C source code (as
   13   %              generated from this CWEB code) conforms to the ANSI standard
   14   %              for the C programming language (which is equivalent to the
   15   %              ISO C90 standard for C).
   16   %
   17   % Copyright (C) 2011, Fredrik Jonsson <http://jonsson.eu>
   18   % Non-commercial copying welcome.
   19   %
   20   \input ext/epsf.tex   % For the inclusion of Encapsulated PostScript graphics
   21   \input ext/eplain.tex % To access equation numbering and other cross-referencing
   22   \def\version{1.6}
   23   \def\lastrevdate{December 4, 2011}
   24   \font\eightcmr=cmr8
   25   \font\tensc=cmcsc10
   26   \font\eightcmssq=cmssq8
   27   \font\eightcmssqi=cmssqi8
   28   \font\twentycmcsc=cmcsc10 at 20 truept
   29   \def\sgfilter{{\eightcmr SGFILTER\spacefactor1000}}
   30   \def\ANSI{{\eightcmr ANSI\spacefactor1000}} % The language standard we stick to
   31   \def\GNU{{\eightcmr GNU\spacefactor1000}}       % GNU is Not Unix
   32   \def\GCC{{\eightcmr GCC\spacefactor1000}}       % The GNU C-compiler
   33   \def\CEE{{\eightcmr C\spacefactor1000}}         % The C programming language
   34   \def\FORTRAN{{\eightcmr Fortran\spacefactor1000}} % The Fortran language
   35   \def\CPP{\CEE{\tt++}}                           % The C++ programming language
   36   \def\ANSICEE{{\eightcmr ANSI~C\spacefactor1000}}% The ANSI-C standard language
   37   \def\CWEB{{\eightcmr CWEB\spacefactor1000}}     % The CWEB programming language
   38   \def\UNIX{{\eightcmr UNIX\spacefactor1000}}
   39   \def\CYGWIN{{\eightcmr CYGWIN\spacefactor1000}}
   40   \def\GNUPLOT{{\eightcmr GNUPLOT\spacefactor1000}}
   41   \def\CTANGLE{{\eightcmr CTANGLE\spacefactor1000}}
   42   \def\CWEAVE{{\eightcmr CWEAVE\spacefactor1000}}
   43   \def\mod{\mathop{\rm mod}\nolimits}  % The real part of a complex number
   44   \def\re{\mathop{\rm Re}\nolimits}  % The real part of a complex number
   45   \def\im{\mathop{\rm Im}\nolimits}  % The imaginary part of a complex number
   46   \def\backslash{\char'134}          % The `\' character
   47   \def\vertbar{\char'174}            % The `|' character
   48   \def\dollar{\char'044}             % The `$' character
   49   \def\tilde{\char'176}              % The `~' character
   50   \def\tothepower{\char'136}         % The `^' character
   51   \def\onehalf{{\textstyle{{1}\over{2}}}}
   52   \def\threefourth{{\textstyle{{3}\over{4}}}}
   53   \def\endalg{\vrule height 1.4ex width .6ex depth .4ex} % Rectangular bullet
   54   \def\ie{i.\thinspace{e.}~\ignorespaces}
   55   \def\eg{e.\thinspace{g.}~\ignorespaces}
   56   \let\,\thinspace
   57   %
   58   % Define a handy macro for listing the steps of an algorithm.
   59   %
   60   \newdimen\aitemindent \aitemindent=26pt
   61   \newdimen\aitemleftskip \aitemleftskip=36pt
   62   \def\aitem[#1]{\smallbreak\noindent\hbox to 10pt{}%
   63     \hbox to\aitemindent{\bf #1\hfill}%
   64     \hangindent\aitemleftskip\ignorespaces}
   65   %
   66   % Define a handy macro for bibliographic references.
   67   %
   68   \newdimen\refitemindent \refitemindent=18pt
   69   \def\refitem[#1]{\smallbreak\noindent%
   70     \hbox to\refitemindent{[#1]\hfill}%
   71     \hangindent\refitemindent\ignorespaces}
   72   \def\references{\medskip\goodbreak\noindent{\bf References}\nobreak\smallskip}
   73   %
   74   % Define a handy macro for nicely typeset variable descriptions.
   75   %
   76   \newdimen\varitemindent \varitemindent=100pt
   77   \def\varitem[#1]#2{\smallbreak\noindent\hbox to 20pt{}%
   78     \hbox to\varitemindent{#1\hfill}%
   79     \hangindent 120pt\ignorespaces#2\par}
   80   %
   81   % Define a handy macro for nicely typeset descriptions of command line options.
   82   %
   83   \newdimen\optitemindent \optitemindent=80pt
   84   \def\optitem[#1]#2{\smallbreak\noindent\hbox to\hsize{\hskip40pt{#1}\hfill}
   85     \hangindent\optitemindent\ignorespaces{\hskip\optitemindent{#2}}\smallskip}
   86   %
   87   % Define a handy macro for the list of program revisions.
   88   %
   89   \newdimen\citemindent \citemindent=80pt
   90   \newdimen\citemleftskip \citemleftskip=90pt
   91   \def\citem[#1]{\smallbreak\noindent\hbox to 10pt{}%
   92     \hbox to\citemindent{\bf #1\hfill}%
   93     \hangindent\citemleftskip\ignorespaces}
   94   \def\citindent{\smallbreak\noindent\hbox to 10pt{}%
   95     \hbox to\citemindent{\hfil}%
   96     \hangindent\citemleftskip\ignorespaces}
   97   %
   98   % Define a handy macro for the typesetting of quotes of text.
   99   %
  100   \newdimen\quoteindent \quoteindent=40pt
  101   \def\quote#1{\par{\advance\leftskip by\quoteindent
  102     \advance\rightskip by\quoteindent
  103     \medskip\noindent{``#1''}\medskip}\par}
  104   %
  105   % Define a handy macro for the typesetting of figure captions. The double
  106   % '{{' and '}}' are here in place to prevent the increased \leftskip and
  107   % \rightskip to apply globally after the macro has been expanded.
  108   %
  109   \newdimen\figcapindent \figcapindent=\quoteindent
  110   \def\figcaption[#1]#2{{\advance\leftskip by\figcapindent
  111     \advance\rightskip by\figcapindent
  112     \smallskip\noindent{\bf Figure #1.} {#2}\smallskip}}
  113   %
  114   % Define the \beginvrulealign and \endvrulealign macros as described in
  115   % Donald Knuth's The TeXbook, Appendix D: Dirty Tricks. These macros are
  116   % used in typesetting nicely looking tables.
  117   %
  118   \def\beginvrulealign{\setbox0=\vbox\bgroup}
  119   \def\endvrulealign{\egroup % now \box0 holds the entire alignment
  120     \setbox0=\vbox{\setbox2=\hbox{\vrule height\ht0 depth\dp0 width0pt}
  121       \unvbox0 \setbox0=\lastbox % now \box0 is the bottom row
  122       \nointerlineskip \copy0 % put it back
  123       \global\setbox1=\hbox{} % initialize box that will contain rules
  124       \setbox4=\hbox{\unhbox0 % now open up the bottom row
  125         \loop \skip0=\lastskip \unskip % remove tabskip glue
  126         \advance\skip0 by-.4pt % rules are .4 pt wide
  127         \divide\skip0 by 2
  128         \global\setbox1=\hbox{\hskip\skip0\vrule\hskip\skip0
  129           \unhbox2\unhbox1}%
  130         \setbox2=\lastbox % remove alignment entry
  131         \ifhbox2 \setbox2=\hbox{\kern\wd2}\repeat}}%
  132     \hbox{\rlap{\box0}\box1}} % superimpose the alignment on the rules
  133   %
  134   % Make sure that the METAPOST logo can be loaded even in plain TeX.
  135   %
  136   \ifx\MP\undefined
  137       \ifx\manfnt\undefined
  138               \font\manfnt=logo10
  139       \fi
  140       \ifx\manfntsl\undefined
  141               \font\manfntsl=logosl10
  142       \fi
  143       \def\MP{{\ifdim\fontdimen1\font>0pt \let\manfnt = \manfntsl \fi
  144         {\manfnt META}\-{\manfnt POST}}\spacefactor1000 }%
  145   \fi
  146   \ifx\METAPOST\undefined \let\METAPOST=\MP \fi
  147   \datethis
  148   
  149   @*Introduction.
  150   \vskip 120pt
  151   \centerline{\twentycmcsc SGFilter}
  152   \vskip 20pt
  153   \centerline{A stand-alone implementation of the Savitzky--Golay smoothing
  154     filter.}
  155   \vskip 2pt
  156   \centerline{(Version \version\ of \lastrevdate)}
  157   \vskip 10pt
  158   \centerline{Written by Fredrik Jonsson}
  159   \vskip 10pt
  160   \centerline{\epsfbox{fig0.1}}
  161   \vfill\eject
  162   \noindent
  163   This document was automatically extracted from the \CWEB\ master source code
  164   for the \sgfilter\ program and typeset in the Computer Modern typeface using
  165   plain \TeX\ and \MP.
  166   The source code and documentation of this
  167   program is electronically available at
  168   \.{http://jonsson.eu/programs/cweb/sgfilter/}.
  169   \vfill
  170   \line{Copyright \copyright\ Fredrik Jonsson \the\year\hfil}
  171   \medskip
  172   \noindent{All rights reserved. No part of this publication may be reproduced,
  173   stored in a retrieval system, or transmitted, in any form, or by any means,
  174   electronic, mechanical, photo-copying, recording, or otherwise, without the
  175   prior consent of the author. Non-commercial copying welcome.\hfil}
  176   \medskip
  177   \line{Printed on \today, at \hours\hfil}
  178   \medskip
  179   \line{\TeX\ is a trademark of the American Mathematical Society\hfil}
  180   
  181   @*The Savitzky--Golay smoothing filter.
  182   The Savitzky--Golay smoothing filter was originally presented in 1964 by
  183   Abraham Savitzky\footnote{$\dagger$}{Abraham Savitzky (1919--1999) was
  184   an American analytical chemist. (Wikipedia)} and Marcel J.~E.
  185   Golay\footnote{$\ddagger$}{Marcel J.~E.~Golay (1902--1989) was
  186   a Swiss-born mathematician, physicist, and information theorist, who applied
  187   mathematics to real-world military and industrial problems. (Wikipedia)}
  188   in their paper ``Smoothing and Differentiation of Data by Simplified Least
  189   Squares Procedures'', Anal.~Chem., {\bf 36}, 1627--1639 (1964).
  190   % \.{http://dx.doi.org/10.1021/ac60214a047}
  191   Being chemists and physicists, at the time of publishing associated with
  192   the Perkin--Elmer Corporation (still today a reputable manufacturer of
  193   equipment for spectroscopy), they found themselves often encountering noisy
  194   spectra where simple noise-reduction techniques, such as running averages,
  195   simply were not good enough for extracting well-determined characteristica
  196   of spectral peaks.
  197   In particular, any running averaging tend to flatten and widening peaks in
  198   a spectrum, and as the peak width is an important parameter when determining
  199   relaxation times in molecular systems, such noise-reduction techniques are
  200   clearly non-desirable.
  201   
  202   The main idea presented by Savitzky and Golay was a work-around avoiding the
  203   problems encountered with running averages, while still maintaining the
  204   smoothing of data and preserving features of the distribution such as relative
  205   maxima, minima and width.
  206   To quote the original paper on the target and purpose:
  207   \quote{This paper is
  208   concerned with computational methods for the removal of the random noise from
  209   such information, and with the simple evaluation of the first few derivatives
  210   of the information with respect to the graph abscissa.
  211   [$\ldots$]
  212   The objective here is to present specific methods for handling current
  213   problems in the processing of such tables of analytical data. The methods
  214   apply as well to the desk calculator, or to simple paper and pencil
  215   operations for small amounts of data, as they do to the digital computer for
  216   large amounts of data, since their major utility is to simplify and speed up
  217   the processing of data.}
  218   
  219   @ The work-around presented by Savitzky and Golay for avoiding distortion of
  220   peaks or features in their spectral data is essentially based on the
  221   idea to perform a linear regression of some polynomial {\sl individually for
  222   each sample}, followed by the {\sl evaluation of that polynomial exactly at the
  223   very position of the sample}. While this may seem a plausible idea, the actual
  224   task of performing a separate regression for each point easily becomes a very
  225   time-consuming task unless we make a few observations about this problem.
  226   However (and this is the key point in the method), for the regression of a
  227   polynomial of a finite power, say of an order below 10, the coefficients
  228   involved in the actual regression may be computed once and for all in an
  229   early stage, followed by performing a {\sl convolution} of the discretely
  230   sampled input data with the coefficient vector. As the coefficient vector
  231   is significantly shorter than the data vector, this convolution is fast and
  232   straightforward to implement.
  233   
  234   The starting point in deriving the algorithm for the Savitzky--Golay smoothing
  235   filter is to consider a smoothing method in which an equidistantly spaced
  236   set of $M$ samples $f_n$, $n=1,\ldots,M$ are linearly combined to form a
  237   filtered value $h_k$ according to\footnote{$\dagger$}{More generally,
  238   Eq.~\eqref{eq:10} can be interpreted as a discretized version of the
  239   convolution between a kernel $c(x)$ and a function $f(x)$,
  240   $$
  241     h(x)=\int^{\Delta_R}_{-\Delta_L}c(s) f(x+s)\,ds
  242         \rightarrow\sum^{n_R}_{j=-n_L} c_j f_{k+j}.
  243   $$
  244   Notice, however, the different sign of $s$ as compared to the standard form
  245   of the convolution integral in mathematics, where the argument of the
  246   function usually yields ``$f(x-s)$''.}
  247   $$
  248     h_k=\sum^{n_R}_{j=-n_L} c_j f_{k+j},
  249     \eqdef{eq:10}
  250   $$
  251   where $n_L$ is the number of samples ``to the left'' and $n_R$ the number of
  252   samples ``to the right'' of the centre index $k$.
  253   Notice that the running average smoothing corresponds to the case where
  254   all coefficients $c_n$ are equal, with $c_n=1/(n_L+n_R+1)$.
  255   The idea of the Savitzky--Golay filter is, however, to find the set of $c_n$
  256   which better preserves the shape of features present in the sampled profile.
  257   The approach is to make a linear regression of a polynomial to the $n_L+n_R+1$
  258   samples in the window around sample $k$, and then evaluating this polynomial
  259   at that very sample, for all $k$ from $1$ to $M$.
  260   
  261   We consider the regression of an $m$:th degree polynomial
  262   $$
  263     p(x)=a_0+a_1(x-x_k)+a_2(x-x_k)^2+\ldots+a_m (x-x_k)^m
  264     \eqdef{eq:20}
  265   $$
  266   to the set of $n_L$ samples to the left and $n_R$ to the right of sample $f_k$,
  267   including the sample inbetween. Evaluating this polynomial at $x=x_k$ is
  268   particularly easy, as we at that point simply have $p(x_k)=a_0$.
  269   The linear regression of this polynomial to the samples $f_{k+j}$, with
  270   $-n_L\le j\le n_R$, means that we look for the least-square approximated
  271   solution to the linear system of $n_L+n_R+1$ equations
  272   $$
  273     \eqalign{
  274       p(x_{k-n_L})&=a_0+a_1(x_{k-n_L}-x_k)+a_2(x_{k-n_L}-x_k)^2
  275                       +\ldots+a_m (x_{k-n_L}-x_k)^m\approx f_{k-n_L},\cr
  276                 &\hskip60pt\vdots\hskip100pt\vdots\cr
  277       p(x_{k-1})&=a_0+a_1(x_{k-1}-x_k)+a_2(x_{k-1}-x_k)^2
  278                       +\ldots+a_m (x_{k-1}-x_k)^m\approx f_{k-1},\cr
  279       p(x_{k})&=a_0\approx f_{k}\cr
  280       p(x_{k+1})&=a_0+a_1(x_{k+1}-x_k)+a_2(x_{k+1}-x_k)^2
  281                       +\ldots+a_m (x_{k+1}-x_k)^m\approx f_{k+1},\cr
  282                 &\hskip60pt\vdots\hskip100pt\vdots\cr
  283       p(x_{k+n_R})&=a_0+a_1(x_{k+n_R}-x_k)+a_2(x_{k+n_R}-x_k)^2
  284                       +\ldots+a_m (x_{k+n_R}-x_k)^m\approx f_{k+n_R},\cr
  285     }
  286     \eqdef{eq:30}
  287   $$
  288   which (under the assumption that $n_L+n_R+1>m+1$) provides an overdetermined
  289   system for the $m+1$ coefficients $a_j$ which can be expressed in matrix
  290   form as
  291   $$
  292     {\bf A}\cdot{\bf a}\equiv
  293     \underbrace{
  294       \pmatrix{
  295         1&(x_{k-n_L}-x_k)&(x_{k-n_L}-x_k)^2&\cdots&(x_{k-n_L}-x_k)^m\cr
  296         \vdots& &\vdots& &\vdots\cr
  297         1&(x_{k-1}-x_k)&(x_{k-1}-x_k)^2&\cdots&(x_{k-1}-x_k)^m\cr
  298         1&0&0&\cdots&0\cr
  299         1&(x_{k+1}-x_k)&(x_{k+1}-x_k)^2&\cdots&(x_{k+1}-x_k)^m\cr
  300         \vdots& &\vdots& &\vdots\cr
  301         1&(x_{k+n_R}-x_k)&(x_{k+n_R}-x_k)^2&\cdots&(x_{k+n_R}-x_k)^m\cr
  302       }
  303     }_{[(n_L+n_R+1)\times(m+1)]}
  304     \underbrace{
  305       \pmatrix{
  306         a_0\cr
  307         a_1\cr
  308         a_2\cr
  309         \vdots\cr
  310         a_m\cr
  311       }
  312       \vphantom{
  313           \pmatrix{
  314             ()^m\cr
  315             \vdots\cr
  316             ()^m\cr
  317             1\cr
  318             ()^m\cr
  319             \vdots\cr
  320             ()^m\cr
  321           }
  322       }
  323     }_{[(m+1)\times1]}
  324     =
  325     \underbrace{
  326       \pmatrix{
  327         f_{k-n_L}\cr
  328         \vdots\cr
  329         f_{k-1}\cr
  330         f_{k}\cr
  331         f_{k+1}\cr
  332         \vdots\cr
  333         f_{k+n_R}\cr
  334       }
  335       \vphantom{
  336           \pmatrix{
  337             ()^m\cr
  338             \vdots\cr
  339             ()^m\cr
  340             1\cr
  341             ()^m\cr
  342             \vdots\cr
  343             ()^m\cr
  344           }
  345       }
  346     }_{[(m+1)\times1]}
  347     \equiv{\bf f}.
  348     \eqdef{eq:40}
  349   $$
  350   The least squares solution to Eq.~\eqref{eq:40} is obtained by multiplying its
  351   left- and right-hand sides by the transpose of the system matrix ${\bf A}$,
  352   followed by solving the resulting $[(m+1)\times(m+1)]$-system of linear
  353   equations for ${\bf a}$,
  354   $$
  355     {\bf A}^{\rm T}\cdot({\bf A}\cdot{\bf a})=
  356       ({\bf A}^{\rm T}\cdot{\bf A})\cdot{\bf a}=
  357       {\bf A}^{\rm T}\cdot{\bf f}
  358     \qquad\Leftrightarrow\qquad
  359     {\bf a}=
  360       ({\bf A}^{\rm T}\cdot{\bf A})^{-1}\cdot
  361       ({\bf A}^{\rm T}\cdot{\bf f})
  362     \eqdef{eq:50}
  363   $$
  364   Recapitulate that what we here target is the evaluation of $p(x_k)=a_0$, which
  365   according to Eqs.~\eqref{eq:20} and~\eqref{eq:50} is equivalent to evaluating
  366   the first (``zero:th'') row of the solution for ${\bf a}$, or
  367   $$
  368     p(x_k)=a_0=
  369       \big[
  370         \underbrace{
  371           ({\bf A}^{\rm T}\cdot{\bf A})^{-1}
  372         }_{[(m+1)\times(m+1)]}
  373         \cdot
  374         \underbrace{
  375           ({\bf A}^{\rm T}\cdot{\bf f})
  376         }_{[(m+1)\times1]}
  377       \big]_{\{{\rm row\ }0\}}
  378     \eqdef{eq:60}
  379   $$
  380   So, having arrived at this result for the regression, we clearly have a
  381   solution for $a_0$ depending on the actual function values $f_{k+j}$ in the
  382   vicinity of sample $f_k$. Doesn't this mean that we nevertheless need to
  383   repeat the regression for every single sample to be included in the smoothing?
  384   Moreover, how does the result presented in Eq.~\eqref{eq:60} relate to the
  385   original mission, which we recapitulate was to find a general way of computing
  386   the coefficients $c_j$ in the kernel of the convolution in Eq.~\eqref{eq:10}?
  387   
  388   The answer to these questions lies in expanding Eq.~\eqref{eq:60} to yield the
  389   expression for the first (``zero:th'') row as
  390   $$
  391     \eqalign{
  392     p(x_k)=a_0
  393          &=\big[
  394           ({\bf A}^{\rm T}\cdot{\bf A})^{-1}\cdot
  395           ({\bf A}^{\rm T}\cdot{\bf f})\big]_{\{{\rm row\ }0\}}\cr
  396          &=\sum^{m+1}_{j=1}
  397           \big[({\bf A}^{\rm T}\cdot{\bf A})^{-1}\big]_{0j}
  398           \big[{\bf A}^{\rm T}\cdot{\bf f}\big]_{j}\cr
  399          &=\sum^{m+1}_{j=1}
  400           \big[({\bf A}^{\rm T}\cdot{\bf A})^{-1}\big]_{0j}
  401           \sum^{m+1}_{k=1}\big[{\bf A}^{\rm T}\big]_{jk}\big[{\bf f}\big]_{k}\cr
  402     }
  403     \eqdef{eq:70}
  404   $$
  405   and observing that the $n$:th coefficient $c_n$ is obtained as equal to the
  406   coefficient $a_0$ whenever ${\bf f}$ is replaced by the unit vector
  407   ${\bf e}_n$ (with all elements zero except for the unitary $n$:th element).
  408   Hence
  409   $$
  410     \eqalign{
  411       c_n&=\sum^{m+1}_{j=1}
  412         \big[({\bf A}^{\rm T}\cdot{\bf A})^{-1}\big]_{0j}
  413         \sum^{m+1}_{k=1}\big[{\bf A}^{\rm T}\big]_{jk}\big[{\bf e}_n\big]_{k}\cr
  414       &=\sum^{m+1}_{j=1}
  415         \big[({\bf A}^{\rm T}\cdot{\bf A})^{-1}\big]_{0j}
  416         \sum^{m+1}_{k=1}\big[{\bf A}^{\rm T}\big]_{jk}\delta_{nk}\cr
  417       &=\sum^{m+1}_{j=1}
  418         \big[({\bf A}^{\rm T}\cdot{\bf A})^{-1}\big]_{0j}
  419         \big[{\bf A}^{\rm T}\big]_{jn}\cr
  420       &=\sum^{m+1}_{j=1}
  421         \big[({\bf A}^{\rm T}\cdot{\bf A})^{-1}\big]_{0j}A_{nj}\cr
  422     }
  423     \eqdef{eq:80}
  424   $$
  425   This concludes the derivation of the coefficients in the Savitzky--Golay
  426   filter, which may be computed once and for all in the beginning and afterwards
  427   used as the (for any practical purposes short) kernel $c_j$ in the convolution
  428   described by Eq.~\eqref{eq:10}.
  429   
  430   @ In 2000, editors James Riordon, Elizabeth Zubritsky, and Alan Newman of
  431   {\sl Analytical Chemistry} made a review article of what they had identified
  432   as the top-ten seminal papers in the history of the journal, based on the
  433   number of citations.
  434   Among the listed papers, of which some were written by Nobel-laureates-to-be,
  435   the original paper by Savitzky and Golay makes a somewhat odd appearance, as
  436   it not only concerns mainly numerical analysis, but also because it actually
  437   includes \FORTRAN\ code for the implementation.
  438   The review article\footnote{$\dagger$}{Available at
  439   \.{http://pubs.acs.org/doi/pdf/10.1021/ac002801q}} concludes the discussion
  440   of the Savitzky--Golay smoothing filter with a reminiscence by Abraham Savitzky
  441   about the work:
  442   \quote{In thinking about why the technique has been so widely used, I've come
  443   to the following conclusions. First, it solves a common problem--the reduction
  444   of random noise by the well-recognized least-squares technique. Second, the
  445   method was spelled out in detail in the paper, including tables and a sample
  446   computer subroutine. Third, the mathematical basis for the technique, although
  447   explicitly and rigorously stated in the article, was separated from a
  448   completely nonmathematical explanation and justification. Finally, the
  449   technique itself is simple and easy to use, and it works.}
  450   
  451   @ As a final remark on the Savitzky--Golay filtering algorithm, a few points
  452   on the actual implementation of the convolution need to be made.
  453   While the \sgfilter\ program relies on the method for computation of the
  454   Savitzky--Golay coefficients as presented in {\sl Numerical Recipes in C},
  455   2nd Edn~(Cambridge University Press, New York, 1994), it must be emphasized
  456   that the suggestion there presented for the convolution, which is to apply
  457   the |convlv| routine of {\sl Numerical Recipes in C}, is significantly
  458   increasing the complexity and memory consumption in the filtering.
  459   In particular, the |convlv| routine in turn relies on consistent calls to the
  460   |twofft| routine, which in order to deliver proper data needs to be supplied
  461   with a return vector of {\sl twice} the size of the input vector.
  462   In addition, |convlv| requires the size of the input vector to be of an
  463   integer power of two (say, $M=1024$, $4096$, etc.), which may be acceptable
  464   for one-off tests but is a rather inconvenient limitation for any more general
  465   applications.
  466   
  467   Whether the \sgfilter\ program should employ the convolution engine supplied by
  468   the |convlv| routine (recommended in {\sl Numerical Recipes in C}) or the
  469   direct convolution as implemented in the |sgfilter| routine (recommended by
  470   me) is controlled by the |CONVOLVE_WITH_NR_CONVLV| definition in the
  471   \.{sgfilter.h} header file. With reference to the above issues with |convlv|,
  472   I strongly advise keeping the default (\.{0}) setting for
  473   |CONVOLVE_WITH_NR_CONVLV|.
  474   
  475   @*Revision history of the program.
  476   \medskip
  477   \citem[2006-01-18]{[v.1.0]} {\tt <fj@@phys.soton.ac.uk>}\hfill\break
  478   First properly working version of the \sgfilter\ program.
  479   
  480   \citem[2006-01-20]{[v.1.1]} {\tt <fj@@phys.soton.ac.uk>}\hfill\break
  481   Added the test case for Savitzky--Golay filtering, modeling an underlying
  482   function $g(x)$ with superimposed Gaussian noise as
  483   $$
  484     f(x)=\underbrace{
  485            \cos(3x)\sin^2(x^3)+4\sum^4_{k=1}\exp(-(x-x_k)^2/w^2_k)
  486          }_{g(x)} + \underbrace{Vu(x)\vphantom{\sum^4_{k=1}}}_{\rm noise}
  487   $$
  488   where $u(x)$ is a normally distributed stochastic variable of mean zero and
  489   unit variance, $V$ is the local variance as specified arbitrarily, and the
  490   remaining parameters $(x_k,w_k)$ are the positions and widths of four
  491   Gaussian peaks superimposed onto the otherwise trigonometric expression
  492   for the underlying function.
  493   
  494   \citem[2006-01-21]{[v.1.2]} {\tt <fj@@phys.soton.ac.uk>}\hfill\break
  495   Changed the streaming of output (filtered) data so that the |stdout| stream
  496   now is directed to file whenever \.{-o} or \.{\--outputfile} options are
  497   present at the calling command line.
  498   
  499   \citem[2006-05-06]{[v.1.3]} {\tt <fj@@phys.soton.ac.uk>}\hfill\break
  500   Added an introductory section documenting the derivation of the
  501   Savitzky--Golay filter as such. Always nice to have at hand when it comes
  502   to actually understanding why certain parameters in the filtering need to
  503   be in certain ranges. Also added automatic support for extracting the number
  504   of input samples automatically from the input file, hence making the \.{-M}
  505   and \.{--num\_samples} option obsolete.
  506   
  507   \citem[2006-05-06]{[v.1.4]} {\tt <fj@@phys.soton.ac.uk>}\hfill\break
  508   Replaced the convolution from the previously used |convlv| routine to a
  509   brute-force but more economical one, which now does not rely on the input
  510   data being a set of $2^N$ samples for some $N$. However, I have chosen to
  511   keep the old implementation, which can be re-applied simply by changing
  512   the definition of |CONVOLVE_WITH_NR_CONVLV| to ``(\.{1}).''
  513   
  514   \citem[2009-11-01]{[v.1.5]} {\tt <http://jonsson.eu>}\hfill\break
  515   Included the \.{example.c} file in the \CWEB\ source of \sgfilter, for
  516   automatically is generated from the master \CWEB\ source when passed through
  517   \CTANGLE. A very nifty way indeed for keeping updated test cases.
  518   
  519   
  520   \citem[2011-12-04]{[v.1.6]} {\tt <http://jonsson.eu>}\hfill\break
  521   In the block concerned with the redirection of |stdout| when delivering
  522   filtered data, I found strange behaviour when executing \sgfilter\ under
  523   Windows. What happens is that any redirection of |stdout| back to terminal
  524   output in Windows naturally must be done in a different way than with the
  525   |freopen("/dev/tty", "a", stdout)| which is accepted by OS X (Free BSD),
  526   Linux, or any other \UNIX-like platforms. Hence, I simply added a (primitive)
  527   check on the platform type in the header file.
  528   
  529   @*Compiling the source code. The program is written in \CWEB, generating
  530   \ANSICEE\ (ISO C99) conforming source code and documentation as plain
  531   \TeX-source, and is to be compiled using the sequences as outlined in the
  532   \.{Makefile} listed below.
  533   For general information on literate programming, \CTANGLE, or \CWEAVE, see
  534   \.{http://www.literateprogramming.com}.
  535   \bigskip
  536   {\obeylines\obeyspaces\tt
  537   \#
  538   \# Makefile designed for use with ctangle, cweave, gcc, and plain TeX.
  539   \#
  540   \# Copyright (C) 2002-2011, Fredrik Jonsson <http://jonsson.eu>
  541   \#
  542   \# The CTANGLE program converts a CWEB source document into a C program
  543   \# which may be compiled in the usual way. The output file includes \#line
  544   \# specifications so that debugging can be done in terms of the CWEB source
  545   \# file.
  546   \#
  547   \# The CWEAVE program converts the same CWEB file into a TeX file that may
  548   \# be formatted and printed in the usual way. It takes appropriate care of
  549   \# typographic details like page layout and the use of indentation, italics,
  550   \# boldface, etc., and it supplies extensive cross-index information that it
  551   \# gathers automatically.
  552   \#
  553   \# CWEB allows you to prepare a single document containing all the informa-
  554   \# tion that is needed both to produce a compilable C program and to produce
  555   \# a well-formatted document describing the program in as much detail as the
  556   \# writer may desire.  The user of CWEB ought to be familiar with TeX as well
  557   \# as C.
  558   \#
  559   PROJECT  = sgfilter
  560   CTANGLE  = ctangle
  561   CWEAVE   = cweave
  562   CC       = gcc
  563   CCOPTS   = -O2 -Wall -ansi -std=iso9899:1999 -pedantic
  564   LNOPTS   = -lm
  565   TEX      = tex
  566   DVIPS    = dvips
  567   DVIPSOPT = -ta4 -D1200
  568   PS2PDF   = ps2pdf
  569   METAPOST = mpost
  570   ~   ~
  571   all: \$(PROJECT) \$(PROJECT).pdf
  572   ~   ~
  573   \$(PROJECT): \$(PROJECT).o
  574   ~        \$(CC) \$(CCOPTS) -o \$(PROJECT) \$(PROJECT).o \$(LNOPTS)
  575   ~   ~
  576   \$(PROJECT).o: \$(PROJECT).c
  577   ~        \$(CC) \$(CCOPTS) -c \$(PROJECT).c
  578   ~   ~
  579   \$(PROJECT).c: \$(PROJECT).w
  580   ~        \$(CTANGLE) \$(PROJECT) \$(PROJECT).c
  581   ~   ~
  582   \$(PROJECT).pdf: \$(PROJECT).ps
  583   ~        \$(PS2PDF) \$(PROJECT).ps \$(PROJECT).pdf
  584   ~   ~
  585   \$(PROJECT).ps: \$(PROJECT).dvi
  586   ~        \$(DVIPS) \$(DVIPSOPT) \$(PROJECT).dvi -o \$(PROJECT).ps
  587   ~   ~
  588   \$(PROJECT).dvi: \$(PROJECT).tex
  589   ~        \$(TEX) \$(PROJECT).tex
  590   ~   ~
  591   \$(PROJECT).tex: \$(PROJECT).w
  592   ~        \$(CWEAVE) \$(PROJECT)
  593   ~   ~
  594   clean:
  595   ~        -rm -Rf \$(PROJECT) *~ *.c *.h *.o *.exe *.dat *.pdf *.mp *.trj *.mpx
  596   ~        -rm -Rf *.tex *.aux *.log *.toc *.idx *.scn *.dvi *.ps *.1 *.eps
  597   ~   ~
  598   archive:
  599   ~        make -ik clean
  600   ~        tar --gzip --directory=../ -cf ../\$(PROJECT).tar.gz \$(PROJECT)
  601   }
  602   \bigskip
  603   This \.{Makefile} essentially executes two major calls. First, the \CTANGLE\
  604   program parses the \CWEB\ source file \.{sgfilter.w} to extract a \CEE\ source
  605   file \.{sgfilter.c} which may be compiled into an executable program using any
  606   \ANSICEE\ conformant compiler. The output source file \.{sgfilter.c} includes
  607   \.{\#line} specifications so that any debugging conveniently can be done in
  608   terms of line numbers in the original \CWEB\ source file \.{sgfilter.w}.
  609   Second, the \CWEAVE\ program parses the same \CWEB\ source file \.{sgfilter.w}
  610   to extract a plain \TeX\ file \.{sgfilter.tex} which may be compiled into a
  611   PostScript or PDF document.
  612   The document file \.{sgfilter.tex} takes appropriate care of typographic
  613   details like page layout and text formatting, and supplies extensive
  614   cross-indexing information which is gathered automatically.
  615   In addition to extracting the documentary text, \CWEAVE\ also includes the
  616   source code in cross-referenced blocks corresponding to the descriptors as
  617   entered in the \CWEB\ source code.
  618   
  619   Having executed \.{make} (or \.{gmake} for the \GNU\ enthusiast) in the same
  620   directory where the files \.{sgfilter.w} and \.{Makefile} are located, one is
  621   left with the executable file \.{sgfilter}, being the ready-to-use compiled
  622   program, and the PostScript file \.{sgfilter.ps} (or PDF file \.{sgfilter.pdf})
  623   which contains the full documentation of the program, that is to say the
  624   document you currently are reading.
  625   Notice that on platforms running any operating system by Microsoft, the
  626   executable file will instead automatically be named \.{sgfilter.exe}.
  627   This convention also applies to programs compiled under the \UNIX-like
  628   environment \CYGWIN.
  629   
  630   @*Running the program. The program is entirely controlled by the command
  631   line options supplied when invoking the program. The syntax for executing the
  632   program is\par
  633   \medskip
  634   \line{\hskip 40pt\.{sgfilter [options]}\hfill}\par
  635   \medskip
  636   \noindent
  637   where \.{options} include the following, given in their long as well as
  638   their short forms (with prefixes `\.{--}' and `\.{-}', respectively):
  639   \medskip
  640   \optitem[\.{--inputfile}, \.{-i} $\langle${\it input filename}$\rangle$]%
  641     {Specifies the raw, unfiltered input data to be crunched. The input file
  642     should describe the input as two columns containing $x$- and $y$-coordinates
  643     of the samples.}
  644   \optitem[\.{--outputfile}, \.{-o} $\langle${\it output filename}$\rangle$]%
  645     {Specifies the output file to which the program should write the filtered
  646     profile, again in a two-column format containing $x$- and $y$-coordinates
  647     of the filtered samples. If this option is omitted, the generated
  648     filtered data will instead be written to the console (terminal).}
  649   \optitem[\.{-nl} $\langle n_L\rangle$]%
  650     {Specifies the number of samples $n_L$ to use to the ``left'' of the basis
  651     sample in the regression window (kernel). The total number of samples in
  652     the window will be $nL+nR+1$.}
  653   \optitem[\.{-nr} $\langle n_R\rangle$]%
  654     {Specifies the number of samples $n_R$ to use to the ``right'' of the basis
  655     sample in the regression window (kernel). The total number of samples in
  656     the window will be $nL+nR+1$.}
  657   \optitem[\.{-m} $\langle m\rangle$]%
  658     {Specifies the order $m$ of the polynomial
  659     $p(x)=a_0+a_1(x-x_k)+a_2(x-x_k)^2+\ldots+a_m (x-x_k)^m$ to use in the
  660     regression analysis leading to the Savitzky--Golay coefficients.
  661     Typical values are between $m=2$ and $m=6$. Beware of too high values,
  662     which easily makes the regression too sensitive, with an oscillatory
  663     result.}
  664   \optitem[\.{-ld} $\langle l_D\rangle$]%
  665     {Specifies the order of the derivative to extract from the Savitzky--Golay
  666     smoothing algorithm. For regular Savitzky--Golay smoothing of the input
  667     data as such, use $l_D=0$. For the Savitzky--Golay smoothing and extraction
  668     of derivatives, set $l_D$ to the order of the desired derivative and make
  669     sure that you correctly interpret the scaling parameters as described in
  670     {\sl Numerical Recipes in C}, 2nd Edn~(Cambridge University Press, New
  671     York, 1994).}
  672   \optitem[\.{--help}, \.{-h}]%
  673     {Displays a brief help message and terminates the \sgfilter\ program clean
  674     from any error codes.}
  675   \optitem[\.{--verbose}, \.{-v}]%
  676     {Toggle verbose mode. (Default: Off.) This option should always be omitted
  677     whenever no output file has been specified (that is to say, omit any
  678     \.{--verbose} or \.{-v} option whenever \.{--outputfile} or \.{-o} has
  679     been omitted), as the verbose logging otherwise will contaminate the
  680     filtered data stream written to the console (terminal).}
  681   
  682   @*Example of Savitzky--Golay filtering with the \sgfilter\ program.
  683   It is always a good idea to create a rather ``nasty'' test case for an
  684   algorithm, and in this case it also provides the reader of this code with
  685   a (hopefully) clear example of usage of the \sgfilter\ program.
  686   
  687   We start off with generating a test suite of noisy data, in this case
  688   modeling an underlying function $g(x)$ with superimposed Gaussian noise as
  689   $$
  690     f(x)=\underbrace{
  691            \cos(3x)\sin^2(x^3)+4\sum^4_{k=1}\exp(-(x-x_k)^2/w^2_k)
  692          }_{g(x)} + \underbrace{Vu(x)\vphantom{\sum^4_{k=1}}}_{\rm noise}
  693   $$
  694   where $u(x)$ is a normally distributed stochastic variable of mean zero and
  695   unit variance, $V$ is the local variance as specified arbitrarily, and the
  696   remaining parameters $(x_k,w_k)$ are the positions and widths of four
  697   Gaussian peaks superimposed onto the otherwise trigonometric expression
  698   for the underlying function. For the current test suite we use
  699   $(x_1,w_1)=(0.2,0.007)$, $(x_2,w_2)=(0.4,0.01)$, $(x_3,w_3)=(0.6,0.02)$,
  700   and $(x_4,w_4)=(0.8,0.04)$. These Gaussian peaks serve as to provide various
  701   degrees of rapidly varying data, to check the performance in finding maxima.
  702   Meanwhile, the less rapidly varying domains which are dominated by the
  703   trigonometric expression serves as a test for the capability of the filter
  704   to handle rather moderate variations of low amplitude.
  705   
  706   The underlying test function $g(x)$ is shown\footnote{$\dagger$}{See the
  707   \.{*.eps} blocks in the enclosed \.{Makefile} for details on how \METAPOST
  708   was used in the generation of the encapsulated PostScript images shown in
  709   Figs.~1--6.} in Fig.~1.
  710   \medskip
  711   \centerline{\epsfbox{fig1.1}}
  712   \figcaption[1]{The underlying test function $g(x)=\cos(3x)\sin^2(x^3)%
  713   +4\sum^4_{k=1}\exp(-(x-x_k)^2/w^2_k)$ without any added noise.
  714   Here the positions and widths of the Gaussian peaks are $(x_1,w_1)=(0.2,0.007)$,
  715   $(x_2,w_2)=(0.4,0.01)$, $(x_3,w_3)=(0.6,0.02)$, and $(x_4,w_4)=(0.8,0.04)$.}
  716   \medskip
  717   \noindent
  718   The generator of artificial test data to be tested by the Savitsky--Golay
  719   filtering algorithm is here included as a simple \CEE\ program, which
  720   automatically is generated from the master \CWEB\ source when passed through
  721   \CTANGLE.
  722   
  723   @(example.c@>=
  724   #include <stdio.h>
  725   #include <stdlib.h>
  726   #include <time.h>
  727   #include <math.h>
  728   #define TWOPI (2.0*3.141592653589793)
  729   @|
  730   float gauss(float x, float w, float xa) {
  731      return(exp(-pow(((x-xa)/w),2)));
  732   }
  733   @|
  734   float func(float x) {
  735      float retval=gauss(x,0.007,0.2); /* $x_1=0.2$, $w_1=0.007$ */
  736      retval+=gauss(x,0.01,0.4);       /* $x_2=0.4$, $w_2=0.01$ */
  737      retval+=gauss(x,0.02,0.6);       /* $x_3=0.6$, $w_3=0.02$ */
  738      retval+=gauss(x,0.04,0.8);       /* $x_4=0.8$, $w_4=0.04$ */
  739      retval*=4.0;
  740      retval+=cos(3.0*x)*pow(sin(pow(x,3)),2);
  741      return(retval);
  742   }
  743   @|
  744   int main(int argc, char *argv[]) {
  745      int k, mm=1024;
  746      float var=1.0,xmax=2.5, x1, x2, u, v, f, z;
  747      if (argc>1) sscanf(argv[1],"%f",&var); /* Read first argument as variance */
  748      srand((unsigned)time(NULL)); /* Initialize random number generator */
  749      for (k=0;k<mm-1;k+=2) {
  750         x1=xmax*k/((double)mm-1);
  751         x2=xmax*(k+1)/((double)mm-1);
  752         u=((float)rand())/RAND_MAX; /* Uniformly distributed over $[0,1]$ */
  753         v=((float)rand())/RAND_MAX; /* Uniformly distributed over $[0,1]$ */
  754         if (u>0.0) { /* Apply the Box--Muller algorithm on |u| and |v| */
  755            f=sqrt(-2*log(u));
  756            z=TWOPI*v;
  757            u=f*cos(z); /* Normally distributed with E(|u|)=0 and Var(|u|)=1 */
  758            v=f*sin(z); /* Normally distributed with E(|u|)=0 and Var(|u|)=1 */
  759            fprintf(stdout,"%1.8f %1.8f\n",x1,func(x1)+var*u); /* $f(x_1)$ */
  760            fprintf(stdout,"%1.8f %1.8f\n",x2,func(x2)+var*v); /* $f(x_2)$ */
  761         }
  762      }
  763      return(0);
  764   }
  765   
  766   @ After having compiled the above code \.{example.c}, simply run
  767   \.{./example }$\langle${\it noise variance}$\rangle$\.{ > }$\langle$%
  768   {\it file name}$\rangle$ in order to generate the test function with a
  769   superimposed normally distributed (gaussian) noise of desired variance.
  770   In particular, we will in this test suite consider variances of 0 (that is to
  771   say, the underlying function without any noise), 0.5, 1.0, and 2.0.
  772   Such data files are simply generated by executing\footnote{$\dagger$}{The
  773   resulting test data, which so far has not been subject to any filtering,
  774   may easily be viewed by running the following script in Octave/Matlab:
  775   \medskip
  776   {\obeylines\obeyspaces\tt
  777   clear all; close all;
  778   hold on;
  779   u=load('example-2.0.dat'); plot(u(:,1),u(:,2),'-b');
  780   u=load('example-1.0.dat'); plot(u(:,1),u(:,2),'-c');
  781   u=load('example-0.5.dat'); plot(u(:,1),u(:,2),'-r');
  782   u=load('example-0.0.dat'); plot(u(:,1),u(:,2),'-k');
  783   legend('var(f(x))=2.0','var(f(x))=1.0','var(f(x))=0.5','var(f(x))=0.0');
  784   hold off;
  785   title('Artificial data generated for tests of Savitzky-Golay filtering');
  786   xlabel('x');
  787   ylabel('f(x)');
  788   }
  789   }
  790   \bigskip
  791   {\obeylines\obeyspaces\tt
  792   ./example 0.0 > example-0.0.dat
  793   ./example 0.5 > example-0.5.dat
  794   ./example 1.0 > example-1.0.dat
  795   ./example 2.0 > example-2.0.dat
  796   }
  797   \bigskip
  798   \noindent
  799   The resulting ``noisified'' suite of test data are in Figs.~2--4 shown for the
  800   respective noise variances of $V=0.5$, $V=1.0$, and $V=2.0$, respectively.
  801   \medskip
  802   \centerline{\epsfbox{fig2.1}}
  803   \figcaption[2]{The test function $g(x)$ with added Gaussian noise of variance
  804     $V=2.0$, as stored in file \.{example-2.0.dat} in the test suite.}
  805   
  806   @ Applying the Savitzky--Golay filter to the test data is now a straightforward
  807   task.
  808   Say that we wish to test filtering with polynomial degree $|m|=4$ and $|ld|=0$
  809   (which is the default value of |ld|, for regular smoothing with a delivered
  810   ``zero:th order derivative'', that is to say the smoothed non-differentiated
  811   function), for the two cases
  812   $|nl|=|nr|=10$ (in total 21 points in the regression kernel) and
  813   $|nl|=|nr|=60$ (in total 121 points in the regression kernel).
  814   Using the previously generated test suite of noisy data, the filtering is then
  815   easily accomplished by executing:
  816   \bigskip
  817   {\obeylines\obeyspaces\tt
  818   ./sgfilter -m 4 -nl 60 -nr 60 -i example-0.0.dat -o example-0.0-f-60.dat
  819   ./sgfilter -m 4 -nl 10 -nr 10 -i example-0.0.dat -o example-0.0-f-10.dat
  820   ./sgfilter -m 4 -nl 60 -nr 60 -i example-0.5.dat -o example-0.5-f-60.dat
  821   ./sgfilter -m 4 -nl 10 -nr 10 -i example-0.5.dat -o example-0.5-f-10.dat
  822   ./sgfilter -m 4 -nl 60 -nr 60 -i example-1.0.dat -o example-1.0-f-60.dat
  823   ./sgfilter -m 4 -nl 10 -nr 10 -i example-1.0.dat -o example-1.0-f-10.dat
  824   ./sgfilter -m 4 -nl 60 -nr 60 -i example-2.0.dat -o example-2.0-f-60.dat
  825   ./sgfilter -m 4 -nl 10 -nr 10 -i example-2.0.dat -o example-2.0-f-10.dat
  826   }
  827   \bigskip
  828   \noindent
  829   The resulting filtered data sets are shown in Figs.~3--6 for noise variances
  830   of $V=2.0$, $V=1.0$, $V=0.5$, and $V=0$, respectively. The final case
  831   corresponds to the interesting case of filtering the underlying function
  832   $g(x)$ without any added noise whatsoever, which corresponds to performing
  833   a local regression of the regression polynomial to the analytical trigonometric
  834   and exponential functions being terms of $g(x)$, a regression where we by no
  835   means should expect a perfect match.
  836   \medskip
  837   \centerline{\epsfbox{fig3.1}}
  838   \figcaption[3]{The profiles resulting from Savitzky--Golay-filtering of the
  839     test function $g(x)$ with added Gaussian noise of variance
  840     $V=2.0$.}
  841   \medskip
  842   As can be seen in Fig.~3, the results from filtering the ``worst-case'' set
  843   (with noise variance $V=2.0$) with $|nl|=|nr|=10$ ($|nl|+|nr|+1=21$ samples
  844   in the regression kernel) and $|m|=4$ (curve in blue) yield a rather good
  845   tracking of the narrow Gaussian peaks, meanwhile performing rather poor in
  846   the low-amplitude and rather slowly varying trigonometric hills in the
  847   right-hand side of the graph.
  848   As a reference, the underlying function $g(x)$ is mapped in dashed black.
  849   On the other hand, with $|nl|=|nr|=60$ ($|nl|+|nr|+1=121$ samples in the
  850   regression kernel) and keeping the same degree of the
  851   regression polynomial (curve in red), the narrow peaks are barely followed,
  852   meanwhile having a rather poor performance in the slowly varying hills as
  853   well (albeit of a very poor signal-to-noise ratio).
  854   \vfill\eject
  855   
  856   However, with a lower variance of the superimposed noise, we at $V=1$ have a
  857   nice tracking of the slowly varying hills by the 121-sample regression kernel,
  858   as shown in Fig.~4, meanwhile having a good tracking of the narrow Gaussian
  859   peaks by the 21-sample regression kernel.
  860   That the $|nl|+||=21$-sample window is noisier than the 121-sample one should
  861   not really be surprising, as the higher number of samples in the regression
  862   window tend to smoothen out the regression even further.
  863   
  864   \centerline{\epsfbox{fig4.1}}
  865   \figcaption[4]{The profiles resulting from Savitzky--Golay-filtering of the
  866     test function $g(x)$ with added Gaussian noise of variance
  867     $V=1.0$.}
  868   \medskip
  869   \centerline{\epsfbox{fig5.1}}
  870   \figcaption[5]{The profiles resulting from Savitzky--Golay-filtering of the
  871     test function $g(x)$ with added Gaussian noise of variance
  872     $V=0.5$.}
  873   \vfill\eject
  874   
  875   \centerline{\epsfbox{fig6.1}}
  876   \figcaption[6]{The profiles resulting from Savitzky--Golay-filtering of the
  877     test function $g(x)$ with added Gaussian noise of variance
  878     $V=0$, that is to say a direct regression against the underlying function
  879     $g(x)$. This figure illustrates the limitations of the attempts of
  880     linear regression of a polynomial to an underlying function which
  881     clearly cannot be approximated by a simple polynomial expressin in
  882     certain domains. One should always keep this limitation in mind before
  883     accepting (or discarding) data filtered by the Savitzky--Golay algorithm,
  884     despite the many cases in which it performs exceptionally well.}
  885   \vfill\eject
  886   
  887   @*Header files. We put the interface part of our routines in the header file
  888   \.{sgfilter.h}, and we restrict ourselves to a lowercase file name to
  889   maintain portability among operating systems with case-insensitive file names.
  890   There are just a few global definitions present in the \sgfilter\ program:\par
  891   \varitem[|VERSION|]{The current program revision number.}
  892   \varitem[|COPYRIGHT|]{The copyright banner.}
  893   \varitem[|DEFAULT_NL|]{The default value to use for the |nl| variable unless
  894     stated otherwise at the command line during startup of the program.
  895     This parameter specifies the number of samples $n_L$ to use to the ``left''
  896     of the basis sample in the regression window (kernel). The total number of
  897     samples in the window will be $nL+nR+1$.}
  898   \varitem[|DEFAULT_NR|]{The default value to use for the |nr| variable unless
  899     stated otherwise at the command line during startup of the program.
  900     This parameter specifies the number of samples $n_R$ to use to the ``right''
  901     of the basis sample in the regression window (kernel). The total number of
  902     samples in the window will be $nL+nR+1$.}
  903   \varitem[|DEFAULT_M|]{The default value to use for the |m| variable unless
  904     stated otherwise at the command line during startup of the program.
  905     This parameter specifies the order $m$ of the polynomial
  906     $p(x)=a_0+a_1(x-x_k)+a_2(x-x_k)^2+\ldots+a_m (x-x_k)^m$ to use in the
  907     regression analysis leading to the Savitzky--Golay coefficients.
  908     Typical values are between $m=2$ and $m=6$. Beware of too high values,
  909     which easily makes the regression too sensitive, with an oscillatory
  910     result.}
  911   \varitem[|DEFAULT_LD|]{The default value to use for the |ld| variable unless
  912     stated otherwise at the command line during startup of the program.
  913     This parameter specifies the order of the derivative to extract from the
  914     Savitzky--Golay smoothing algorithm.
  915     For regular Savitzky--Golay smoothing of the input
  916     data as such, use $l_D=0$. For the Savitzky--Golay smoothing and extraction
  917     of derivatives, set $l_D$ to the order of the desired derivative and make
  918     sure that you correctly interpret the scaling parameters as described in
  919     {\sl Numerical Recipes in C}, 2nd Edn~(Cambridge University Press, New
  920     York, 1994).}
  921   \varitem[|NCHMAX|]{The maximum number of characters allowed in strings for
  922     storing file names, including path.}
  923   \varitem[|log(...)|]{The |log()| macro forms the user interface to run-time
  924     error messaging and console logging activities, by invoking the
  925     |log_printf()| routine. The |log()| macro is preferrably used rather than
  926     direct calls to |log_printf()|, as it automatizes the extraction of the
  927     calling routine, via the |__func__| macro. Notice that the clean construction
  928     of the macro using variable-length arguments (via |__VA_ARGS__|) {\sl only} is
  929     supported at ISO 9899:1999 (ISO C99) standard level or above (see, for
  930     example, \.{http://en.wikipedia.org/wiki/C99}), which is the default
  931     compliance level used in the enclosed \.{Makefile}. (In fact, this is here
  932     {\sl the} very reason for using ISO 9899:1999 rather than ISO 9899:1990).}
  933   
  934   @(sgfilter.h@>=
  935   #define VERSION "1.6"
  936   #define COPYRIGHT "Copyright (C) 2006-2011, Fredrik Jonsson"
  937   #define DEFAULT_NL (15)
  938   #define DEFAULT_NR (15)
  939   #define DEFAULT_M (4)
  940   #define DEFAULT_LD (0)
  941   #define EPSILON ((double)(1.0e-20))
  942   #define NCHMAX (256)
  943   #define CONVOLVE_WITH_NR_CONVLV (0)
  944   #define log(...) log_printf(__func__, __VA_ARGS__)
  945   
  946   @|
  947   #if defined(_CYGWIN_SIGNAL_H)||defined(__APPLE__)||defined(__unix__)||defined(__linux)
  948   #define UNIX_LIKE_OS (1)
  949   #endif
  950   
  951   @|
  952   void log_printf(const char *function_name, const char *format, ...);
  953   int *ivector(long nl, long nh);
  954   double *dvector(long nl, long nh);
  955   double **dmatrix(long nrl, long nrh, long ncl, long nch);
  956   void free_ivector(int *v, long nl, long nh);
  957   void free_dvector(double *v, long nl, long nh);
  958   void free_dmatrix(double **m,long nrl,long nrh,long ncl,long nch);
  959   void lubksb(double **a, int n, int *indx, double b[]);
  960   void ludcmp(double **a, int n, int *indx, double *d);
  961   void four1(double data[], unsigned long nn, int isign);
  962   void twofft(double data1[], double data2[], double fft1[], double fft2[],
  963      unsigned long n);
  964   void realft(double data[], unsigned long n, int isign);
  965   char convlv(double data[], unsigned long n, double respns[], unsigned long m,
  966               int isign, double ans[]);
  967   char sgcoeff(double c[], int np, int nl, int nr, int ld, int m);
  968   char sgfilter(double yr[], double yf[], int mm, int nl, int nr, int ld, int m);
  969   char *strip_away_path(char filename[]);
  970   long int num_coordinate_pairs(FILE *file);
  971   
  972   @*The main program. Here follows the general outline of the |main| routine of
  973   the \sgfilter\ program. This is where it all starts.
  974   
  975   @c
  976   @<Library inclusions@>@;
  977   @<Global variables@>@;
  978   @<Definitions of routines@>@;
  979   
  980   int main(int argc, char *argv[])
  981   {
  982      @<Definition of variables@>@;
  983      @<Parse command line for options and parameters@>@;
  984      @<Allocate memory and read |M| samples of unfiltered data from file@>@;
  985      @<Filter raw data through the Savitzky--Golay smoothing filter@>@;
  986      @<Write filtered data to file or terminal output@>@;
  987      @<Deallocate memory@>@;
  988      return(EXIT_SUCCESS);
  989   }
  990   
  991   @ Library dependencies.
  992   The standard \ANSICEE\ libraries included in this program are:\par
  993   \varitem[\.{math.h}]{For access to common mathematical functions.}
  994   \varitem[\.{stdio.h}]{For file access and any block involving |fprintf|.}
  995   \varitem[\.{stdlib.h}]{For access to the |exit| function.}
  996   \varitem[\.{stdarg.h}]{For access to |vsprintf|, |vfprintf| etc.}
  997   \varitem[\.{string.h}]{For string manipulation, |strcpy|, |strcmp| etc.}
  998   \varitem[\.{ctype.h}]{For access to the |isalnum| function.}
  999   \varitem[\.{time.h}]{For access to |ctime|, |clock| and time stamps.}
 1000   \varitem[\.{sys/time.h}]{For access to millisecond-resolved timer.}
 1001   
 1002   @<Library inclusions@>=
 1003   #include <math.h>
 1004   #include <stdlib.h>
 1005   #include <stdarg.h>
 1006   #include <stdio.h>
 1007   #include <string.h>
 1008   #include <ctype.h>
 1009   #include <time.h>
 1010   #include <sys/time.h>
 1011   #include "sgfilter.h"
 1012   
 1013   @ Declaration of global variables. The only global variables allowed in
 1014   my programs are |optarg|, which is a pointer to the the string of characters
 1015   that specified the call from the command line, and |progname|, which simply
 1016   is a pointer to the string containing the name of the program, as it was
 1017   invoked from the command line.
 1018   Typically, the |progname| string is the result remaining after having passed
 1019   |argv[0]| through the |strip_away_path| routine, just before parsing the
 1020   command line for parameters.
 1021   
 1022   @<Global variables@>=
 1023   extern char *optarg;
 1024   char *progname;
 1025   
 1026   @ Declaration of local variables of the |main| program.
 1027   In \CWEB\ one has the option of adding variables along the program, for example
 1028   by locally adding temporary variables related to a given sub-block of code.
 1029   However, the philosophy in the \sgfilter\ program is to keep all variables of
 1030   the |main| section collected together, so as to simplify tasks as, for example,
 1031   tracking down a given variable type definition.
 1032   The local variables of the program are as follows:
 1033   \medskip
 1034   
 1035   \varitem[|*x|, |*yr|, |*yf|]{Pointers to the double precision vectors keeping
 1036     the abscissa, unfiltered (raw) ordinata, and the and the resulting filtered
 1037     ordinata, respectively.}
 1038   \varitem[|mm|]{Keeps track of the number of samples to analyze, $|mm|\equiv M$.}
 1039   \varitem[|*file|]{Dummy file pointer used for scanning through input data and
 1040     for writing data to an output file.}
 1041   \varitem[|input_filename|]{String for keeping the file name of the input data.}
 1042   \varitem[|output_filename|]{Ditto for the output data.}
 1043   \varitem[|verbose|]{Determines whether the program should run in verbose mode
 1044     (logging various activities along the execution) or not. Default is off.}
 1045   
 1046   @<Definition of variables@>=
 1047   int no_arg;
 1048   int nl=DEFAULT_NL;
 1049   int nr=DEFAULT_NR;
 1050   int ld=DEFAULT_LD;
 1051   int m=DEFAULT_M;
 1052   long int k, mm=0;
 1053   double *x, *yr, *yf;
 1054   char input_filename[NCHMAX]="",output_filename[NCHMAX]="";
 1055   char verbose=0;
 1056   FILE *file;
 1057   
 1058   @ Parsing command line options. All input parameters are passed to the
 1059   program through command line options to the \sgfilter\ program.
 1060   The syntax of the accepted command line options is listed whenever the
 1061   program is invoked without any options, or whenever any of the \.{--help}
 1062   or \.{-h} options are specified at startup.
 1063   
 1064   @<Parse command line for options and parameters@>=
 1065   progname=strip_away_path(argv[0]);
 1066   no_arg=argc;
 1067   while(--argc) {
 1068      if(!strcmp(argv[no_arg-argc],"-o") ||
 1069         !strcmp(argv[no_arg-argc],"--outputfile")) {
 1070         --argc;
 1071         strcpy(output_filename,argv[no_arg-argc]);
 1072      } else if(!strcmp(argv[no_arg-argc],"-i") ||
 1073         !strcmp(argv[no_arg-argc],"--inputfile")) {
 1074         --argc;
 1075         strcpy(input_filename,argv[no_arg-argc]);
 1076      } else if(!strcmp(argv[no_arg-argc],"-h") ||
 1077                !strcmp(argv[no_arg-argc],"--help")) {
 1078         showsomehelp();
 1079         exit(0);
 1080      } else if(!strcmp(argv[no_arg-argc],"-v") ||
 1081         !strcmp(argv[no_arg-argc],"--verbose")) {
 1082         verbose=(verbose?0:1);
 1083      } else if(!strcmp(argv[no_arg-argc],"-nl")) {
 1084         --argc;
 1085         if (!sscanf(argv[no_arg-argc],"%d",&nl)) {
 1086            log("Error in '-nl' option.");
 1087            exit(1);
 1088         }
 1089      } else if(!strcmp(argv[no_arg-argc],"-nr")) {
 1090         --argc;
 1091         if (!sscanf(argv[no_arg-argc],"%d",&nr)) {
 1092            log("Error in '-nr' option.");
 1093            exit(1);
 1094         }
 1095      } else if(!strcmp(argv[no_arg-argc],"-ld")) {
 1096         --argc;
 1097         if (!sscanf(argv[no_arg-argc],"%d",&ld)) {
 1098            log("Error in '-ld' option.");
 1099            exit(1);
 1100         }
 1101      } else if(!strcmp(argv[no_arg-argc],"-m")) {
 1102         --argc;
 1103         if (!sscanf(argv[no_arg-argc],"%d",&m)) {
 1104            log("Error in '-m' option.");
 1105            exit(1);
 1106         }
 1107      } else {
 1108         log("Unrecognized option '%s'.",argv[no_arg-argc]);
 1109         showsomehelp();
 1110         exit(1);
 1111      }
 1112   }
 1113   if (verbose) fprintf(stdout,"This is %s v.%s. %s\n",progname,VERSION,COPYRIGHT);
 1114   
 1115   @ This is the where the raw (unfiltered) data is loaded into the computer
 1116   memory. This block assumes that the number of data points, $M$, has {\sl not}
 1117   previously been determined, neither by analysis of the input file nor
 1118   explicitly stated via command line options, and hence the number of samples
 1119   is automatically extracted using the |num_coordinate_pairs| routine.
 1120   
 1121   @<Allocate memory and read |M| samples of unfiltered data from file@>=
 1122   if (!strcmp(input_filename,"")) {
 1123      log("No input file specified! (Please use the '--inputfile' option.)");
 1124      log("Execute '%s --help' for help.",progname);
 1125      exit(1);
 1126   }
 1127   if ((file=fopen(input_filename,"r"))==NULL) {
 1128      log("Could not open %s for loading raw data!",input_filename);
 1129      exit(1);
 1130   }
 1131   mm=num_coordinate_pairs(file);
 1132   if (mm<nl+nr+1) {
 1133      log("Error: The number M=%ld of data points must be at least nl+nr+1=%d",
 1134         mm, nl+nr+1);
 1135      log("Please check your -nl or -nr options.");
 1136      exit(1);
 1137   }
 1138   if (verbose) {
 1139      log("Loading %ld unfiltered samples from %s...", mm, input_filename);
 1140      log(" ... allocating memory for storage ...");
 1141   }
 1142   x=dvector(1,mm);
 1143   yr=dvector(1,mm);
 1144   #if CONVOLVE_WITH_NR_CONVLV
 1145   yf=dvector(1,2*mm);
 1146   #else
 1147   yf=dvector(1,mm);
 1148   #endif
 1149   if (verbose)
 1150      log(" ... scanning %s for input data ...",input_filename);
 1151   for (k=1;k<=mm;k++) {
 1152      fscanf(file,"%lf",&x[k]); /* Scan $x$-coordinate */
 1153      fscanf(file,"%lf",&yr[k]); /* Scan unfiltered $y$-coordinate */
 1154   }
 1155   fclose(file);
 1156   if (verbose)
 1157      log(" ... done. Input now residing in RAM.");
 1158   
 1159   @ Filter data. This is simple. One single call to the |sgfilter| routine and
 1160   we are done.
 1161   
 1162   @<Filter raw data through the Savitzky--Golay smoothing filter@>=
 1163   if (!sgfilter(yr, yf, mm, nl, nr, ld, m)) {
 1164       if (verbose)
 1165         log("Successfully performed Savitzky-Golay filtering.");
 1166   } else {
 1167      if (verbose)
 1168         log("Error: Could not perform Savitzky-Golay filtering.");
 1169   }
 1170   
 1171   @ Write the filtered data to file or |stdout|, depending on whether or not the
 1172   command line options \.{-o} or \.{-outputfile} were present when starting the
 1173   \sgfilter\ program. We here use a redirection of the |stdout| stream using
 1174   |freopen(output_filename, "w", stdout)) == NULL)| and (in the case of unix-like
 1175   systems) a re-redirection back to terminal output again with
 1176   |freopen("/dev/tty", "a", stdout)|.
 1177   
 1178   @<Write filtered data to file or terminal output@>=
 1179   if (!strcmp(output_filename,"")) { /* No filename specified */
 1180      if (verbose)
 1181         log("Writing %ld filtered samples to console...", mm);
 1182   } else { /* If file name specified $\Rightarrow$ redirect |stdout| to file */
 1183      if (verbose)
 1184         log("Writing %ld filtered samples to %s...", mm, output_filename);
 1185      if((file = freopen(output_filename, "w", stdout)) == NULL) {
 1186         log("Error: Unable to redirect stdout stream to file %s.",
 1187            output_filename);
 1188         exit(1);
 1189      }
 1190   }
 1191   for (k=1;k<=mm;k++) fprintf(stdout,"%1.8f %1.8f\n",x[k],yf[k]);
 1192   #ifdef UNIX_LIKE_OS
 1193   freopen("/dev/tty", "a", stdout); /* Redirect |stdout| back to console output */
 1194   #endif
 1195   if (verbose) log(" ... done.");
 1196   
 1197   @ Deallocate memory.
 1198   
 1199   @<Deallocate memory@>=
 1200   free_dvector(x,1,mm);
 1201   free_dvector(yr,1,mm);
 1202   #if CONVOLVE_WITH_NR_CONVLV
 1203   free_dvector(yf,1,2*mm);
 1204   #else
 1205   free_dvector(yf,1,mm);
 1206   #endif
 1207   
 1208   @*Routines used by the program.
 1209   
 1210   @<Definitions of routines@>=
 1211      @<Routine for error messaging@>@;
 1212      @<Routine for allocation of integer precision vectors@>@;
 1213      @<Routine for allocation of double floating-point precision vectors@>@;
 1214      @<Routine for allocation of double floating-point precision matrices@>@;
 1215      @<Routine for deallocation of integer precision vectors@>@;
 1216      @<Routine for deallocation of double floating-point precision vectors@>@;
 1217      @<Routine for deallocation of double floating-point precision matrices@>@;
 1218      @<Routine for solving systems of linear equations by LU decomposition@>@;
 1219      @<Routine for performing LU decomposition of a matrix@>@;
 1220      @<Routine for discrete Fourier transformation of real-valued data@>@;
 1221      @<Routine for simultaneous fast Fourier transformation of two data sets@>@;
 1222      @<Routine for Fourier transformation of real-valued data@>@;
 1223      @<Routine for numerical convolution@>@;
 1224      @<Routine for computation of coefficients for Savitzky--Golay filtering@>@;
 1225      @<Routine for Savitzky--Golay filtering@>@;
 1226      @<Routine for removing preceding path of filenames@>@;
 1227      @<Routine for displaying a brief help message on usage@>@;
 1228      @<Routine for obtaining the number of coordinate pairs in a file@>@;
 1229   
 1230   @ The |void log_printf(const char *function_name,const char *format, ...)|
 1231   routine writes formatted entries to standard output, displaying time and
 1232   calling routine in a coherent manner. Notice that although the |log_printf()|
 1233   routine is the one which performs the actual messaging, the |log()| macro
 1234   (defined in the header file) is the preferred way of accessing this routine,
 1235   as it provides a more compact notation and automatically takes care of supplying
 1236   the reference to the name of the calling function.
 1237   
 1238   Also notice that the |const char| type of the last two input pointer arguments
 1239   here is absolutely essential in order to pass strict pedantic compilation with
 1240   \GCC.
 1241   
 1242   The routine accepts two input parameters. First, |function_name| which should
 1243   be the name of the calling function. This is to ensure that any displayed
 1244   error messages are properly matched to the issuing routines. Notice, however,
 1245   that the |log()| macro (which is the preferred way of displaying error messages)
 1246   automatically takes care of supplying the proper function name.
 1247   Second, |format|, which simply is the format and message string to be
 1248   displayed, formatted in the \CEE-standard |printf()| or |fprintf()| syntax.
 1249   
 1250   @<Routine for error messaging@>=
 1251   void log_printf(const char *function_name, const char *format, ...) {
 1252      va_list args;
 1253      time_t time0;
 1254      struct tm lt;
 1255      struct timeval tv;
 1256      char logentry[1024];
 1257   
 1258      gettimeofday(&tv, NULL);
 1259      time(&time0);
 1260      lt = *localtime(&time0);
 1261      sprintf(logentry, "%02u%02u%02u %02u:%02u:%02u.%03d ",
 1262         lt.tm_year-100, lt.tm_mon+1, lt.tm_mday,
 1263         lt.tm_hour, lt.tm_min, lt.tm_sec, tv.tv_usec/1000);
 1264      sprintf(logentry+strlen(logentry),"(%s) ",function_name);
 1265      va_start(args, format); /* Initialize args by the |va_start()| macro */
 1266      vsprintf(logentry+strlen(logentry), format, args);
 1267      va_end(args); /* Terminate the use of args by the |va_end()| macro */
 1268      sprintf(logentry+strlen(logentry), "\n"); /* Always append newline */
 1269      fprintf(stdout, "%s", logentry);
 1270      return;
 1271   }
 1272   
 1273   @ The |int *ivector(long nl, long nh)| routine allocates a real-valued vector
 1274   of integer precision, with vector index ranging from |nl| to |nh|.
 1275   
 1276   @<Routine for allocation of integer precision vectors@>=
 1277   int *ivector(long nl, long nh) {
 1278      int *v;
 1279      v=(int *)malloc((size_t) ((nh-nl+2)*sizeof(int)));
 1280      if (!v) {
 1281         log("Error: Allocation failure.");
 1282         exit(1);
 1283      }
 1284      return v-nl+1;
 1285   }
 1286   
 1287   @ The |double *dvector(long nl, long nh)| routine allocates a real-valued
 1288   vector of double precision, with vector index ranging from |nl| to |nh|.
 1289   
 1290   @<Routine for allocation of double floating-point precision vectors@>=
 1291   double *dvector(long nl, long nh) {
 1292      double *v;
 1293      long k;
 1294      v=(double *)malloc((size_t) ((nh-nl+2)*sizeof(double)));
 1295      if (!v) {
 1296         log("Error: Allocation failure.");
 1297         exit(1);
 1298      }
 1299      for (k=nl;k<=nh;k++) v[k]=0.0;
 1300      return v-nl+1;
 1301   }
 1302   
 1303   @ The |double **dmatrix(long nrl, long nrh, long ncl, long nch)| routine
 1304   allocates an array of double floating-point precision, with row index
 1305   ranging from |nrl| to |nrh| and column index ranging from |ncl| to |nch|.
 1306   
 1307   @<Routine for allocation of double floating-point precision matrices@>=
 1308   double **dmatrix(long nrl, long nrh, long ncl, long nch)
 1309   {
 1310      long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
 1311      double **m;
 1312      m=(double **) malloc((size_t)((nrow+1)*sizeof(double*)));
 1313      if (!m) {
 1314         log("Allocation failure 1 occurred.");
 1315         exit(1);
 1316      }
 1317      m += 1;
 1318      m -= nrl;
 1319      m[nrl]=(double *) malloc((size_t)((nrow*ncol+1)*sizeof(double)));
 1320      if (!m[nrl]) {
 1321         log("Allocation failure 2 occurred.");
 1322         exit(1);
 1323      }
 1324      m[nrl] += 1;
 1325      m[nrl] -= ncl;
 1326      for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
 1327      return m;
 1328   }
 1329   
 1330   @ The |void free_ivector(int *v, long nl, long nh)| routine release the
 1331   memory occupied by the real-valued vector |v[nl...nh]|.
 1332   
 1333   @<Routine for deallocation of integer precision vectors@>=
 1334   void free_ivector(int *v, long nl, long nh) {
 1335      free((char*) (v+nl-1));
 1336   }
 1337   
 1338   @ The |free_dvector| routine release the memory occupied by the
 1339   real-valued vector |v[nl...nh]|.
 1340   
 1341   @<Routine for deallocation of double floating-point precision vectors@>=
 1342   void free_dvector(double *v, long nl, long nh) {
 1343      free((char*) (v+nl-1));
 1344   }
 1345   
 1346   @ The |void free_dmatrix(double **m,long nrl,long nrh,long ncl,long nch)|
 1347   routine releases the memory occupied by the double floating-point precision
 1348   matrix |v[nl...nh]|, as allocated by |dmatrix()|.
 1349   
 1350   @<Routine for deallocation of double floating-point precision matrices@>=
 1351   void free_dmatrix(double **m,long nrl,long nrh,long ncl,long nch) {
 1352      free((char*) (m[nrl]+ncl-1));
 1353      free((char*) (m+nrl-1));
 1354   }
 1355   
 1356   @ The |lubksb(double **a, int n, int *indx, double b[])| routine solves the set
 1357   of |n| linear equations ${\bf A}\cdot{\bf x}={\bf b}$ where ${\bf A}$ is a
 1358   real-valued $[n\times n]$-matrix and ${\bf x}$ and ${\bf b}$ are real-valued
 1359   $[n\times1]$-vectors. Here |a[1...n][1...n]| is input, however {\sl not} as the
 1360   matrix ${\bf A}$ but rather as its corresponding LU decomposition as determined
 1361   by the |ludcmp| routine. Here |indx[1...n]| is input as the permutation vector
 1362   returned by |ludcmp|, |b[1...n]| is input as the right-hand side vector
 1363   ${\bf b}$, and returns with the solution vector ${\bf x}$.
 1364   The parameters |a|, |n|, and |indx| are not modified by this routine and
 1365   can be left in place for successive calls with different right-hand sides
 1366   ${\bf b}$. This routine takes into account the possibility that ${\bf b}$
 1367   will begin with many zero elements, so it is efficient for use in matrix
 1368   inversion.
 1369   The |lubksb| routine is adopted from {\sl Numerical Recipes in C}, 2nd Edn
 1370   (Cambridge University Press, New York, 1994).
 1371   
 1372   @<Routine for solving systems of linear equations by LU decomposition@>=
 1373   void lubksb(double **a, int n, int *indx, double b[])
 1374   {
 1375      int i,ii=0,ip,j;
 1376      double sum;
 1377   
 1378      for (i=1;i<=n;i++) {
 1379         ip=indx[i];
 1380         sum=b[ip];
 1381         b[ip]=b[i];
 1382         if (ii)
 1383            for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j];
 1384         else if (sum) ii=i;
 1385         b[i]=sum;
 1386      }
 1387      for (i=n;i>=1;i--) {
 1388         sum=b[i];
 1389         for (j=i+1;j<=n;j++) sum -= a[i][j]*b[j];
 1390         b[i]=sum/a[i][i];
 1391      }
 1392   }
 1393   
 1394   @ Given a square and real-valued matrix |a[1...n][1...n]|, the
 1395   |ludcmp(float **a, int n, int *indx, float *d)| routine replaces it by the
 1396   corresponding LU decomposition of a rowwise permutation
 1397   of itself. Entering the routine, the square matrix |a| and its number of
 1398   columns (or rows) |n| are inputs. On return, |a| is arranged as in Eq.~(2.3.14)
 1399   of {\sl Numerical Recipes in C}, 2nd edition, while |indx[1...n]| is an output
 1400   vector that records the row permutation effected by the partial pivoting, and
 1401   |d| is output as $\pm1$ depending on whether the number of row interchanges
 1402   was even or odd, respectively. This routine is commonly used in combination
 1403   with |lubksb| to solve linear equations or invert a matrix.
 1404   The |ludcmp| routine is adopted from {\sl Numerical Recipes in C}, 2nd Edn
 1405   (Cambridge University Press, New York, 1994).
 1406   
 1407   @<Routine for performing LU decomposition of a matrix@>=
 1408   void ludcmp(double **a, int n, int *indx, double *d)
 1409   {
 1410      int i,imax=0,j,k;
 1411      double big,dum,sum,temp;
 1412      double *vv;
 1413   
 1414      vv=dvector(1,n);
 1415      *d=1.0;
 1416      for (i=1;i<=n;i++) {
 1417         big=0.0;
 1418         for (j=1;j<=n;j++)
 1419            if ((temp=fabs(a[i][j])) > big) big=temp;
 1420         if (big == 0.0) {
 1421            log("Error: Singular matrix found in routine ludcmp()");
 1422            exit(1);
 1423         }
 1424         vv[i]=1.0/big;
 1425      }
 1426      for (j=1;j<=n;j++) {
 1427         for (i=1;i<j;i++) {
 1428            sum=a[i][j];
 1429            for (k=1;k<i;k++) sum -= a[i][k]*a[k][j];
 1430            a[i][j]=sum;
 1431         }
 1432         big=0.0;
 1433         for (i=j;i<=n;i++) {
 1434            sum=a[i][j];
 1435            for (k=1;k<j;k++)
 1436               sum -= a[i][k]*a[k][j];
 1437            a[i][j]=sum;
 1438            if ( (dum=vv[i]*fabs(sum)) >= big) {
 1439               big=dum;
 1440               imax=i;
 1441            }
 1442         }
 1443         if (j != imax) {
 1444            for (k=1;k<=n;k++) {
 1445               dum=a[imax][k];
 1446               a[imax][k]=a[j][k];
 1447               a[j][k]=dum;
 1448            }
 1449            *d = -(*d);
 1450            vv[imax]=vv[j];
 1451         }
 1452         indx[j]=imax;
 1453         if (a[j][j] == 0.0) a[j][j]=EPSILON;
 1454         if (j != n) {
 1455            dum=1.0/(a[j][j]);
 1456            for (i=j+1;i<=n;i++) a[i][j] *= dum;
 1457         }
 1458      }
 1459      free_dvector(vv,1,n);
 1460   }
 1461   
 1462   @ The |four1(double data[], unsigned long nn, int isign)| routine replaces
 1463   |data[1...2*nn]| by its discrete Fourier transform, if |isign| is input
 1464   as $+1$; or replaces |data[1..2*nn]| by |nn| times its inverse discrete Fourier
 1465   transform, if |isign| is input as $-1$. Here |data| is a complex-valued array
 1466   of length |nn| or, equivalently, a real array of length |2*nn|, wher |nn|
 1467   {\sl must} be an integer power of 2 (this is not checked for).
 1468   The |four1| routine is adopted from {\sl Numerical Recipes in C}, 2nd Edn
 1469   (Cambridge University Press, New York, 1994).
 1470   
 1471   @<Routine for discrete Fourier transformation of real-valued data@>=
 1472   #include <math.h>
 1473   #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
 1474   void four1(double data[], unsigned long nn, int isign)
 1475   {
 1476      unsigned long n,mmax,m,j,istep,i;
 1477      double wtemp,wr,wpr,wpi,wi,theta;
 1478      double tempr,tempi;
 1479   
 1480      n=nn << 1;
 1481      j=1;
 1482      for (i=1;i<n;i+=2) {
 1483         if (j > i) {
 1484            SWAP(data[j],data[i]);
 1485            SWAP(data[j+1],data[i+1]);
 1486         }
 1487         m=n >> 1;
 1488         while (m >= 2 && j > m) {
 1489            j -= m;
 1490            m >>= 1;
 1491         }
 1492         j += m;
 1493      }
 1494      mmax=2;
 1495      while (n > mmax) {
 1496         istep=mmax << 1;
 1497         theta=isign*(6.28318530717959/mmax);
 1498         wtemp=sin(0.5*theta);
 1499         wpr = -2.0*wtemp*wtemp;
 1500         wpi=sin(theta);
 1501         wr=1.0;
 1502         wi=0.0;
 1503         for (m=1;m<mmax;m+=2) {
 1504            for (i=m;i<=n;i+=istep) {
 1505               j=i+mmax;
 1506               tempr=wr*data[j]-wi*data[j+1];
 1507               tempi=wr*data[j+1]+wi*data[j];
 1508               data[j]=data[i]-tempr;
 1509               data[j+1]=data[i+1]-tempi;
 1510               data[i] += tempr;
 1511               data[i+1] += tempi;
 1512            }
 1513            wr=(wtemp=wr)*wpr-wi*wpi+wr;
 1514            wi=wi*wpr+wtemp*wpi+wi;
 1515         }
 1516         mmax=istep;
 1517      }
 1518   }
 1519   #undef SWAP
 1520   
 1521   
 1522   @ The |twofft(double data1[], double data2[], double fft1[], double fft2[],
 1523   unsigned long n)| routine takes two real-valued arrays |data1[1...n]|
 1524   and |data2[1...n]| as input, makes a call to the |four1| routine and
 1525   returns two complex-valued output arrays |fft1[1...2*n]|
 1526   and |fft2[1...2*n]|, each of complex length |n| (that is to say, a real
 1527   length of |2*n| elements), which contain the discrete Fourier transforms
 1528   of the respective data arrays. Here |n| {\sl must} be an integer power of 2.
 1529   The |twofft| routine is adopted from {\sl Numerical Recipes in C}, 2nd Edn
 1530   (Cambridge University Press, New York, 1994).
 1531   
 1532   @<Routine for simultaneous fast Fourier transformation of two data sets@>=
 1533   void twofft(double data1[], double data2[], double fft1[], double fft2[],
 1534      unsigned long n)
 1535   {
 1536      void four1(double data[], unsigned long nn, int isign);
 1537      unsigned long nn3,nn2,jj,j;
 1538      double rep,rem,aip,aim;
 1539   
 1540      nn3=1+(nn2=2+n+n);
 1541      for (j=1,jj=2;j<=n;j++,jj+=2) {
 1542         fft1[jj-1]=data1[j];
 1543         fft1[jj]=data2[j];
 1544      }
 1545      four1(fft1,n,1);
 1546      fft2[1]=fft1[2];
 1547      fft1[2]=fft2[2]=0.0;
 1548      for (j=3;j<=n+1;j+=2) {
 1549         rep=0.5*(fft1[j]+fft1[nn2-j]);
 1550         rem=0.5*(fft1[j]-fft1[nn2-j]);
 1551         aip=0.5*(fft1[j+1]+fft1[nn3-j]);
 1552         aim=0.5*(fft1[j+1]-fft1[nn3-j]);
 1553         fft1[j]=rep;
 1554         fft1[j+1]=aim;
 1555         fft1[nn2-j]=rep;
 1556         fft1[nn3-j] = -aim;
 1557         fft2[j]=aip;
 1558         fft2[j+1] = -rem;
 1559         fft2[nn2-j]=aip;
 1560         fft2[nn3-j]=rem;
 1561      }
 1562   }
 1563   
 1564   @ The |realft(double data[], unsigned long n, int isign)| routine
 1565   calculates the Fourier transform of a set of |n| real-valued data points.
 1566   On return the routine replaces this data (which is stored in array
 1567   |data[1...n]|) by the positive frequency half of its complex Fourier
 1568   transform. The real-valued first and last components of the complex
 1569   transform are returned in elements |data[1]| and |data[2]|, respectively.
 1570   Here |n| must be a power of 2. This routine also calculates the inverse
 1571   transform of a complex data array if it is the transform of real data.
 1572   (In this case, the result must be multiplied by |2/n|.)
 1573   The |realft| routine is adopted from {\sl Numerical Recipes in C}, 2nd Edn
 1574   (Cambridge University Press, New York, 1994).
 1575   
 1576   @<Routine for Fourier transformation of real-valued data@>=
 1577   void realft(double data[], unsigned long n, int isign)
 1578   {
 1579      void four1(double data[], unsigned long nn, int isign);
 1580      unsigned long i,i1,i2,i3,i4,np3;
 1581      double c1=0.5,c2,h1r,h1i,h2r,h2i;
 1582      double wr,wi,wpr,wpi,wtemp,theta;
 1583   
 1584      theta=3.141592653589793/(double) (n>>1);
 1585      if (isign == 1) {
 1586         c2 = -0.5;
 1587         four1(data,n>>1,1);
 1588      } else {
 1589         c2=0.5;
 1590         theta = -theta;
 1591      }
 1592      wtemp=sin(0.5*theta);
 1593      wpr = -2.0*wtemp*wtemp;
 1594      wpi=sin(theta);
 1595      wr=1.0+wpr;
 1596      wi=wpi;
 1597      np3=n+3;
 1598      for (i=2;i<=(n>>2);i++) {
 1599         i4=1+(i3=np3-(i2=1+(i1=i+i-1)));
 1600         h1r=c1*(data[i1]+data[i3]);
 1601         h1i=c1*(data[i2]-data[i4]);
 1602         h2r = -c2*(data[i2]+data[i4]);
 1603         h2i=c2*(data[i1]-data[i3]);
 1604         data[i1]=h1r+wr*h2r-wi*h2i;
 1605         data[i2]=h1i+wr*h2i+wi*h2r;
 1606         data[i3]=h1r-wr*h2r+wi*h2i;
 1607         data[i4] = -h1i+wr*h2i+wi*h2r;
 1608         wr=(wtemp=wr)*wpr-wi*wpi+wr;
 1609         wi=wi*wpr+wtemp*wpi+wi;
 1610      }
 1611      if (isign == 1) {
 1612         data[1] = (h1r=data[1])+data[2];
 1613         data[2] = h1r-data[2];
 1614      } else {
 1615         data[1]=c1*((h1r=data[1])+data[2]);
 1616         data[2]=c1*(h1r-data[2]);
 1617         four1(data,n>>1,-1);
 1618      }
 1619   }
 1620   
 1621   @ The |convlv(double data[], unsigned long n, double respns[],
 1622   unsigned long m, int isign, double ans[])| routine convolves or deconvolves
 1623   a real data set |data[1...n]| (including any user-supplied zero padding)
 1624   with a response function |respns[1...n]|. The response function must be
 1625   stored in wrap-around order in the first |m| elements of respns, where |m|
 1626   is an odd integer less than or equal to |n|. Wrap-around order means
 1627   that the first half of the array respns contains the impulse response
 1628   function at positive times, while the second half of the array contains
 1629   the impulse response function at negative times, counting down from the
 1630   highest element |respns[m]|. On input, |isign| is $+1$ for convolution, and
 1631   $-1$ for deconvolution. The answer is returned in the first |n| components
 1632   of |ans|. However, |ans| {\sl must} be supplied in the calling program with
 1633   dimensions |[1...2*n]|, for consistency with the |twofft| routine.
 1634   Here |n| {\sl must} be an integer power of two.
 1635   The |convlv| routine is adopted from {\sl Numerical Recipes in C}, 2nd Edn
 1636   (Cambridge University Press, New York, 1994).
 1637   
 1638   @<Routine for numerical convolution@>=
 1639   char convlv(double data[], unsigned long n, double respns[], unsigned long m,
 1640               int isign, double ans[])
 1641   {
 1642      void realft(double data[], unsigned long n, int isign);
 1643      void twofft(double data1[], double data2[], double fft1[], double fft2[],
 1644         unsigned long n);
 1645      unsigned long i,no2;
 1646      double dum,mag2,*fft;
 1647   
 1648      fft=dvector(1,n<<1);
 1649      for (i=1;i<=(m-1)/2;i++)
 1650         respns[n+1-i]=respns[m+1-i];
 1651      for (i=(m+3)/2;i<=n-(m-1)/2;i++)
 1652         respns[i]=0.0;
 1653      twofft(data,respns,fft,ans,n);
 1654      no2=n>>1;
 1655      for (i=2;i<=n+2;i+=2) {
 1656         if (isign == 1) {
 1657            ans[i-1]=(fft[i-1]*(dum=ans[i-1])-fft[i]*ans[i])/no2;
 1658            ans[i]=(fft[i]*dum+fft[i-1]*ans[i])/no2;
 1659         } else if (isign == -1) {
 1660            if ((mag2=ans[i-1]*ans[i-1]+ans[i]*ans[i]) == 0.0) {
 1661               log("Attempt of deconvolving at zero response in convlv().");
 1662               return(1);
 1663            }
 1664            ans[i-1]=(fft[i-1]*(dum=ans[i-1])+fft[i]*ans[i])/mag2/no2;
 1665            ans[i]=(fft[i]*dum-fft[i-1]*ans[i])/mag2/no2;
 1666         } else {
 1667            log("No meaning for isign in convlv().");
 1668            return(1);
 1669         }
 1670      }
 1671      ans[2]=ans[n+1];
 1672      realft(ans,n,-1);
 1673      free_dvector(fft,1,n<<1);
 1674      return(0);
 1675   }
 1676   
 1677   @ The |void sgcoeff(double c[], int np, int nl, int nr, int ld, int m)| routine
 1678   computes the coefficients |c[1...np]| for Savitzky--Golay filtering.
 1679   The coefficient vector |c[1...np]| is returned in {\sl wrap-around order}
 1680   consistent with the argument |respns| in the {\sl Numerical Recipes in C}
 1681   routine |convlv|.
 1682   ``Wrap-around order'' means that the first half of the array |respns|
 1683   contains the impulse response function at positive times, while the second
 1684   half of the array contains the impulse response function at negative times,
 1685   counting down from the highest element |respns[m]|.
 1686   The Savitzky--Golay filter coefficients are computed for |nl|
 1687   leftward (past) data points and |nr| rightward (future) data points, making
 1688   the total number of data points used in the window as $|np|=|nl|+|nr|+1$.
 1689   |ld| is the order of the derivative desired (for example, $|ld|=0$ for smoothed
 1690   function). Here $m$ is the order of the smoothing polynomial, also equal to the
 1691   highest conserved moment; usual values are $|m|=2$ or $|m|=4$.
 1692   The |sgcoeff| routine is adopted from {\sl Numerical Recipes in C}, 2nd Edn
 1693   (Cambridge University Press, New York, 1994).
 1694   
 1695   @<Routine for computation of coefficients for Savitzky--Golay filtering@>=
 1696   char sgcoeff(double c[], int np, int nl, int nr, int ld, int m)
 1697   {
 1698      void lubksb(double **a, int n, int *indx, double b[]);
 1699      void ludcmp(double **a, int n, int *indx, double *d);
 1700      int imj,ipj,j,k,kk,mm,*indx;
 1701      double d,fac,sum,**a,*b;
 1702   
 1703      if (np < nl+nr+1 || nl < 0 || nr < 0 || ld > m || nl+nr < m) {
 1704         log("Inconsistent arguments detected in routine sgcoeff.");
 1705         return(1);
 1706      }
 1707      indx=ivector(1,m+1);
 1708      a=dmatrix(1,m+1,1,m+1);
 1709      b=dvector(1,m+1);
 1710      for (ipj=0;ipj<=(m << 1);ipj++) {
 1711         sum=(ipj ? 0.0 : 1.0);
 1712         for (k=1;k<=nr;k++) sum += pow((double)k,(double)ipj);
 1713         for (k=1;k<=nl;k++) sum += pow((double)-k,(double)ipj);
 1714         mm=(ipj<2*m-ipj?ipj:2*m-ipj);
 1715         for (imj = -mm;imj<=mm;imj+=2) a[1+(ipj+imj)/2][1+(ipj-imj)/2]=sum;
 1716      }
 1717      ludcmp(a,m+1,indx,&d);
 1718      for (j=1;j<=m+1;j++) b[j]=0.0;
 1719      b[ld+1]=1.0;
 1720      lubksb(a,m+1,indx,b);
 1721      for (kk=1;kk<=np;kk++) c[kk]=0.0;
 1722      for (k = -nl;k<=nr;k++) {
 1723         sum=b[1];
 1724         fac=1.0;
 1725         for (mm=1;mm<=m;mm++) sum += b[mm+1]*(fac *= k);
 1726         kk=((np-k) % np)+1;
 1727         c[kk]=sum;
 1728      }
 1729      free_dvector(b,1,m+1);
 1730      free_dmatrix(a,1,m+1,1,m+1);
 1731      free_ivector(indx,1,m+1);
 1732      return(0);
 1733   }
 1734   
 1735   @ The |sgfilter(double yr[],double yf[],int mm, int nl, int nr, int ld, int m)|
 1736   routine provides the interface for the actual Savitzky--Golay filtering of
 1737   data. As input, this routine takes the following parameters:
 1738   \varitem[|yr[1...mm]|]%
 1739     {A vector containing the raw, unfiltered data}
 1740   \varitem[|mm|]%
 1741     {The number of data points in the input vector}
 1742   \varitem[|nl|]%
 1743     {The number of samples $n_L$ to use to the ``left'' of the basis
 1744     sample in the regression window (kernel). The total number of samples in
 1745     the window will be $nL+nR+1$.}
 1746   \varitem[|nr|]%
 1747     {The number of samples $n_R$ to use to the ``right'' of the basis
 1748     sample in the regression window (kernel). The total number of samples in
 1749     the window will be $nL+nR+1$.}
 1750   \varitem[|m|]%
 1751     {The order $m$ of the polynomial
 1752     $p(x)=a_0+a_1(x-x_k)+a_2(x-x_k)^2+\ldots+a_m (x-x_k)^m$ to use in the
 1753     regression analysis leading to the Savitzky--Golay coefficients.
 1754     Typical values are between $m=2$ and $m=6$. Beware of too high values,
 1755     which easily makes the regression too sensitive, with an oscillatory
 1756     result.}
 1757   \varitem[|ld|]%
 1758     {The order of the derivative to extract from the Savitzky--Golay
 1759     smoothing algorithm. For regular Savitzky--Golay smoothing of the input
 1760     data as such, use $l_D=0$. For the Savitzky--Golay smoothing and extraction
 1761     of derivatives, set $l_D$ to the order of the desired derivative and make
 1762     sure that you correctly interpret the scaling parameters as described in
 1763     {\sl Numerical Recipes in C}, 2nd Edn~(Cambridge University Press, New
 1764     York, 1994).}
 1765   \medskip
 1766   \noindent
 1767   On return, the Savitzky--Golay-filtered profile is contained in the vector
 1768   |yf[1...mm]|.
 1769   Notice the somewhat peculiar accessing of the coefficients $c_j$ via the
 1770   |c[(j>=0?j+1:nr+nl+2+j)]| construction in this routine. This reflects the
 1771   {\sl wrap-around storage} of the coefficients, where $c[1]=c_0$, $c[2]=c_1$,
 1772   $\ldots$, $c[|nr|+1]=c_{n_R}$, $c[|nr|+2]=c_{-n_L}$, $c[|nr|+3]=c_{-n_L+1}$,
 1773   $\ldots$, $c[|nr|+|nl|+1]=c_{-1}$.
 1774   
 1775   @<Routine for Savitzky--Golay filtering@>=
 1776   char sgfilter(double yr[], double yf[], int mm, int nl, int nr, int ld, int m) {
 1777      int np=nl+1+nr;
 1778      double *c;
 1779      char retval;
 1780   
 1781   #if CONVOLVE_WITH_NR_CONVLV /* Please do not use this $\ldots$ */
 1782      c=dvector(1,mm); /* Size required by the {\sl NR in C} |convlv| routine */
 1783      retval=sgcoeff(c,np,nl,nr,ld,m);
 1784      if (retval==0)
 1785         convlv(yr,mm,c,np,1,yf);
 1786      free_dvector(c,1,mm);
 1787   #else /* $\ldots$ use this instead. (Strongly recommended.) */
 1788      int j;
 1789      long int k;
 1790      c=dvector(1,nl+nr+1); /* Size required by direct convolution */
 1791      retval=sgcoeff(c,np,nl,nr,ld,m);
 1792      if (retval==0) {
 1793         for (k=1;k<=nl;k++) { /* The first |nl| samples */
 1794            for (yf[k]=0.0,j=-nl;j<=nr;j++) {
 1795               if (k+j>=1) {
 1796                  yf[k]+=c[(j>=0?j+1:nr+nl+2+j)]*yr[k+j];
 1797               }
 1798            }
 1799         }
 1800         for (k=nl+1;k<=mm-nr;k++) { /* Samples $|nl|+1 \ldots |mm|-|nr|$ */
 1801            for (yf[k]=0.0,j=-nl;j<=nr;j++) {
 1802               yf[k]+=c[(j>=0?j+1:nr+nl+2+j)]*yr[k+j];
 1803            }
 1804         }
 1805         for (k=mm-nr+1;k<=mm;k++) { /* The last |nr| samples */
 1806            for (yf[k]=0.0,j=-nl;j<=nr;j++) {
 1807               if (k+j<=mm) {
 1808                  yf[k]+=c[(j>=0?j+1:nr+nl+2+j)]*yr[k+j];
 1809               }
 1810            }
 1811         }
 1812      }
 1813      free_dvector(c,1,nr+nl+1);
 1814   #endif
 1815      return(retval); /* Returning |0| if successful filtering */
 1816   }
 1817   
 1818   @ Routines for removing preceding path of filenames.
 1819   In this block all routines related to removing preceding path strings go.
 1820   Not really fancy programming, and no contribution to any increase of numerical
 1821   efficiency or precision; just for the sake of keeping a tidy terminal output
 1822   of the program. The |strip_away_path()| routine is typically called when
 1823   initializing the program name string |progname| from the command line string
 1824   |argv[0]|, and is typically located in the blocks related to parsing of the
 1825   command line options.
 1826   The |strip_away_path()| routine takes a character string |filename| as
 1827   argument, and returns a pointer to the same string but without any
 1828   preceding path segments.
 1829   
 1830   @<Routine for removing preceding path of filenames@>=
 1831   @<Routine for checking for a valid path character@>@;
 1832   char *strip_away_path(char filename[]) {
 1833      int j,k=0;
 1834      while (pathcharacter(filename[k])) k++;
 1835      j=(--k); /* this is the uppermost index of the full path+file string */
 1836      while (isalnum((int)(filename[j]))) j--;
 1837      j++; /* this is the lowermost index of the stripped file name */
 1838      return(&filename[j]);
 1839   }
 1840   
 1841   @ In this program, valid path characters are any alphanumeric character or
 1842   `\.{.}', `\.{/}', `\.{\\}', `\.{\_}', `\.{-}', or `\.{+}'.
 1843   
 1844   @<Routine for checking for a valid path character@>=
 1845   short pathcharacter(int ch) {
 1846      return(isalnum(ch)||(ch=='.')||(ch=='/')||(ch=='\\')||(ch=='_')||
 1847          (ch=='-')||(ch=='+'));
 1848   }
 1849   
 1850   @ The |void showsomehelp(void)| displays a brief help message to standard
 1851   terminal output. This is where the caller always should end up in case
 1852   anything is wrong in the input parameters.
 1853   
 1854   @<Routine for displaying a brief help message on usage@>=
 1855   @<Routine for displaying a single line of the help message@>@;
 1856   void showsomehelp(void) {
 1857      hl("Usage: %s [options]",progname);
 1858      hl("Options:");
 1859      hl(" -h, --help");
 1860      hl("    Display this help message and exit clean.");
 1861      hl(" -i, --inputfile <str>");
 1862      hl("    Specifies the file name from which unfiltered data is to be read.");
 1863      hl("    The input file should describe the input as two columns contain-");
 1864      hl("    ing $x$- and $y$-coordinates of the samples.");
 1865      hl(" -o, --outputfile <str>");
 1866      hl("    Specifies the file name to which filtered data is to be written,");
 1867      hl("    again in a two-column format containing $x$- and $y$-coordinates");
 1868      hl("    of the filtered samples. If this option is omitted, the generated");
 1869      hl("    filtered data will instead be written to the console (terminal).");
 1870      hl(" -nl <nl>");
 1871      hl("    Specifies the number of samples nl to use to the 'left' of the");
 1872      hl("    basis sample in the regression window (kernel). The total number");
 1873      hl("    of samples in the window will be nL+nR+1.");
 1874      hl(" -nr <nr>");
 1875      hl("    Specifies the number of samples nr to use to the 'right' of the");
 1876      hl("    basis sample in the regression window (kernel). The total number");
 1877      hl("    of samples in the window will be nL+nR+1.");
 1878      hl(" -m <m>");
 1879      hl("    Specifies the order m of the polynomial to use in the regression");
 1880      hl("    analysis leading to the Savitzky-Golay coefficients. Typical");
 1881      hl("    values are between m=2 and m=6. Beware of too high values, which");
 1882      hl("    easily makes the regression too sensitive, with an oscillatory");
 1883      hl("    result.");
 1884      hl(" -ld <ld>");
 1885      hl("    Specifies the order of the derivative to extract from the ");
 1886      hl("    Savitzky--Golay smoothing algorithm. For regular Savitzky-Golay");
 1887      hl("    smoothing of the input data as such, use ld=0. For the Savitzky-");
 1888      hl("    Golay smoothing and extraction of derivatives, set ld to the");
 1889      hl("    order of the desired derivative and make sure that you correctly");
 1890      hl("    interpret the scaling parameters as described in 'Numerical");
 1891      hl("    Recipes in C', 2nd Edn (Cambridge University Press, New York,");
 1892      hl("    1994).");
 1893      hl(" -v, --verbose");
 1894      hl("    Toggle verbose mode. (Default: Off.)  This option should always");
 1895      hl("    be omitted whenever no  output file has been specified (that is");
 1896      hl("    to say, omit any --verbose or -v option whenever --outputfile or");
 1897      hl("    -o has been omitted), as the verbose logging otherwise will");
 1898      hl("    contaminate the filtered data stream written to the console");
 1899      hl("    (terminal).");
 1900   }
 1901   
 1902   @ In order to simplify the messaging, the |hl(const char *format, ...)| routine
 1903   acts as a simple front-end merely for compactifying the code by successive
 1904   calls to |hl(...)| rather than the full |fprintf(stderr, ...)|, still
 1905   maintaining all the functionality of string formatting in the regular
 1906   |printf()| or |fprintf()| syntax.
 1907   
 1908   @<Routine for displaying a single line of the help message@>=
 1909   void hl(const char *format, ...) {
 1910      va_list args;
 1911      char line[1024];
 1912   
 1913      va_start(args, format); /* Initialize args by the |va_start()| macro */
 1914      vsprintf(line, format, args);
 1915      va_end(args); /* Terminate the use of args by the |va_end()| macro */
 1916      sprintf(line+strlen(line), "\n"); /* Always append newline */
 1917      fprintf(stdout, "%s", line);
 1918      return;
 1919   }
 1920   
 1921   @ Routine for obtaining the number of coordinate pairs in a file. This
 1922   routine is called prior to loading the input data, in order to automatically
 1923   extract the size needed for allocating the memory for the storage.
 1924   
 1925   @<Routine for obtaining the number of coordinate pairs in a file@>=
 1926   long int num_coordinate_pairs(FILE *file) {
 1927      double tmp;
 1928      int tmpch;
 1929      long int mm=0;
 1930      fseek(file,0L,SEEK_SET); /* rewind file to beginning */
 1931      while((tmpch=getc(file))!=EOF) {
 1932         ungetc(tmpch,file);
 1933         fscanf(file,"%lf",&tmp); /* Read away the $x$ coordinate */
 1934         fscanf(file,"%lf",&tmp); /* Read away the $y$ coordinate */
 1935         mm++;
 1936         tmpch=getc(file); /* Read away any blanks or linefeeds */
 1937         while ((tmpch!=EOF)&&(!isdigit(tmpch))) tmpch=getc(file);
 1938         if (tmpch!=EOF) ungetc(tmpch,file);
 1939      }
 1940      fseek(file,0L,SEEK_SET); /* rewind file to beginning */
 1941      return(mm);
 1942   }
 1943   
 1944   @*Index.
 1945   

Return to previous page

Generated by ::viewsrc::

Last modified Saturday 17 Dec 2011