Fast Compact Algorithms and Software for Spline Smoothing

Fast Compact Algorithms and Software for Spline Smoothing investigates algorithmic alternatives for computing cubic smoothing splines when the amount of smoothing is determined automatically by minimizing the generalized cross-validation score. These algo

  • PDF / 200,424 Bytes
  • 7 Pages / 439.37 x 666.142 pts Page_size
  • 60 Downloads / 200 Views

DOWNLOAD

REPORT


FFT Algorithm

For a given l, the cubic spline smoother is a time-varying linear filter, but it can be approximated by a time-invariant linear filter. It will then be amenable to frequency domain analysis and implementation using the FFT. The FFT algorithm will be compared to the Cholesky algorithm in terms of execution time, accuracy, and memory use. For digital signal processing background, see [1]. Other results on splines in the frequency domain can be found in [2].

4.1

Frequency Response of the Spline Smoother

Eq. (1.7) can be viewed as a fourth-order difference equation ci þ ðl=6  4Þciþ1 þ ð2l=3 þ 6Þciþ2 þ ðl=6  4Þciþ3 þ ciþ4 ¼ yiþ2  2yiþ3 þ yiþ4 ;

(4.1)

for 1  i  n  6; with left boundary conditions ð2l=3 þ 6Þc1 þ ðl=6  4Þc2 þ c3 ¼ y1  2y2 þ y3 ; ðl=6  4Þc1 þ ð2l=3 þ 6Þc2 þ ðl=6  4Þc3 þ c4 ¼ y2  2y3 þ y4 ;

(4.2)

and right boundary conditions cn5 þ ðl=6  4Þcn4 þ ð2l=3 þ 6Þcn3 þ ðl=6  4Þcn2 ¼ yn3  2yn2 þ yn1 ; cn4 þ ðl=6  4Þcn3 þ ð2l=3 þ 6Þcn2 ¼ yn2  2yn1 þ yn :

(4.3)

Suppose we ignore the boundary conditions and consider (4.1) to be valid for all i. We will also ignore the first two and last two equations in (1.8) and use

H.L. Weinert, Fast Compact Algorithms and Software for Spline Smoothing, SpringerBriefs in Computer Science, DOI 10.1007/978-1-4614-5496-0_4, # The Author(s) 2013

29

30

4 FFT Algorithm

si ¼ yi  ci2 þ 2ci1  ci ;

(4.4)

for all i. Taking bilateral z-transforms of (4.1) and (4.4), 

 z4 þ ðl=6  4Þz3 þ ð2l=3 þ 6Þz2 þ ðl=6  4Þz þ 1 CðzÞ ¼ z2 ðz  1Þ2 YðzÞ; SðzÞ ¼ YðzÞ  z2 ðz  1Þ2 CðzÞ:

(4.5)

The transfer function HðzÞ of the spline smoother is thus HðzÞ ¼

SðzÞ lz2 ðz þ 4 þ z1 Þ ¼ : YðzÞ 6ðz  1Þ4 þ lz2 ðz þ 4 þ z1 Þ

(4.6)

The denominator polynomial is a scaled version of (2.15). Since there are no poles (denominator roots) on the unit circle, a frequency response HðoÞ exists and can be obtained by replacing z in (4.6) with ejo . After some simplification, HðoÞ ¼

lð2 þ cos oÞ 12ð1  cos oÞ2 þ lð2 þ cos oÞ

:

(4.7)

This frequency response is real and periodic with period 2p. It is also even so we can restrict attention to the interval o 2 ½0; p. Fig. 4.1 shows plots of HðoÞ for various l. We see that HðoÞ is nonnegative, achieves a maximum value of one at the origin, and decreases monotonically with increasing frequency (no ripples). Furthermore, lim HðoÞ ¼ 1; 8o; ( 1; o ¼ 0 ; lim HðoÞ ¼ l!0 0; o 6¼ 0 l!1

HðpÞ ¼

l : l þ 48

(4.8)

(4.9)

For small o, H ðoÞ ffi 1  l1 o4 :

(4.10)

which implies that the first three derivatives of HðoÞ are zero at the origin. Consequently, the smoother passes cubic polynomials unchanged. The 3 dB cutoff frequency (or bandwidth) in pi units, plotted in Fig. 4.2, is   pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi o0 1 ¼ p1 cos1 1 þ 24 l  lðl þ 144Þ ;  ¼ 2  1; p

(4.11)

4.1 Frequency Response of the Spline Smoother

Fig. 4.1 Frequency response

Fig. 4.2 Bandwidth

31

32

4 FFT Algorithm

pffiffiffi  which is well-defined as long as l  48 2 þ 1 , in which case the smoother is a (zero-phase) lowpass filte