Contents of file 'sgfilter/sgfilter.tex':
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
Generated by ::viewsrc::