Fitting Sine and Cosine Models to Periodic Data

From a Periodic Pattern to a Linear Model

Periodic measurements can be modelled using sine and cosine functions at a chosen angular frequency. The times at which measurements were taken determine the values of these functions, so the model becomes a matrix equation in its unknown coefficients. This is useful even when the observations are noisy: a least-squares fit selects coefficients whose predicted values are as close as possible to the observed data in the sum-of-squared-errors sense.

Suppose measurements \(y_1,\dots,y_m\) are recorded at times \(t_1,\dots,t_m\), and choose an angular frequency \(\omega\in\mathbb{R}\). We will fit an intercept together with a sine and cosine term. The intercept allows the average level of the data to differ from zero. The coefficients \(a\) and \(b\) determine the size and timing of the periodic variation.

The frequency is chosen before estimating \(c,a,b\). With \(\omega\ne0\), the sine and cosine components repeat after the positive period \(T=2\pi/|\omega|\). If \(\omega=0\), their values are constant in time, so the three columns in this model cannot in general be separately identified. Even for nonzero \(\omega\), the chosen observation times may fail to distinguish the columns; we will see this directly in an example.

Finding the Least-Squares Coefficients

For a matrix \(\mathbf{A}\) and data vector \(\mathbf{y}\), write the residual as \(\mathbf{r}=\mathbf{y}-\mathbf{A}\mathbf{x}\). The least-squares condition can be expressed as the normal equations \(\mathbf{A}^{\top}\mathbf{A}\mathbf{x}=\mathbf{A}^{\top}\mathbf{y}\). In this application, each row of \(\mathbf{A}\) gives the intercept, sine and cosine values at one observation time. Thus the coefficients are found by solving a linear system, even though the model is periodic as a function of time.

Proof. Set \(\mathbf{r}=\mathbf{y}-\mathbf{A}\mathbf{x}\). First suppose that \(\mathbf{x}\) is a least-squares coefficient vector. Let \(\mathbf{a}_j\) be column \(j\) of \(\mathbf{A}\), and change only coefficient \(x_j\) by an amount \(s\in\mathbb{R}\). The resulting residual is \(\mathbf{r}-s\mathbf{a}_j\), so its squared norm is \(\|\mathbf{r}-s\mathbf{a}_j\|^2=\|\mathbf{r}\|^2-2s\mathbf{a}_j^{\top}\mathbf{r}+s^2\|\mathbf{a}_j\|^2\). This expression has a minimum at \(s=0\). If \(\mathbf{a}_j^{\top}\mathbf{r}\ne0\), choosing \(s\) with the same sign as \(\mathbf{a}_j^{\top}\mathbf{r}\) and sufficiently small nonzero magnitude makes the linear term negative with larger magnitude than the quadratic term; the expression then decreases. Therefore \(\mathbf{a}_j^{\top}\mathbf{r}=0\). This holds for every column, so \(\mathbf{A}^{\top}\mathbf{r}=\mathbf{0}\). Conversely, suppose \(\mathbf{A}^{\top}\mathbf{r}=\mathbf{0}\). For any change \(\mathbf{d}\in\mathbb{R}^p\), the new residual is \(\mathbf{r}-\mathbf{A}\mathbf{d}\), and \(\|\mathbf{r}-\mathbf{A}\mathbf{d}\|^2=\|\mathbf{r}\|^2-2\mathbf{r}^{\top}\mathbf{A}\mathbf{d}+\|\mathbf{A}\mathbf{d}\|^2=\|\mathbf{r}\|^2+\|\mathbf{A}\mathbf{d}\|^2\ge\|\mathbf{r}\|^2\). Thus \(\mathbf{x}\) minimizes the squared residual norm. Finally, expanding \(\mathbf{A}^{\top}(\mathbf{y}-\mathbf{A}\mathbf{x})=\mathbf{0}\) gives the stated normal equations. \(\square\)

Least-squares solutions exist for finite data; this is the projection result from the earlier material on residual geometry. The criterion above says exactly which coefficient vectors are solutions. A further useful distinction is whether the coefficients are unique. For a fixed design matrix, they are unique precisely when its columns are linearly independent.

Proof. Suppose the columns are linearly independent, so the only \(\mathbf{v}\in\mathbb{R}^p\) satisfying \(\mathbf{A}\mathbf{v}=\mathbf{0}\) is \(\mathbf{v}=\mathbf{0}\). If \(\mathbf{x}\) and \(\mathbf{z}\) are both least-squares solutions, the normal-equation criterion gives \(\mathbf{A}^{\top}(\mathbf{y}-\mathbf{A}\mathbf{x})=\mathbf{0}\). The residual comparison in that criterion, applied to the change \(\mathbf{d}=\mathbf{z}-\mathbf{x}\), shows that \(\|\mathbf{y}-\mathbf{A}\mathbf{z}\|^2=\|\mathbf{y}-\mathbf{A}\mathbf{x}\|^2+\|\mathbf{A}(\mathbf{z}-\mathbf{x})\|^2\). Both solutions have the minimum residual norm, so \(\|\mathbf{A}(\mathbf{z}-\mathbf{x})\|^2=0\). Hence \(\mathbf{A}(\mathbf{z}-\mathbf{x})=\mathbf{0}\), and independence implies \(\mathbf{z}-\mathbf{x}=\mathbf{0}\). Thus the solution is unique. Conversely, if the columns are dependent, there is a nonzero \(\mathbf{v}\) such that \(\mathbf{A}\mathbf{v}=\mathbf{0}\). For any least-squares solution \(\mathbf{x}\), the predictions from \(\mathbf{x}+\mathbf{v}\) equal those from \(\mathbf{x}\), since \(\mathbf{A}(\mathbf{x}+\mathbf{v})=\mathbf{A}\mathbf{x}\). Therefore \(\mathbf{x}+\mathbf{v}\) is another solution and is different from \(\mathbf{x}\). The solution cannot be unique. \(\square\)

Worked Fits at a Chosen Frequency

Interpreting Coefficients and Sampling Limits

The coefficients \(a\) and \(b\) can also be expressed as amplitude and phase. The identity \(R\sin(\omega t+\phi)=R\cos(\phi)\sin(\omega t)+R\sin(\phi)\cos(\omega t)\) shows that, when \(R=\sqrt{a^2+b^2}>0\), the fitted periodic part has amplitude \(R\) and phase \(\phi\) satisfying \(\cos(\phi)=a/R\) and \(\sin(\phi)=b/R\). Equivalently, \(\phi\) is the angle with cosine and sine given by those ratios. The intercept \(c\) is the vertical shift. If \(a=b=0\), the fitted periodic part is constant and its phase is not determined.

Common pitfall: using a signed period. A tempting but incorrect calculation is to say that \(\omega=-2\) gives period \(2\pi/\omega=-\pi\). A period is positive, and the sine and cosine components repeat after \(T=2\pi/|\omega|=2\pi/2=\pi\). In general, use the absolute value in the denominator. Also, do not infer that a visually plausible fitted curve proves the coefficients are unique: uniqueness depends on whether the design-matrix columns are independent at the actual observation times.

A fitted model is an approximation tied to the frequency and sampling scheme that were chosen. Its residuals record the differences between observed and fitted values; the normal equations ensure that those residuals have zero dot product with each design column. The model therefore balances errors against the constant, sine and cosine patterns represented by the columns. Changing the chosen frequency changes those columns and can change both the fitted coefficients and the quality of the approximation.