Researchers in a variety of biomedical fields have utilized frequency domain properties of heart rate variability (HRV), or the elapsed time between consecutive heart beats. HRV is measured from the electrocardiograph signal through the interbeat interval series. Popular approaches for estimating power spectra from these interval data apply common spectral analysis methods that are designed for the analysis of evenly sampled time series. The application of these methods to the interbeat interval series, which is indexed over an uneven time grid, requires a bias-inducing transformation. The goal of this article is to explore the use of penalized sum of squares for nonparametric estimation of the spectrum of HRV directly from the interbeat intervals. A novel cross-validation procedure is introduced for smoothing parameter selection. Empirical properties of the proposed estimation procedure are explored and compared with popular methods in a simulation study. The proposed method is used in an analysis of data from an insomnia study, which seeks to illuminate the association between the power spectrum of HRV during different periods of sleep with response to behavioral therapy.