In the context of functional data analysis, probability density functions as non-negative functions are characterized by specific properties of scale invariance and relative scale which enable to represent them with the unit integral constraint without loss of information. On the other hand, all these properties are a challenge when the densities need to be approximated with spline functions, including construction of the respective spline basis. The Bayes space methodology of density functions enables to express them as real functions in the standard $L^2$ space using the centered log-ratio transformation. The resulting functions satisfy the zero integral constraint. This is a key to propose a new spline basis, holding the same property, and consequently to build a new class of spline functions, called compositional splines, which can approximate probability density functions in a consistent way. The paper provides also construction of smoothing compositional splines and possible orthonormalization of the spline basis which might be useful in some applications. Finally, statistical processing of densities using the new approximation tool is demonstrated in case of simplicial functional principal component analysis with anthropometric data.