The cubic smoothing spline as a digital filter

LAMONT DOHERTY GEOLOGICAL OBSERVATORY

OF COLUMBIA UNIVERSITY

TREE-RING LABORATORY
Palisades, New York 10964

*

Technical Report No. CU-1-81/TR1

This material is based upon work supported by the
National Science Foundation under Grant No. ATM79-19223,
Gordon C. Jacoby, Principal Investigator.

*

THE CUBIC SMOOTHING SPLINE AS A DIGITAL FILTER

by

Kenneth Peters and Edward R. Cook
October 1981

*

Reproduction in whole or in part is permitted for any purpose of the
United States Government. Others must secure the author's
permission for use of this material.

★

Distribution of this document is unlimited.

LAMONT- DOHERTY GEOLOGICAL OBSERVATORY
of Columbia University
TREE-RING LABORATORY
Palisades, New York 10964

Technical Report No. CU-1-81/TR1

*

This material is based upon work supported by the
National Science Foundation under Grant No. ATM79-19223,
Gordon C. Jacoby, Principal Investigator.

*

THE CUBIC SMOOTHING SPLINE AS A DIGITAL FILTER

by

Kenneth Peters and Edward R. Cook
October 1981

*

Reproduction in whole or in part is permitted for any purpose of
the United States Government. Others must secure the author's
permission for use of this material.

*

Distribution of this document is unlimited.

■k

Any opinions, findings, and conclusions or recommendations expressed
in this publication are those of the authors and do not necessarily
reflect the views of the National Science Foundation.

Digitized by the Internet Archive
in 2020 with funding from
Columbia University Libraries

https://archive.org/details/cubicsmoothingspOOpete

THE CUBIC SMOOTHING SPLINE AS A DIGITAL FILTER

Kenneth Peters and Edward R. Cook

ABSTRACT

Reinsch (1967) introduced the cubic smoothing spline as the solu¬
tion to a certain variational problem in the theory of data smoothing.

This smoothing operator has good transient response, a minimal curva¬
ture property, and compensation for end effects so that all the data
is used. In addition, the data can have arbitrary spacing and weighting.

In this paper we represent Reinsch1 s equations in terms of convolu¬
tions for the special case of equally-spaced and equally-weighted data.

We then use a Fourier transform to obtain a frequency domain represen¬
tation, solve the characteristic polynomial of the system, and apply an
inverse transform to the frequency response function to obtain the
impulse response function. The Lagrangian multiplier of the variational
problem parameterizes the family of digital filters derived in this man¬
ner. This representation which is exact for finite data sets extends the
range of application of the smoothing spline and, incidentally, simplifies
the algebra and reduces the amount of computer time and storage required,
by allowing one to select the degree of smoothing on the basis of frequency
domain considerations. One can do time domain filtering without end
effects in a computationally efficient manner.

Smoothings of four tree-ring width sequences from forest interior
sites are presented to illustrate the particular suitability of these
filters for data exhibiting both non-stationary and episodic behavior.

1

THE CUBIC SMOOTHING SPLINE AS A DIGITAL FILTER

I INTRODUCTION

Reinsch [7] discusses the problem of smoothing noisy data when the
underlying function is not known. We begin with a brief statement of
his main results in his own notation and terminology. He suggests that
the fitting function, g(x), minimize the total square curvature,

2

dx

(la)

xonder the constraint that a weighted sum of squares of the residuals
be bounded above ,

n

"g(Xi)-yil 2 < S;

6Yi

(lb)

where the y^ represent the noisy data, the 6y^ are a set of weights
and S is a scaling parameter. Equivalently, introducing a Lagrangian
multiplier, p, and an auxiliary variable, z, the functional

is to be minimized.

The solution is a cubic spline, i.e., a concatenation of cubic
segments :

g^(x) = ai + (x-x^) + c^ (x-x^) 2 + d^ (x-x^) 3, x^ ^ x < x^+1, (2)

such that g(x) , dg(x)/dx, and d2g(x)/dx- are continuous at x^ for
i=l, n-1.

2

If the Lagrangian multiplier, p, is known, the values can be
computed from ^

(QD Q + pT ) c = pQy, (3)

where the circumflex denotes transposition, and the a^ from

Then

a = y -

p-1D2Qc.

(4)

di = (ci+1 - ci) / 3hi, i = 0, n - 1 (5)

