The central theme of this pair of papers (Parts I and II in this issue) is self-similarity, which is used as a bridge for connecting splines and fractals. The first part of the investigation is deterministic, and the context is that of L-splines; these are defined in the following terms: s (t) is a cardinal L-spline iff L{s(t)} = Sigma(k is an element of Z) a[k]delta(t - k), where L is a suitable pseudodifferential operator. Our starting point for the construction of "self-similar" splines is the identification of the class of differential operators L that are both translation and scale invariant. This results into a two-parameter family of generalized fractional derivatives, partial derivative(gamma)(tau), where gamma is the order of the derivative and tau is an additional phase factor. We specify the corresponding L-splines, which yield an extended class of fractional splines. The operator partial derivative(gamma)(tau),Y is used to define a scale-invariant energy measure-the squared L-2-norm of the gamma th derivative of the signal-which provides a regularization functional for interpolating or fitting the noisy samples of a signal. We prove that the corresponding variational (or smoothing) spline estimator is a cardinal fractional spline of order 2 gamma, which admits a stable representation in a B-spline basis. We characterize the equivalent frequency response of the estimator and show that it closely matches that of a classical Butterworth filter of order 2 gamma. We also establish a formal link between the regularization parameter lambda and the cutoff frequency of the smoothing spline filter: omega(0) approximate to lambda(-2 gamma). Finally, we present an efficient computational solution to the fractional smoothing spline problem: It uses the fast Fourier transform and takes advantage of the multiresolution properties of the underlying basis functions.