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
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
Data Loading...