where, for the sake of uniformity of notation, Reinsch chooses to set
Cq and c^ equal to 0 ; and

bi = E(ai+1 ~ a±) /hi3 - cihi - dihi , i = 0, n - 1.
The notation used is Reinsch' s:

hi = *

i+1 " xi;

c (c]_, • • • 'cn-l^ '

(6)

(7a)
(7b)
(7c)
( 7d)
(7e)
(7f)

y = (y0# yi . yn.r yn) *

a ^ a0 ' al' • • • '^-1' an^ '

D = diag (6yQ, . . . ,6yn) ;

T = a positive-definite tridiagonal matrix of order n-1,

where t^ = 2 (hi_1 + hi)/3, and ti>i+1 = h±/3 = ti+1/i;

Q = a tridiagonal matrix with n+1 rows and n-1 columns (7g)

where qi_1^i = l/hi_1, qi;j_ = - (l/hj^ + l/hi) , and

gi+l,i " 1/bi*

For Reinsch1 s intended application to data smoothing it is natural
to specify the weights 5y^, usually as the square root of the error
variance at x=x^ if it is known or can be estimated, and S, as the value

of a chi-square variable with n degrees of freedom. He then proceeds to

3

and may be of interest to the reader. Reinsch t 8] recasts the problem of
reference £71 in much more general terms and shows that similar results hold.
Demmler and Reinsch £41 achieve by advanced methods a very general result of
which formula (16) of this paper is a special case. Wahba £9] discusses the
problem of optimizing the degree of smoothing given certain assumptions
regarding the random component of the data and gives explicit formulas
which can be used in conjunction with the heuristic considerations which
motivated this paper. Wahba and Wold £10] and Craven and Wahba [3] apply a
method of cross-validation to select the smoothing parameter and in £11]
Wahba and Wold discuss the use of periodic splines in spectral density
estimation which is related to a suggestion for natural splines in this
paper in Section V. Wold £12] is concerned with the suitability of splines
of various kinds from a data analytic point of view. In this paper we are
concerned with a discrete problem — the relationship between the vectors y
and a. In the limit of large p there are no digital filtering effects.
Horowitz C6] considers the continuous problem of the filtering effects of
common interpolating splines. His analysis may be of interest if the dis¬
crete data arise from accumulation or from sampling continuous records.

The references given here should provide an adequate introduction to
smoothing splines.

The rest of the paper consists of four sections and two appendices.
Section II contains the derivation of the convolutional representation
and of the frequency response function. Section III contains the solution
of the characteristic equation of the system. Section IV contains the deri¬
vation of the impulse response function. Section V includes a discussion
of how to interpret and apply the representation derived in sections II,

III and IV, and an example which illustrates the time domain behavior of

5

the cubic smoothing spline. Appendix A is a discussion of this behavior in

the limit of small curvature, and Appendix B contains an analysis of the
coupled ordinary differential equations of which the cubic smoothing spline
is an interpolated finite difference solution.

The reader interested only in practical applications should read
sections I and V and consider the main practical results which are in
equations (16) and (33) and figures 1 and 5. The critical reader will
want to include the beginning of section II. The appendices and the
parts of the text relating to the correction terms (indicated by primes)
and the small curvature limit, can be ignored, but clues to understanding
the Fourier representation are contained therein. Finally, those filters
generated by the real negative roots of the characteristic polynomial
are of little practical interest since they result in very little
smoothing.

II THE CONVOLUTION REPRESENTATION AND THE FREQUENCY RESPONSE FUNCTIONS

Consider the special case of equations (3) and (4) for which
h^ = 1.0 = Sy^ for all i. In this case the matrices Q and T and R =

QQ + pT will be Toeplitz, i.e. , all the terms on the same diagonal will
be equal, and symmetric. The matrix- vector product Rc will be equivalent
to convolution with a symmetric (phase-free) filter except for the first
and the last two points. In fact these restrictions are necessary
because convolution is defined as a time-invariant linear process, and
it is impractical to do Fourier transforms on unequally spaced data.

The particular numerical values of the spacing, in this case unity,

-3

serves to scale or normalize the parameter p which has the units of x
With these restrictions then the desired representation can be achieved
as follows.

6

