Return to previous page

Contents of file 'sgfilter/sgfilter.tex':




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

Return to previous page

Generated by ::viewsrc::

Last modified Saturday 17 Dec 2011