Lecture notes. Written from the material of Berkeley Stat 153, licensed CC BY 4.0. These are notes, not a transcript: the material has been reorganised and rewritten. This adaptation carries the same licence, and the original is linked above.
94. Frequency Estimation and Aliasing¶
What this covers¶
This chapter works through a lab problem: given data assumed to follow a single sinusoid plus noise, how do you estimate the frequency, and how much does it matter that the observation times are equally spaced? It builds a least-squares (frequentist) estimate of the frequency and a Bayesian posterior for it, both by grid search rather than a closed-form formula, and then examines aliasing — the fact that equally spaced sampling makes some frequencies indistinguishable from others. It assumes familiarity with ordinary least squares regression; the specific posterior formula used below is quoted as something already derived in lecture, and this chapter does not re-derive it.
The model¶
The data \(y_1,\dots,y_n\) are generated from
for fixed but unknown \(\beta_0,\beta_1,\beta_2,f\) and noise level \(\sigma\), with \(\epsilon_t\) independent mean-zero noise. The lab simulates such a series with \(f = 0.2\), \(n = 400\), \(\beta_0 = 0\), \(\beta_1 = 3\), \(\beta_2 = 5\), \(\sigma = 10\), and the task is to recover \(f\) — together with uncertainty about it — from \(y\) alone.
The awkward feature of this model is that it is linear in \(\beta_0,\beta_1,\beta_2\) but nonlinear in \(f\): once \(f\) is fixed, \(\cos(2\pi f t)\) and \(\sin(2\pi f t)\) are just two more columns of a design matrix, and estimating the \(\beta\)'s is ordinary linear regression with \(p = 3\) parameters. There is no equivalent shortcut for \(f\) itself, so it has to be searched for.
Estimating f by least squares¶
Profile out the linear parameters: for each candidate value of \(f\), fit \(\beta_0,\beta_1,\beta_2\) by OLS on the design matrix with columns \(1\), \(\cos(2\pi f t)\), \(\sin(2\pi f t)\), and record the resulting residual sum of squares. This defines a criterion function
Minimizing \(\text{crit}(f)\) over \(f\) gives the least-squares (and, under Gaussian noise, maximum likelihood) estimate \(\hat f\). Because \(RSS(f)\) has no closed form minimizer in \(f\), this minimization is done by evaluating \(\text{crit}(f)\) on a fine grid over \(f\) and taking the arg-min — a grid of \(100{,}000\) points on \([0, 0.5]\) in the lab. Once \(\hat f\) is found, the estimates of \(\beta_0,\beta_1,\beta_2\) and \(\sigma\) are the ordinary regression estimates obtained by plugging \(f = \hat f\) back into the design matrix and refitting.
For the simulated series above this gives \(\hat f \approx 0.19990\), extremely close to the true \(f = 0.2\). Refitting at \(\hat f\) gives coefficient estimates \((\hat\beta_0,\hat\beta_1,\hat\beta_2) \approx (0.98, 3.63, 4.11)\) — all within the regression's own 95% confidence intervals of the true values \((0, 3, 5)\) — and \(\hat\sigma \approx 9.50\), close to the true \(\sigma = 10\).
Why the search only needs to run over [0, 1/2]¶
The grid above only covers \(f \in [0, 0.5]\), and that is not an arbitrary choice of range — it is forced by the fact that \(t\) only ever takes integer values \(1,\dots,n\). For any integer \(k\),
whenever \(t\) is an integer, since \(2\pi k t\) is a whole multiple of \(2\pi\). And because cosine is even and sine is odd,
Putting the two together: for any integer \(k\), the frequencies \(f\) and \(k - f\) produce identical \(\cos\) columns and a sign-flipped \(\sin\) column in the design matrix, at every integer \(t\). So they give exactly the same residual sum of squares (with \(\beta_2 \mapsto -\beta_2\)), and the least squares criterion cannot tell them apart. This is aliasing: \(\text{crit}(f)\) is periodic with period 1 and symmetric about every half-integer, so its shape on \([0, \tfrac12]\) already determines it everywhere else, and it is enough to search that one interval.
The lab makes this concrete with \(k = 2\): the frequency \(f = 1.8\) and \(f_0 = 2 - 1.8 = 0.2\) are aliases of one another, in the sense that the two sinusoids \(s_t = R\cos(2\pi f t + \phi)\) and \(s_t = R\cos(2\pi f_0 t - \phi)\) are identical whenever \(t\) is an integer. Concretely: if the true frequency is changed from \(f = 0.2\) to \(f = 1.8\) (keeping everything else, including the equally spaced sample times, the same), the least-squares and Bayesian analyses both still return an estimate close to \(0.2\) — because that is the only frequency in \([0, 0.5]\) that the data, sampled at integer times, can express. Consistent with the sign flip above, the estimate of \(\beta_2\) comes out with the opposite sign from what it would have had at \(f = 0.2\).
Breaking aliasing with unequally spaced samples¶
The derivation above depended entirely on \(t\) being an integer — that is what let \(2\pi k t\) collapse to a multiple of \(2\pi\). If the observation times are not equally spaced, that cancellation no longer happens for a generic (non-integer) \(t\), and frequencies outside \([0, 0.5]\) become identifiable again.
The lab demonstrates this by re-simulating the same model with \(f = 1.8\), \(n = 400\), and the same \(\beta_0,\beta_1,\beta_2,\sigma\), but this time drawing the \(n\) observation times as \(n\) points sampled uniformly at random from \([1, 400]\) and sorted, rather than as \(1,\dots,n\). With time no longer restricted to integers, the criterion function is evaluated on a grid over the wider range \([0, 10]\) instead of \([0, 0.5]\), since there is no longer a reason to expect the answer to lie in \([0, 0.5]\). The resulting least-squares estimate is \(\hat f \approx 1.7996\), correctly recovering the true frequency rather than its alias, and the Bayesian posterior (below) is likewise tightly concentrated near \(1.8\) rather than near \(0.2\).
A Bayesian posterior for f¶
The course had already derived, for this model, the (unnormalized) posterior density of \(f\):
where \(X_f\) is the design matrix with columns \(1, \cos(2\pi f t), \sin(2\pi f t)\) at the candidate frequency \(f\), \(p = 3\) is the number of regression coefficients, and \(S(\hat\beta(f), f)\) is the residual sum of squares at \(f\) with \(\beta\) profiled out — exactly the \(\text{crit}(f) = RSS(f)\) from the least-squares section. (The indicator restricting \(f\) to \([0, 1/2]\) is exactly the aliasing restriction derived above; when sampling is non-uniform, as in the previous section, the same argument extends the range the grid is evaluated over.)
Because \(RSS(f)\) is raised to the power \(-(n-p)/2\) — a large negative power for \(n = 400\) — the posterior is extremely sensitive to small differences in \(RSS(f)\), which is also why the posterior mode turns out to coincide with the least-squares/MLE estimate \(\hat f\): both are found by locating where \(RSS(f)\) is smallest.
Evaluating this on a grid raises two numerical issues, both handled by the lab:
- Work with the log posterior. Directly exponentiating can overflow or underflow, so the grid
search computes \(\log|X_f^\top X_f|^{-1/2} + \tfrac{p-n}{2}\log RSS(f)\) instead, using
np.linalg.slogdetfor the log-determinant term. - Exclude frequencies very close to the boundary. As \(f \to 0\) or \(f \to 1/2\), the columns of \(X_f\) become nearly collinear (with the constant column, or with each other) and \(X_f^\top X_f\) becomes nearly singular, so \(|X_f^\top X_f|^{-1/2}\) blows up. The lab trims the first and last 100 points of a \(100{,}000\)-point grid before evaluating the log posterior, rather than let those spurious spikes dominate.
Once the log posterior values are in hand, exponentiate after subtracting the maximum value (so the largest term becomes \(1\) rather than something that might overflow), then normalize by the sum over the grid to get posterior probabilities, and rescale by the number of grid points per unit interval to get an (approximate) posterior density.
From here, the posterior mean and mode of \(f\) can be read off directly, and a credible interval is built by growing a window symmetrically around the posterior mode: starting from a window of width \(0\), widen it by one grid point at a time until the posterior probability captured inside it first reaches \(95\%\).
For the original \(f = 0.2\) simulation, this gives a posterior mean \(\approx 0.19990\) and posterior mode exactly equal to \(\hat f \approx 0.19990\), with a \(95\%\) credible interval of roughly \([0.19957, 0.20022]\) — a very narrow interval that comfortably contains the true value \(0.2\). For the non-uniformly sampled \(f = 1.8\) simulation, the same procedure gives a posterior mean \(\approx 1.79959\) and mode \(\approx 1.79962\), again close to the true frequency.
Sources¶
- Berkeley STAT 153 (Spring 2025), Lab 3, "More on fitting sinusoidal models," part 1: Estimating Frequency in a Simulated Dataset — the model, the least-squares/MLE grid search for \(\hat f\), and the Bayesian posterior computation, credible interval, and numerical results quoted above.
- Same lab, part 2: Frequency Aliasing — the aliasing identity, the \(f = 1.8\) example, the non-uniform-time simulation and its recovered estimate.
- The posterior formula used in the Bayesian section is quoted in the lab as "the formula for the Bayesian posterior that we derived in class"; the derivation itself (the prior assumed, and how \(\beta\) and \(\sigma\) are integrated out) is not part of the supplied material and is not reconstructed here.
← 93. Anatomy of a Regression Fit (part 2) · Contents · 95. Inference in Nonlinear Regression Models →