@@ -5,14 +5,14 @@ Cubic{BC<:Flag}(::BC) = Cubic{BC}()
55Assuming uniform knots with spacing 1, the `i`th piece of cubic spline
66implemented here is defined as follows.
77
8- y_i(x) = ϕm p(x-i) + ϕ q(x-i) + ϕp q(1- (x-i)) + ϕpp p(1 - (x-i))
8+ y_i(x) = cm p(x-i) + c q(x-i) + cp q(1- (x-i)) + cpp p(1 - (x-i))
99
1010where
1111
12- p(☆ ) = 1/6 * (1-☆ )^3
13- q(☆ ) = 2/3 - ☆ ^2 + 1/2 ☆ ^3
12+ p(δx ) = 1/6 * (1-δx )^3
13+ q(δx ) = 2/3 - δx ^2 + 1/2 δx ^3
1414
15- and the values `ϕX ` for `X ⋹ {m, _, p, pp}` are the pre-filtered coefficients.
15+ and the values `cX ` for `X ∈ {m, _, p, pp}` are the pre-filtered coefficients.
1616
1717For future reference, this expands out to the following polynomial:
1818
@@ -24,6 +24,11 @@ When we derive boundary conditions we will use derivatives `y_0'(x)` and
2424"""
2525Cubic
2626
27+ """
28+ `define_indices_d` for a cubic b-spline calculates `ix_d = floor(x_d)` and
29+ `fx_d = x_d - ix_d` (corresponding to `i` `and `δx` in the docstring for
30+ `Cubic`), as well as auxiliary quantities `ixm_d`, `ixp_d` and `ixpp_d`
31+ """
2732function define_indices_d {BC} (:: Type{BSpline{Cubic{BC}}} , d, pad)
2833 symix, symixm, symixp = symbol (" ix_" ,d), symbol (" ixm_" ,d), symbol (" ixp_" ,d)
2934 symixpp, symx, symfx = symbol (" ixpp_" ,d), symbol (" x_" ,d), symbol (" fx_" ,d)
@@ -39,6 +44,14 @@ function define_indices_d{BC}(::Type{BSpline{Cubic{BC}}}, d, pad)
3944 end
4045end
4146
47+ """
48+ `define_indices_d` for a cubic, periodic b-spline calculates `ix_d = floor(x_d)`
49+ and `fx_d = x_d - ix_d` (corresponding to `i` and `δx` in the docstring entry
50+ for `Cubic`), as well as auxiliary quantities `ixm_d`, `ixp_d` and `ixpp_d`.
51+
52+ If any `ixX_d` for `x ∈ {m, p, pp}` (note: not `c_d`) should fall outside of
53+ the data interval, they wrap around.
54+ """
4255function define_indices_d (:: Type{BSpline{Cubic{Periodic}}} , d, pad)
4356 symix, symixm, symixp = symbol (" ix_" ,d), symbol (" ixm_" ,d), symbol (" ixp_" ,d)
4457 symixpp, symx, symfx = symbol (" ixpp_" ,d), symbol (" x_" ,d), symbol (" fx_" ,d)
@@ -55,15 +68,16 @@ padding{BC<:Flag}(::Type{BSpline{Cubic{BC}}}) = Val{1}()
5568padding (:: Type{BSpline{Cubic{Periodic}}} ) = Val {0} ()
5669
5770"""
58- In this function we assume that `fx_d = x-ix_d` and we produce `cX_d` for
59- `X ⋹ {m, _, p, pp}` such that
71+ In `coefficients` for a cubic b-spline we assume that `fx_d = x-ix_d`
72+ and we define `cX_d` for `X ⋹ {m, _, p, pp}` such that
6073
6174 cm_d = p(fx_d)
6275 c_d = q(fx_d)
6376 cp_d = q(1-fx_d)
6477 cpp_d = p(1-fx_d)
6578
66- where `p` and `q` are defined in the docstring entry for `Cubic`.
79+ where `p` and `q` are defined in the docstring entry for `Cubic`, and
80+ `fx_d` in the docstring entry for `define_indices_d`.
6781"""
6882function coefficients {C<:Cubic} (:: Type{BSpline{C}} , N, d)
6983 symm, sym = symbol (" cm_" ,d), symbol (" c_" ,d)
@@ -82,15 +96,16 @@ function coefficients{C<:Cubic}(::Type{BSpline{C}}, N, d)
8296end
8397
8498"""
85- In this function we assume that `fx_d = x-ix_d` and we produce `cX_d` for
86- `X ⋹ {m, _, p, pp}` such that
99+ In `gradient_coefficients` for a cubic b-spline we assume that `fx_d = x-ix_d`
100+ and we define `cX_d` for `X ⋹ {m, _, p, pp}` such that
87101
88102 cm_d = p'(fx_d)
89103 c_d = q'(fx_d)
90104 cp_d = q'(1-fx_d)
91105 cpp_d = p'(1-fx_d)
92106
93- where `p` and `q` are defined in the docstring for `Cubic`.
107+ where `p` and `q` are defined in the docstring entry for `Cubic`, and
108+ `fx_d` in the docstring entry for `define_indices_d`.
94109"""
95110function gradient_coefficients {C<:Cubic} (:: Type{BSpline{C}} , d)
96111 symm, sym, symp, sympp = symbol (" cm_" ,d), symbol (" c_" ,d), symbol (" cp_" ,d), symbol (" cpp_" ,d)
147162# Prefiltering #
148163# ------------ #
149164
165+ """
166+ `Cubic`: continuity in function value, first and second derivatives yields
167+
168+ 2/3 1/6
169+ 1/6 2/3 1/6
170+ 1/6 2/3 1/6
171+ ⋱ ⋱ ⋱
172+ """
150173function inner_system_diags {T,C<:Cubic} (:: Type{T} , n:: Int , :: Type{C} )
151174 du = fill (convert (T, SimpleRatio (1 , 6 )), n- 1 )
152175 d = fill (convert (T, SimpleRatio (2 , 3 )), n)
@@ -155,8 +178,8 @@ function inner_system_diags{T,C<:Cubic}(::Type{T}, n::Int, ::Type{C})
155178end
156179
157180"""
158- `Flat` `OnGrid` amounts to setting `y_0 '(x) = 0` at `x = 0`. Applying this
159- condition gives:
181+ `Cubic{ Flat} ` `OnGrid` amounts to setting `y_1 '(x) = 0` at `x = 1`.
182+ Applying this condition yields
160183
161184 -cm + cp = 0
162185"""
@@ -173,10 +196,18 @@ function prefiltering_system{T,TC}(::Type{T}, ::Type{TC}, n::Int,
173196end
174197
175198"""
176- `Flat` `OnCell` amounts to setting `y_0'(x) = 0` at `x = -1/2`. Applying this
177- condition gives:
199+ `Cubic{Flat}`, `OnCell` amounts to setting `y_1'(x) = 0` at `x = 1/2`.
200+ Applying this condition yields
201+
202+ -9/8 cm + 11/8 c - 3/8 cp + 1/8 cpp = 0
203+
204+ or, equivalently,
178205
179206 -9 cm + 11 c -3 cp + 1 cpp = 0
207+
208+ (Note that we use `y_1'(x)` although it is strictly not valid in this domain; if we
209+ were to use `y_0'(x)` we would have to introduce new coefficients, so that would not
210+ close the system. Instead, we extend the outermost polynomial for an extra half-cell.)
180211"""
181212function prefiltering_system {T,TC} (:: Type{T} , :: Type{TC} , n:: Int ,
182213 :: Type{Cubic{Flat}} , :: Type{OnCell} )
@@ -198,10 +229,14 @@ function prefiltering_system{T,TC}(::Type{T}, ::Type{TC}, n::Int,
198229end
199230
200231"""
201- `Line` `OnCell` amounts to setting `y_0 ''(x) = 0` at `x = - 1/2`. Applying this
202- condition gives:
232+ `Cubic{ Line} ` `OnCell` amounts to setting `y_1 ''(x) = 0` at `x = 1/2`.
233+ Applying this condition yields
203234
204235 3 cm -7 c + 5 cp -1 cpp = 0
236+
237+ (Note that we use `y_1'(x)` although it is strictly not valid in this domain; if we
238+ were to use `y_0'(x)` we would have to introduce new coefficients, so that would not
239+ close the system. Instead, we extend the outermost polynomial for an extra half-cell.)
205240"""
206241function prefiltering_system {T,TC} (:: Type{T} , :: Type{TC} , n:: Int ,
207242 :: Type{Cubic{Line}} , :: Type{OnCell} )
@@ -223,13 +258,10 @@ function prefiltering_system{T,TC}(::Type{T}, ::Type{TC}, n::Int,
223258end
224259
225260"""
226- `Line` `OnGrid` amounts to setting `y_0 ''(x) = 0` at `x = 0 `. Applying this
261+ `Cubic{ Line} ` `OnGrid` amounts to setting `y_1 ''(x) = 0` at `x = 1 `. Applying this
227262condition gives:
228263
229- 1 cm -2 c + 1 cp = 0 = 0
230-
231- This is the same system as `Quadratic{Line}`, `OnGrid` so we reuse the
232- implementation
264+ 1 cm -2 c + 1 cp = 0
233265"""
234266function prefiltering_system {T,TC} (:: Type{T} , :: Type{TC} , n:: Int ,
235267 :: Type{Cubic{Line}} , :: Type{OnGrid} )
@@ -247,6 +279,14 @@ function prefiltering_system{T,TC}(::Type{T}, ::Type{TC}, n::Int,
247279 Woodbury (lufact! (Tridiagonal (dl, d, du), Val{false }), specs... ), zeros (TC, n)
248280end
249281
282+ """
283+ `Cubic{Periodic}` `OnGrid` closes the system by looking at the coefficients themselves
284+ as periodic, yielding
285+
286+ c0 = c(N+1)
287+
288+ where `N` is the number of data points.
289+ """
250290function prefiltering_system {T,TC,GT<:GridType} (:: Type{T} , :: Type{TC} , n:: Int ,
251291 :: Type{Cubic{Periodic}} , :: Type{GT} )
252292 dl, d, du = inner_system_diags (T,n,Cubic{Periodic})
@@ -260,14 +300,11 @@ function prefiltering_system{T,TC,GT<:GridType}(::Type{T}, ::Type{TC}, n::Int,
260300end
261301
262302"""
263- The free boundary condition for either `OnGrid` or `OnCell` makes sure the
264- interpoland has a continuous third derivative at the second-to-outermost cell
265- boundary: `y_0'''(1) = y_1'''(1)` and `y_{n-1}'''(n) = y_n'''(n)`. Applying this
266- condition gives:
303+ `Cubic{Free}` `OnGrid` and `Cubic{Free}` `OnCell` amount to requiring an extra
304+ continuous derivative at the second-to-last cell boundary; this means
305+ `y_1'''(2) = y_2'''(2)`, yielding
267306
268307 1 cm -3 c + 3 cp -1 cpp = 0
269-
270- This is the same system as `Quadratic{Free}` so we reuse the implementation
271308"""
272309function prefiltering_system {T,TC,GT<:GridType} (:: Type{T} , :: Type{TC} , n:: Int ,
273310 :: Type{Cubic{Free}} , :: Type{GT} )
0 commit comments