First, extend a typical row of the matrix R by appending an infinite
number of zeroes on each end, and do the same for Q and Q. This exten¬
sion is necessary because the limits on the Fourier integral are infinite.
The particular form of the extension (zeroes) is dictated by the band
form of the corresponding matrices. Thus define

r = {...,0,1, (-4+p/3) , (6+4p/3) , (-4+p/3) ,1,0, . . . } (8)

and

q = {... ,0,1, -2, 1,0,.. .} (9)

in accordance with (7f) and (7g) . The subscript p makes the para¬
meterization explicit. The vectors c, y and a can be extended in a
similar manner, i.e.,

c = {.. . ,0,cQ,c1,... ,cn_1,cn,0, .. .}, c0=cn=0, (7b')

y = {...,0,y0,y1,...,yn_1,yn,0,...}, (7c1)

and

a= { . . . ,0,aQ,a1, . . . ,an_1,an,0, . . . }; (/a')

but this particular extension for c, y and a is not dictated by

Reinsch's equations (except for the end conditions Cg=cn=0) which

entail no assumptions regarding the behavior of y outside the data

interval. We use it only for definiteness and convenience. More

general extensions will be discussed in section V. The main point

here is to represent Reinsch's computations as exactly as possible.

Next notice that while (3) equates two vectors of length n-1,

r * c, where * denotes convolution, and q * y will in general have n + 3
J?

non-zero elements. To represent (3) accurately the extra terms gener¬
ated by the convolutions must be subtracted from both sides before they
are set equal. Thus, (3) becomes

rp*c - = p (q * y-y') ( 3 ’ )

7

with

c ' = { . . . ,0,0^ , i (-4+p/3) c^+c^ ,0 , . . . ,0, Ccn_2+ (-4+p/3) cn_12 ,
cn_]_,0, . . . } (10)

and

y' = { • • • /0,y0, (-2y0+y1) ,0, . . . ,0, (-2yn_1+yn) ,yn,0, . . . }. (11)

Note that the prime will not be used to denote differentiation in this
paper.

In the case of equation (4) the extra terms are retained. It
becomes

y - — q
P

c, p > 0,

(4')

Finally, to achieve symmetry in the transforms, and for representa¬
tional purposes, the discrete vectors will be represented as weighted
sums of delta functions. Letting a denote any of them define

a(t) 7 a. 6(t-i). (12)

*• — l

i=-°°

The explicit dependence on t will identify continuous time domain
functions.

Fourier transforming (3') and (4') modified according to (12)

gives

r (f)c(f) - Cp ( f ) = p Cq(f)y(f) - y'(f)]

(13)

and

a (f ) = y(f) - |q(f)c(f) , p > 0.

ir

(14)

Solving

(13)

for c(f), substituting in (14) and rearranging gives

a (f)

= u (f)y(f) + u ' (f) Cy' (f) - ic • (f)] , p > 0,
p P P P

(15)

where

u (f)

= (p/6) cos2TTf+ (p/3)

(16)

P

cos^2'n'f+ ( (p/6) -2) cos2nf+ ( (p/3) +1)

and

n , _ 2 (cos2TTf-l)

(17)

cos^2iTf+ ( (p/6) -2) cos27rf+ ( (p/3) +1

which are the

frequency response functions.

The

functions u (f) and u '(f) are shown in figures 1 'and 2 over
P P

the

useful

range

of p.

8

Ill THE ROOTS OF THE CHARACTERISTIC POLYNOMIAL

The characteristic polynomial of the system is the (common) denomi¬
nator of u (f) or u ' (f) expressed in terms of z = exp(27rif) . The roots
I? P

of this function will be needed in section IV. Setting the denominator
of Up(f) equal to zero and solving for cos2irf gives

and

Zp = 1- (p/12) + i(p(72-p))1/2/12

z" = 1- (p/12) - i (p (72-p) ) 1/2/12

(18)

(19)

The principle root is used here and always. These roots lie on the

circle of radius 3 centered on the real axis at -2 for 0<p$72.

For p£72 Im z =0. Also, to avoid ambiguity, lim z as p-*30 is defined
P P

to be -2.

If z=exp (27rif ) is a root, so is 1/z since cos 2nf = (exp(2irif) +
exp (-27rif ) ) /2 . Using this relationship and its converse, the four roots
are, for 0<p<72,

and

++ + • +2» 1/2

z = z + i ( 1- z )

P P P

H — + . , +2 . 1/ 2

z = z - i(l-z )

p p p

-+ - -2 1/2
Zp = zp + i(l-zp )

-2 1/2

z = z - i(lz )

P P P

(20)

(21)

(22)

(23)

For p^72 there are four distinct negative real roots. Two of

these, z++ and z +, remain within the unit circle,

P P

and

z +
P

(24)

(25)

9

while the reciprocals, z^+ and z ^ , respectively, remain outside. The

lim z as p-*°° is defined to be 0, while lim z as p-*» is -2 + /3.

P P

The root z^"*" is plotted in figure 3. The other roots are simply
related to z^+. Thus z^ is the reciprocal, z^~ is the complex con¬
jugate, and z is the reciprocal of the complex conjugate.

P

IV THE IMPULSE RESPONSE FUNCTIONS

Using the results of section III, the impulse response functions

can be obtained as the inverse Fourier transforms of the frequency

response functions u (f) and u ' (f) . Decomposing the right side of (16)

P P

into partial fractions gives

u (f) =

Ar

Br

cos2,n'f - z2

cos2iTf - zT

(26)

where

A =

(p/6) Zp + (p/3)

z+ - Z~

and Br

(p/6) Zp + (p/3)

z~ — z+
p p

(27)

Then

u (t) = A ! exp(2TTift) df + B j exP (2irif t)

COs2fTf - z"

P t cos2irf - zx

d p

df

(28)

P

Tri

z^dz

z2-2z+ + 1

P

+ !p

ttT

z^dz

oo

E

5 (t-i)

^ z -2z z + 1

<* p

C *

(29)

the unit circle, C, being traversed in the positive sense for t>0. To
see that (29) follows from (28) consider (12) and the fact that the
series of delta functions used to repeat the principal part of u^(f)
by convolution is the transform of the series in t.

10

The integrals in (29) can be evaluated by the method of residues Cl] .

The roots of the denominator of the first integrand are z++ and z+~ ,

P P

and for the second z and z , which are respectively complex conjugates

P P

of the first (reciprocal) pair. Then, for all p>0, the impulse response
function

u (t)

P

(2K z++|t|+ 2
P P

V

-+ t

CO

) .XI <5 (t-i) =

l=-oo

V (t)

P

.£

l=_oc

6 (t-i) ,

where

K =

++

A / (z
P P

++

Vzp ) ,

and

L„ =

B /(z
P P

-+

Vzp ) .

(30)

(31)

(32)

Using the same conjugacy relationships between the roots, for
0<p<72 ,

Vp(t) = 4Re{Kpz++ I t I } = 4|Kp| |Zp+| I tlcos (Arg Kp + |t|Arg Zp+) (33)

In the limit as p-K), vp(t) -> 0 uniformly everywhere.

For p>72 there are four distinct negative roots occurring in reciprocal
pairs. The impulse response function becomes the sum of two complex numbers
with different amplitudes (decaying exponentially in |t|) rotating in
opposite senses around the t axis, piercing the imaginary-t plane at half
integral values of t and the real-t plane at integral values of t. In the
limit as p-*00, v (t) 1 for t=0 and -*■ 0 elsewhere. That is, u (t) is a

ir ir

spike at the origin.

The function v^ (t) must be evaluated by a limit process. The limits
from the right and the left are the same and thus Vp(t) is continuous at
p=72 for all values of t. The function v (0) in particular is plotted in
figure 4 .

11

The functions v^(t) for log p = -3, -4, -5

The analysis for u ' proceeds as for u

P p

and B ' are different.

P

are plotted in figure 5 .
except that the constants

A ' = (2z+ - 2) / (z+ - z")
P P P P

and

(34)

-3

V = (2zp ' 2,/(zp " ZP (35)

v ' (0) is plotted in figure 4, and the functions v ' (t) for log p =
P P

-4,-5 are plotted in figure 6.

V DISCUSSION AND AN EXAMPLE

Transforming (15) back into the time domain and reverting to
the discrete vector symbolism gives

a = u*y-u'* (y ' - — c' ) , p > 0, (36)

1 p 7 p p ^

and the operation of the smoothing spline is seen to be representable
in the time domain as the sum of two convolutions. The first term
accurately represents the operation when end effects can be ignored;
for instance far enough away from the ends of the data so that the
filter cannot feel them, or in the cases in which a reasonable extra¬
polation of the data can be made. The second term is peculiar to the
smoothing spline and accounts for its apparent lack of end effects.

It is necessary to consider the sense in which the smoothing spline has
no end effects.

This phrase is appropriate if one accepts the global variational
definition of the smoothing spline as implying that its quasi-local
behavior is the same at the ends of the data set as near the center.

12

Indeed the Euler-Lagrange equations must be satisfied at all interior

points. In this sense the second term in (36) "fixes" the computations

so that u alone is an adequate description of its operation, in which
P

case one need only consult figure 1 in selecting an appropriate p value.
We have been using the spline in this sense. It is more appropriate in
this paper, however, to consider the end behavior in terms of the con¬
volution representation.

If a different choice had been made in extending the vectors c, y
and a, the contribution of the second term in (36) would have been dif¬
ferent. In fact there are (an infinite number of) extensions such that
this second term would contribute nothing. They can be computed by con¬
volving Up with y and requiring that the resulting values be equal to a^
for 0<i<n. If one accepts these extensions as "reasonable", then again
u^ by itself adequately describes the filter. This implies that the
smoothing spline can be used as a predictor or to achieve higher reso¬
lution in spectral estimates. If the physical model entails that the
data vary according to some principle of minimum curvature at least at
lower frequencies such an idea is certainly tenable, but as in other
kinds of prediction and spectral enhancement (maximum entropy, for
instance) due consideration must be given to the underlying model.
Prediction is reliable only if the model is appropriate.

Four smoothings of tree-ring width sequences are shown in figure 7.
The same figure appears in [2] . The dashed lines show orthogonal poly¬
nomial fits generated by the INDXA program C53 the orders of which are
determined by stopping when an additional order fails to reduce the
residual variance by 5% or more. The log p value was -4 in each case.

13

Tree-ring width sequences from forest interior sites tend to exhibit
behavior similar to that seen in records of complex physical processes
such as level changes in streams and reservoirs and climatic variations.
These data exhibit sudden changes in level (episodic behavior) over a
wide range of time scales. The spectra tend to be lumpy and have concen¬
trations of power near zero frequency. The particular suitability of the
smoothing spline as a low pass filter for certain tree-ring records, and
perhaps for other similar data, is accounted for by the shape of the
impulse response functions. The spikey exponential envelope with the |t|
dependence ensures that the filter will be responsive to sudden changes
in level. On the other hand, the fact that it is a low-pass filter with
a minimum curvature property ensures that power near zero frequency will
be preserved. The frequency response function is smooth at the origin.
Thus trends and slow wandering in the data will be preserved. Selecting
a value of S (or p) affects a compromise between fitting and smoothing
the data. We see how this compromise manifests itself in the digital
filters .

Programming of equations (31) and (4') can be done in various ways.

We have not tried to optimize the computations beyond using an equation
solver specialized for symmetric band matrices. The user concerned about
storage might try a Toeplitz-specif ic algorithm. The only caution regards
the possible necessity of using double precision accumulation variables.

We found that in single precision (32 binary bits) the computations were

-4

unstable for values of p smaller than 10 roughly.

14

-

APPENDIX A

It is not obvious that the weight function shown in figure 6 will,
in the limit of small p, give a series of values which lie on a straight
line for O^i^n and be zero outside this interval when convolved with
anything at all. The purpose of this appendix is to examine this limit¬
ing case.

For 0<p<72, and thus in the limit as p -*■ 0 ,

II 1 1 1

v ' (t) = 4Re{Kp ' z^+ ' 1 }=4 | Kp ' | | z++ | cos (Arg Kp+|t|Arg z^+) . (Al)

Letting r, 0' and 0 denote small positive numbers, for small p Al can
be written

v ' (t) = 4 Ik ' | (l-r) 1 1 1 cos ( (3tt/4) -0' - 1 1 1 0 ) . (A2)

p P

Expanding in r and 0 ' and 0 ,

v 1 (t) = k [1 - 1 1 1 r + o(r2)] El-0' - 1 1 1 0 - o(02)], (A3)

P

where

k = -4|K ' |/V2.

To first order therefore

v ' (t) = k(l-0‘- s 1 1 1 ) , for s < < 1 1 1 < < 1/s, (A4)

ir

where s = r + 0.

1

Denoting the four non-zero values of (y'-vc ') by w ,

ir P “ 1

wn+l' at least the following four requirements must be met:

k(l-0,)w_1 + k (1-0 ' -s) wQ = 0,

k(l-@'-s)w_1 + k (1-0 ' ) wQ = aQ,

k(l-0')wn + k(l-0'-s)wn+1 = an,

and k(l-0'-s)w + k(l-0')w = 0.

n n+1

w

O'

w

n

(A5)

(A6)

(A7)

(A8)

15

Since in the limit of small p,|K 1 | and thus |k| become very large,

ir

k can be divided out of A5 - A8. Solving the resulting equations for

w and w (A5 and A6) and w and w (A7 and A8) yields to first order
-1 0 n n+1

in s ,

w_x = -aQ (1-0 ' -s) /2s ( 1-0 ' ) ,

(A9)

wQ = a0/2s.

(A10)

w - a /2s,
n nx '

(All)

and

wn+l = -an(l-0,-s)/2s(l-0') .

(A12)

Since s > 0 for all finite p these solutions are defined for the limit.
Also, by the symmetric linear (to first order) dependence of v ' (t)

IT

on t, it follows from A5 that a^ = 0 for all -1/s <<i< -1; from A8,
that a^ = 0 for all n+l<i<<l/s; and from A6 and A7 that a^ varies
linearly with i (to first order) for O^i^n. In the limit as p-*0 the
first order expressions become exact and l/s-*-°°.

16

APPENDIX B

The first terms in the central finite difference expansions of
the second and fourth derivative operators are obvious in the matrices
Q and R. Corresponding to equations (3') and (4') there are the coupled
differential equations

CD4 + (p/3) (D2 + 6)] C(t) = pD2Y(t) (Bl)

and

A ( t) = y - — D2C ( t) , p>0 . (B2)

P

Conversely, a finite difference solution of (Bl) and (B2) will
satisfy equations (3') and (4'). More generally, the theory of Reinsch's
smoothing splines is related to the theory of the interpolated finite
difference solutions of certain ordinary differential equations with
constant coefficients. The remainder of this appendix is a brief
analysis of the equations (Bl) and (B2) , analogous to that done for
the spline equations, which will facilitate some comparisons the reader
may wish to make between the continuous and discrete systems.

Fourier transforming equations (Bl) and (B2) and solving for A(f)
in terms of Y(f) , gives

p(3-27T2f2)

A (f ) = - — - — — Y(f) =U (f)Y(f). (B3)

24irf +p(3-2tt f2)

Note that this same transform would arise from the single equation

(D4 + (p/3) (D2 + 6))A(t) = (p/3) (P2 + 6) Y ( t) .

(B4)

This won't work in the discrete case because of the correction terms.
The roots of the denominator for 0$p-$72 are

fp+ = (P + i (p (72-p) ) 1/2 ) /24tt2

cin d i ^ p

f2_ = (p - i (p (72-p) ) )/24tt

ir

(B5)

(B6)

17

and therefore

and

+<£p+)1/2-

-if2*) 1/2

P

-(f2_)1/2.

P

(B7)

(B8)

(B9)

(BIO)

Corresponding to the constants A and B

P P

there are

2 2+ 4 9+ o_

a = p (3-2tt f ) /24tt ( r + - f ),
P P P P

and

2 2- 4 2- 2+

3 = p(3-2rr f ) /24ir (f - f ) ,

P P ir P

in the inverse transform

(Bll)

(B12)

while

k

P

Vf

++
P '

(B13)

and

£P - vf;+

(B14)

The contour of integration consists of the real axis and the upper
half of the point at infinity, traversed counterclockwise for t>0.

For 0<p<72

U (t) = -2iTexp (-2iTlm{ f^+} 1 1 1 ) Im{k^exp (27riRe{fp+} 1 1 1 ) } . (B15)

For p^72, each contour (t<0) includes two of the real roots and

U (t) = -ir(k sin 2tt f 1 1 1 + 1 sin 2ir f_+ 1 1 1 ) . (B16)

P P P P P

There is no decay in the continuous case for p^72.

Evaluation of the limits as p approaches 72 from the left and right

will show that for any value of t, U (t) is continuous at p=72.

P

18

ACKNOWLEDGEMENTS

The authors wish to thank Paul Stoffa, John Shepherd and Keith
McCamy for many constructive criticisms. Funds for this work were pro¬
vided through the Climate Dynamics Program of the National Science Founda¬
tion, grant no. ATM-77-19217.

19

FIGURE CAPTIONS

Figure 1.

Figure 2.

Figure 3.

Figure 4.

Figure 5.

Figure 6 .

Figure 7.

The gain of the frequency response function u (f) for various

P

values of the parameter p over the range .001^f^.5.

The base ten logarithm of the negative of the frequency response
function u ' (f) for various values of the parameter p over the

.c'

range . 001^f<.5. The dashed line passes through the maximi
(circles) on these curves.

The root z++ of the characteristic polynomial r (z) . Values of
P P

the parameter p are shown along the curve .

The central values vp(0) and v^ 1 (0) of the smooth parts of the
impulse response functions u^ and u ' .

The smooth part v^(t) of the impulse response function u^ for
three values of p most useful for standardizing forest interior
tree-ring series.

The smooth part v ' (t) of the impulse response function u^ '
for three values of p most useful for standardizing forest
interior tree-ring series.

Four examples of the smoothing spline applied to forest interior
ring-width series. Each spline was computed for log p = -4.0.
The solud curve is the spline fit and the dashed curve is the
orthogonal polynomial fit from program INDXA (see ref. 4) . The
four series are from the same site in southeastern New York and
the lower two plots are curves from the same tree.

20

1.0

0.9

0.8

0.7
0.6

0.5'
0.4
0.3

02
0.1

0.0

-3.0 -2.0 -1.0 0.0

Q.

LOG f

Figure 1

2.5

LOG f

Figure 2

LOG [- up '( f )]

++

dZ

ro

Figure

LOG

-25

RINGWIDTHS (mm) RINGWIDTHS (mm)

1.5-

A

E

S

cn 1-0-

X

I-

9

5

o

z

* 0.5-

0.0

r— I— t— — t— i - 1 - r— t— — i - r - r - v. . , - r

1850 1900 1950

I 1 "—'-I 1 ■ " I ■ ' I "I1 ■ — T ' I - 1 I I . n - 1 ' 1 — 1111

1800 1850 1900 1950

— , - i 1 i - 1 i - 1 - 1 - r ■ 11 - 1 - 1 - 1 - 1 - 1 - 1 - - - 1 - 1 - 1 - -

1750 1800 1850 1900 1950

YEARS

Figure 7

REFERENCES

1. Churchill, R.V. : Introduction to complex variables and applications.

New York: McGraw-Hill Book Company, Inc., 1948.

2. Cook, E.R. and Peters, K.: The smoothing spline: a new approach to

standardizing forest interior ring-width series for dendroclimatic
studies. (submitted) .

3. Craven, P.G.: Smoothing noisy data with spline functions, Numer. Math.

31, 377-403 (1979) .

4. Demmler, A., Reinsch, C.: Oscillation matrices with spline smoothing,

Numer. Math. 24, 375-382 (1975).

5. Fritts, H.C., Mosimann, J.E. and Bottorff, C.P.: A revised computer

program for standardizing tree-ring series. Tree-Ring Bulletin
29, 15-20 (1969).

6. Horowitz, L.L. : The effects of spline interpolation on power spectral

density. IEEE Transactions on Acoustics, Speech, and Signal Pro¬
cessing, vol. AAP-22 , number 1, 22-26 (February 1974) .

7. Reinsch, C.H. : Smoothing by spline functions, Numer. Math. 10, 177-183

(1967) .

8. Reinsch, C.H„: Smoothing by spline functions, II. Numer. Math. 16,

451-454 (1971) .

9. Wahba, G. : Smoothing noisy data with spline functions, Numer. Math. 24,

383-393 (1975) .

10. Wahba, G. , Wold, S.: A completely automatic French curve: Fitting

spline functions by cross-validation, Comm. Stat. 4, 1-17 (1975).

11. Wahba, G. , Wold, S.: Periodic splines for spectral density estimation:

The use of cross-validation for determining the degree of smoothing.
Comm. Stat. 4, 125-141 (1975).

12. Wold, S.: Spline functions in data analysis, Technometrics 16, 1-11 (1974) .

21

LAMONT DOHERTY GEOLOGICAL OBSERVATORY

OF COLUMBIA UNIVERSITY