Machine Learning & Signals Learning

Part II Linear Least Squares in Signals and
Signal Modeling

Preface

In recent years, the convergence of machine learning (ML) and signal processing (SP) has gathered growing attention in engineering education. Students are often introduced to ML principles at an early stage, yet many advanced SP topics, ranging from linear systems and time-frequency analysis to probabilistic modeling, traditionally require multiple specialized courses [8]. Although these SP methods yield comprehensive performance insights and rigorous conclusions, teaching them can be both time-consuming and demanding.

A key bridge between basic ML concepts and advanced SP techniques is the least squares (LS) method. LS is grounded in a simple and intuitive idea: minimizing the sum of squared errors. While direct LS computations may be \(\mathcal {O}(N^3)\), and thus less efficient than typical SP methods (\(\mathcal {O}(N \log N)\) to \(\mathcal {O}(N^2)\)), the LS perspective fosters a simpler, data-driven understanding of fundamental SP tasks. For example, the estimation of sinusoidal signal parameters in noise can be introduced by viewing it purely as a regression problem, bypassing the need for more involved probabilistic analyses. Likewise, the discrete Fourier transform (DFT) can be reframed as an extension of sinusoidal parameter estimation, illustrating SP principles with real arithmetic alone.

An LS-centric viewpoint aligns well with the foundational prerequisites of many ML courses and can be integrated at an early stage of engineering or data science programs. It offers an accessible path for teaching core SP ideas to engineering students who might lack extensive mathematical or probabilistic training. Although the underlying techniques are not new, this data-driven, regression-based interpretation may be more intuitive for those already familiar with basic ML concepts, enabling them to explore SP topics with minimal additional theoretical overhead.

19 Sinusoidal Signal Analysis

19.1 Signal Preliminaries

19.1.1 Cosine Signal

A general continuous-time cosine signal can be written as

\begin{equation} \begin{aligned} x(t) &= A\cos (2\pi F_0t + \theta ),\\ &= A\cos (\Omega _0 t + \theta ), \end {aligned} \end{equation}

where

• • \(A>0\) is the amplitude,

• • \(-\pi <\theta \leq \pi \) or \(0\le \theta < 2\pi \) is the phase,

• • \(F_0\) is the frequency in Hz,

• • \(\Omega _0 = 2\pi F_0\) is the radial frequency in rad/sec,

For the further analysis, we use the sampled version \(x[n]\) of the continuous-time signal \(x(t)\), sampled with frequency \(F_s=1/T\),

\begin{equation} \label {eq:base_signal} \begin{aligned} x[n] &= x(nT)\\ &= A\cos \left (\omega _0 n+\theta \right )\quad n=0,\ldots ,L-1, \end {aligned} \end{equation}

where

\begin{equation} \label {eq-signal-omega0} \begin{aligned} \omega _0 &= 2\pi F_0 T=2\pi \frac {F_0}{F_s}\\ \end {aligned} \end{equation}

is the angular frequency (measured in radians per sample) derived from the analog frequency \(F_0\) and \(L\) is the resulting number of samples.

Nyquist criterion

In order to accurately reproduce a cosine signal, the Nyquist criterion demands \(F_0<F_s/2\), which implies \(0\le \omega _0<\pi \). This requirement can be easily illustrated by the following example. Consider two signals:

\[ x_1 (t) = \cos (0.6\pi t ), \ \ \ x_2(t) = \cos (2.6\pi t) \]

Sampling with \(F_s=1\) Hz results in two identical signals,

\begin{align*} x_1[n] &= \cos (0.6\pi n ), \\ x_2[n] &= \cos (2.6\pi n) = \cos (0.6\pi n + 2 \pi n) = x_1[n]. \end{align*} This phenomenon is called aliasing. Therefore, the sampling frequency requirement is \(0\leqslant \omega <\pi \). Note, when \(\omega _0=0\) the signal is the DC level, \(y(t)=y[n]=A\cos (\theta )\). For \(\omega _0=\pi \) the signal is “alternating sign” signal \(y[n] = A(-1)^n\).

This relation holds for all the following discussions and derivations.

19.1.2 Energy & Power

The energy of (any) signal \(x[n]\) is defined as

\begin{equation} E_\bx = \norm {\bx }^2 = \sum _{n=0}^{L-1} x^2[n], \end{equation}

where \(\bx \) is the vector of samples of the signal \(x[n]\). The corresponding power is time-averaged energy,

\begin{equation} \label {eq-signal-power} P_\bx = \frac {1}{L}E_\bx = \frac {1}{L}\norm {\bx }^2. \end{equation}

19.1.3 Noise

We define noisy signal

\begin{equation} y[n] = x[n] + \epsilon [n], \end{equation}

where \(\epsilon (t)\) is additive noise.

The only assumption for the additive noise is that it is zero-mean,

\begin{equation} \bar \epsilon = \frac {1}{L}\sum _{n=0}^{L-1}\epsilon [n] = 0. \end{equation}

No additional assumptions, such as Gaussianity, are applied; however, the special case of additive white Gaussian noise (AWGN) is further refined as tips for selected topics.

The interpretation of the power (19.5) of the zero-mean noise signal is the (biased) noise variance,

\begin{equation} P_\epsilon = \sigma _\epsilon ^2. \end{equation}

19.1.4 Basic Signal Analysis

Given the noise signal \(y[n]\), the goal is to find the parameters of \(x[n]\), as well as noise power, as illustrated in Fig. 19.1.

In the following, for each method, matrices/vectors \(\bX ,\by ,\bw \) and there relation to a sinusoidal signal are to be defined.

19.2 Amplitude estimation

19.2.1 Model

Signal model Given a signal model with \(\theta =0\) and a known frequency \(\omega _0\),

\begin{equation} x[n] = A\cos (\omega _0n) \end{equation}

the goal is to estimate the amplitude \(A\) from

\begin{equation} y[n] = x[n] + \epsilon [n]\quad n=0,\ldots ,L-1 \end{equation}

that best fits a provided model,

\begin{equation} \label {eq-amplitude-model} \hat {y}[n] = \hat {A}\cos (\omega _0n) \end{equation}

LS Formulation Technically, we are looking for the value of \(A\) that minimizes the squared error,

\begin{equation} \begin{aligned} \loss =\sum _n (y[n]-\hat {y}[n])^2. \end {aligned} \end{equation}

In the vector form, \(\bX \) is formed by samples of the signal \(\left \lbrace \cos (\omega _0n)\right \rbrace _{n=0}^{L-1}\) and \(\by \) is formed by samples of the signal \(\left \lbrace y[n]\right \rbrace _{n=0}^{L-1}\),

\begin{equation} \by = \begin{bmatrix} y[0]\\ y[1] \\ y[2]\\ \vdots \\ y[L-1] \end {bmatrix},\quad \bX = \begin{bmatrix} 1 \\ \cos (\omega _0) \\ \cos (2\omega _0) \\ \vdots \\ \cos ((L-1)\omega _0) \end {bmatrix}, \quad \bw = A \end{equation}

such that corresponding vector form of the SSE loss ((19.10)) is

\begin{equation} \loss (A) =\norm {\by -\bX A}^2 \end{equation}

The resulting minimization is straightforward,

\begin{equation} \label {eq:signal:amplitude} \begin{aligned} \hat {A} &= (\bX ^T\bX )^{-1}\bX ^T\by \\[2pt] &=\frac {\sum _n y[n]\cos (\omega _0n)}{\sum _n \cos [2](\omega _0n)} \end {aligned} \end{equation}

Interpretation If we substitute the signal model into the resulting solution,

\begin{equation} \begin{aligned} \hat {A} &=\frac {\sum _n (A\cos (\omega _0n) + \epsilon [n])\cos (\omega _0n)}{\sum _n \cos [2](\omega _0n)}\\[3pt] &=A + \frac {\sum _n \epsilon [n]\cos (\omega _0)}{\sum _n \cos [2](\omega _0n)} \end {aligned} \end{equation}

it produces a true value of \(A\) with some additive noise.

19.2.2 Signal-to-Noise-Ratio

Using the interpretation of \(\norm {\by }^2\) as the signal energy and \(\norm {\hat {\by }}^2\) as an energy of the estimated sinusoidal signal (19.15), the residual error \(\norm {\mathbf {e}}^2\) (19.12) is the estimation of the noise,

\begin{equation} \label {eq-signal-noise-orthogonality} \norm {\by }^2 = \norm {\hat {\by }}^2 + \norm {\mathbf {e}}^2\,\Rightarrow \, \begin{cases} E_\by = E_{\hat {\by }} + E_\be \\ P_\by = P_{\hat {\by }} + P_\be \end {cases} \end{equation}

Moreover, due to zero-mean property of the noise, the estimated power and/or variance of the noise (19.5) is

\begin{equation} \label {eq:signal:Anoise} \hat {P}_\epsilon = \hat {\sigma }_\epsilon ^2 = \frac {1}{L}\norm {\mathbf {e}}^2. \end{equation}

An interesting interpretation of this result is estimated signal to noise ratio (SNR), defined as

\begin{equation} \label {eq:signal:Asnr} \widehat {SNR} = \frac {P_{\hat {\by }}}{P_\be } =\frac {\norm {\hat {\by }}^2}{\norm {\mathbf {e}}^2} \end{equation}

Since the SNR is a ratio of two powers, it is commonly reported on a logarithmic scale in decibels (dB),

\begin{equation} \label {eq:signal:snr_db} SNR_{dB} = 10\log _{10}(SNR), \end{equation}

so that \(SNR=1\) (equal signal and noise power) corresponds to \(0\) dB and every negative value indicates a signal weaker than the noise.

19.2.3 Theoretical performance bounds (*)

• • For AWGN and a known frequency \(\omega _0\), the LS solution (19.19) is also the maximum likelihood (ML) estimator of the amplitude, and it is unbiased, as already shown by the interpretation above [12].

• • Its variance is available in a closed form and attains the Cramér-Rao lower bound (CRLB) with equality,

  \(\seteqnumber{0}{}{24}\)

  \begin{equation} \label {eq-signals-A_cr} \Var [\hat {A}] = \frac {\sigma ^2}{\sum _n \cos [2](\omega _0n)}\approx \frac {2\sigma ^2}{L}\quad [V^2], \end{equation}

  where the approximation uses \(\sum _n \cos [2](\omega _0n)\approx L/2\), valid for any \(\omega _0\ne 0,\pi \).

• • The two excluded frequencies are the DC level \(x[n]=A\) (\(\omega _0=0\)) and the alternating-sign signal \(x[n]=A(-1)^n\) (\(\omega _0=\pi \)). For both, the regressor energy is exactly \(\sum _n \cos [2](\omega _0n)=L\) rather than \(L/2\), so

  \(\seteqnumber{0}{}{25}\)

  \begin{equation} \label {eq-signals-A_cr_dc} \Var [\hat {A}] = \frac {\sigma ^2}{L}\qquad \omega _0=0,\pi . \end{equation}

  These are the most accurate cases: the variance is half of the value at any other frequency, i.e. the standard deviation is \(\sqrt {2}\) times smaller. For \(\omega _0=0\) the estimator (19.19) reduces to the sample mean, \(\hat {A} = \frac {1}{L}\sum _n y[n]\).

• • The practical reading of (19.25) is a relative accuracy of

  \(\seteqnumber{0}{}{26}\)

  \begin{equation} \frac {\sqrt {\Var [\hat {A}]}}{A} \approx \frac {1}{\sqrt {L\cdot SNR}}, \end{equation}

  with \(SNR = A^2/(2\sigma ^2)\). Halving the amplitude error therefore requires four times more samples, and a low SNR can always be compensated by a longer record.

• • As above, for other noise distributions or very short records, the estimation variance is readily approximated by Monte-Carlo simulations (Sec. 19.2.4).

19.2.4 Monte-Carlo Simulation

The accuracy of every estimator of this chapter is characterized by repeating the noisy experiment many times, which is the Monte-Carlo simulation of Sec. 1.6. Only its signal-specific form is recalled here.

A single realization is a noise record \(\epsilon _j[n]\) of length \(L\), added to the deterministic signal \(x[n]\),

\begin{equation} \label {eq-signal-mc-realization} y_j[n] = x[n]+\epsilon _j[n],\qquad n=0,\ldots ,L-1,\quad j=1,\ldots ,K, \end{equation}

where \(K\) is the number of realizations. Each realization yields one estimate of the parameter of interest, say the amplitude \(\hat {A}_j\), and the \(K\) collected estimates characterize the estimator by

\begin{align} \label {eq-signal-mc-summary} \text {bias}\left [\hat {A}\right ] &= \frac 1K\sum _{j=1}^{K}\hat {A}_j - A,\\ \Var [\hat {A}] &\approx \frac {1}{K-1}\sum _{j=1}^{K}\left (\hat {A}_j-\frac 1K\sum _{i=1}^{K}\hat {A}_i\right )^2, \end{align} which are then compared with the theoretical values, whenever such values exist. The identical loop applies to the phase, the frequency, or any other estimated quantity.

The closed-form expressions of the following sections are derived for AWGN and for long records. The loop above makes neither assumption: the noise in (19.28) may follow any zero-mean distribution and \(L\) may be arbitrarily short. This is why the simulation serves as the default validation tool throughout the chapter.

19.3 Amplitude and phase estimation

19.3.1 Model

Signal model The analysis is provided for the more general model,

\begin{equation} x[n] = A\cos \left (\omega _0n + \theta \right )\quad n=0,\ldots ,L-1, \end{equation}

The goal is to estimate amplitude \(A\) and the phase \(\theta \) that best fits a provided model,

\begin{equation} \hat {y}[n] = \hat {A}\cos \left (\omega _0n + \hat {\theta }\right ) \end{equation}

LS Formulation The first step involves the use of trigonometric identities to express the cosine with a phase shift as a linear combination of sine and cosine signals,

\begin{align} A\cos \left (\omega _0n + \theta \right ) &=w_c\cos (\omega _0n) + w_s\sin (\omega _0n), \end{align} where

\begin{equation} \begin{aligned} w_c &= A\cos (\theta )\\ w_s &= -A\sin (\theta ) \end {aligned} \end{equation}

and

\begin{equation} \label {eq-cos-linear-coeff} \begin{aligned} A &= \sqrt {w_c^2 + w_s^2}\\ \theta &= -\arctan (\frac {w_c}{w_s}) \end {aligned} \end{equation}

This transforms the problem into a two-parameter linear LS problem in terms of \(w_c\) and \(w_s\) [1]. In the vector form, \(\by \) is formed by samples of the signal \(\left \lbrace y[n]\right \rbrace _{n=0}^{L-1}\), the matrix \(\bX \) has dimensions \(L\times 2\), and \(\bw \) collects the two linear coefficients,

\begin{equation} \label {eq-X-single-freq} \by = \begin{bmatrix} y[0]\\ y[1] \\ y[2]\\ \vdots \\ y[L-1] \end {bmatrix},\quad \bX = \begin{bmatrix} 1 & 0 \\ \cos (\omega _0) & \sin (\omega _0) \\ \cos (2\omega _0) & \sin (2\omega _0) \\ \vdots & \vdots \\ \cos ((L-1)\omega _0) & \sin ((L-1)\omega _0) \end {bmatrix},\quad \bw = \begin{bmatrix} w_c \\ w_s \end {bmatrix}. \end{equation}

The SSE loss (19.10) is \(\loss (\bw )=\norm {\by -\bX \bw }^2\), with the LS solution given by (19.11).

Interpretation Once \(\hat {\bw }\) has been found, the amplitude and phase can be recovered from (19.35). The SNR and noise variance interpretations follow the orthogonality decomposition (19.12), similar to the previous model in Eqs. (19.23) and (19.22).

19.3.2 Theoretical performance bounds (*)

• • This estimation procedure is optimal in the maximum likelihood (ML) sense under additive white Gaussian noise (AWGN) and achieves the Cramér-Rao lower bound (CRLB) [1, 12].

• • The theoretical lower bound (also termed Cramer-Rao lower bound(CRLB)) on the variance of estimation accuracy of \(w_c,w_s\) is given by [1, Eqs. (5.47-48)]

  \(\seteqnumber{0}{}{36}\)

  \begin{equation} \label {eq-signals-Ap_cr} \Var [\hat {w}_{c,s}]\gtrsim \frac {2\sigma ^2}{L} \end{equation}

  This bound is the tightest for the AWGN case and is less accurate for other noise distributions. A numerical validation is shown in Fig. 19.4.

• • Propagating (19.37) through the parameter transformation (19.35) gives the corresponding bounds on the physical parameters [12],

  \(\seteqnumber{0}{}{37}\)

  \begin{align} \Var [\hat {A}] &\gtrsim \frac {2\sigma ^2}{L}\quad [V^2] \label {eq-signals-Aph_cr_A}\\ \Var [\hat {\theta }] &\gtrsim \frac {2\sigma ^2}{L A^2} = \frac {1}{L\cdot SNR}\quad [rad^2] \label {eq-signals-Aph_cr_ph} \end{align} with \(SNR = A^2/(2\sigma ^2)\). The amplitude bound (19.38) is the same as in the known-phase case (19.25): estimating the phase in addition costs nothing in amplitude accuracy. The phase bound, in contrast, depends on the amplitude itself, since a weak signal is harder to align in time. Both improve as \(1/L\).

• • While the derivation of CRLB is non-trivial, the approximated numerical estimation variance can be easily evaluated by Monte-Carlo simulations (Sec. 19.2.4) for any set of parameters and any distribution of interest.

19.4 Frequency estimation

19.4.1 Model

Signal model Here all of \(A,\theta ,\omega _0\) are unknown. The frequency \(\omega _0\) is estimated by searching for the \(\hat {\omega }_0\) that best fits a sinusoidal model for the observed data, i.e., that minimizes the residual error norm or maximizes the reconstructed signal energy. The corresponding matrix \(\bX \) is parameterized as a frequency-dependent one, \(\bX (\omega )\), so the estimated signal is frequency-dependent

\begin{equation} \label {eq:signal:P_frequency_dependent} \hat {\by }(\omega ) =\bX (\omega )\bw (\omega ), \end{equation}

where \(\bw (\omega )\) are the estimated parameters \(w_c\) and \(w_s\) at that frequency, as defined in (19.36).

LS Formulation The corresponding frequency-dependent residual error is given by (see also (19.21))

\begin{equation} \mathbf {e}(\omega ) = \by - \hat {\by }(\omega ). \end{equation}

Since the error \(\mathbf {e}(\omega )\) is orthogonal to \(\hat {\by }\) (19.12),

\begin{equation} \norm {\by }^2 = \norm {\hat {\by }(\omega )}^2 + \norm {\be (\omega )}^2. \end{equation}

To find the frequency that best represents the data, we seek the one that maximizes the energy of the reconstructed signal (or equivalently minimizes the residual error),

\begin{equation} \label {eq:signal:omega0} \hat {\omega }_0 = \arg \min \limits _{\omega } \norm {\mathbf {e}(\omega )}^2 = \arg \max \limits _{\omega } \norm {\hat {\mathbf {y}}(\omega )}^2. \end{equation}

Note, this optimization problem can be challenging because the objective function may exhibit multiple local maxima/minima. Therefore, an appropriate numerical global optimization method is required.

Interpretation Once \(\hat {\omega }_0\) has been found, the amplitude and phase are estimated using the corresponding linear LS solution \(\bw (\omega _0)\) (Sec. 19.3). This solution also results in SNR and noise variance estimations, as in Eqs. (19.23) and (19.22).

19.4.2 Periodogram

The function

\begin{equation} \label {eq:signal:periodogram} P(\omega ) = \frac {1}{L}\norm {\hat {\by }(\omega )}^2 \end{equation}

as a function of \(\omega \) is termed a periodogram that is a frequency-dependent quantification of signal power that approximates the power spectral density (PSD) of the signal. By scanning over frequencies, the maximum periodogram value is taken as the frequency estimate,

\begin{equation} \hat {\omega }_0 = \arg \min \limits _{\omega } P(\omega ) \end{equation}

The complexity of frequency estimation (without amplitude/phase step) can significantly reduced by noting that

\begin{equation} \arg \max \limits _{\omega } \norm {\hat {\by }(\omega )}^2 \propto \arg \max \limits _{\omega } \abs {\sum _{n=0}^{L-1}{[n]\exp (-j\omega n)}}^2 \end{equation}

However, this definition involves frequencies in the range \(0\le \omega <2\pi \); these has challenging physical interpretation.

19.4.3 Theoretical performance bounds (*)

Under AWGN assumption, theoretical SNR is given by

\begin{equation} SNR = \frac {A^2}{2\sigma ^2} \end{equation}

and the corresponding CRLB on the estimation variances are [12]

\begin{align} \Var [\hat {A}] &\geqslant \frac {2\sigma ^2}{L}\quad [V^2] \\ \Var [\hat {\omega }_0] &\geqslant \frac {12}{SNR\times L(L^2-1)}\approx \frac {12}{SNR\times L^3}\quad \left [\left (\frac {rad}{sample}\right )^2\right ] \\ \Var [\hat {\theta }] &\geqslant \frac {2(2L-1)}{SNR\times L(L+1)}\approx \frac {4}{SNR\times L}\quad [rad^2] \end{align} The amplitude bound is unchanged with respect to the known-frequency case (19.38), while the phase bound is about four times looser than (19.39): the price of estimating \(\omega _0\) is paid by the phase, which is strongly coupled to it.

For analog frequency \(F_0 = \dfrac {\omega _0}{2\pi }F_s\),

\begin{equation} \Var [F_0] = \Var [\omega _0]\left (\frac {F_s}{2\pi }\right )^2 \quad [Hz^2] \end{equation}

19.5 Harmonic Signal Analysis

19.5.1 Model

Signal model A particularly important class of signals encountered in many practical applications is the harmonic or periodic signal (also termed Fourier series). Such a signal can be expressed as a sum of cosine terms whose frequencies are integer multiples (harmonics) of a fundamental frequency \(\omega _0\).

\begin{equation} y[n] = A_0 + \sum _{m=1}^M A_m\cos (m\omega _0 n + \theta _m) = \sum _{m=\mathbf {0}}^M A_m\cos (m\omega _0 n + \theta _m), \end{equation}

where:

• • \(A_0\) is the constant (DC) component,

• • \(A_m\) and \(\theta _m\) represent the amplitude and phase of the \(m\)-th harmonic,

• • \(\omega _0\) is the fundamental angular frequency,

• • \(m\omega _0\) corresponds to the frequency of the \(m\)-th harmonic,

• • and \(M\) is the number of harmonics in the model.

The model order \(M\) (number of harmonics) is a hyper-parameter that should be chosen carefully. Too few harmonics can fail to capture essential signal structure, while too many may overfit noise. The maximum value of \(M\) is bounded by the Nyquist criterion, \(M\omega _0<\pi \).

LS Formulation Given \(\omega _0\), the model is linear in terms of the unknown parameters \(\{A_m, \theta _m\}\) for each harmonic \(m=1,\ldots ,M\). Similar to the single-frequency case, the LS matrix \(\bX \) is constructed with columns corresponding to \(\cos (m\omega _0 n)\) and \(\sin (m\omega _0 n)\) for \(m=1,\ldots ,M\), plus a column of ones for the DC component,

\begin{equation} \bX = \begin{bmatrix} 1 & \cdots & 1 & 0 & \cdots \\ 1 & \cdots & \cos (\omega _0m) & \sin (\omega _0m) & \cdots \\ 1 & \cdots & \cos (2\omega _0m) & \sin (2\omega _0m) & \cdots \\ \vdots & \vdots & \vdots & \vdots & \ddots \\ 1 & \cdots & \cos ((L-1)\omega _0m) & \sin ((L-1)\omega _0m) & \cdots \\ \end {bmatrix}, \quad \bw = \begin{bmatrix} A_0 & \cdots & w_{c,m} & w_{s,m} & \cdots \end {bmatrix}^T. \end{equation}

Each pair \((A_m,\theta _m)\) can be recovered from the LS estimated cosine and sine coefficients in the manner described for single-frequency amplitude-phase estimation from per-harmonic cosine/sine coefficient pairs \(w_{c,m},w_{s,m}\). The resulting SNR and noise variance estimates are similar to those described in the previous sections.

Fundamental frequency estimation If \(\omega _0\) is not known, the approach that is described in the frequency estimation section can also be applied here. Once \(\hat {\omega }_0\) has been determined from a maximum of the harmonic periodogram,

\begin{equation} P_h(\omega ) = \frac {1}{L}\sum _{m=1}^M\norm {\hat {\by }(m\omega )}^2, \end{equation}

the harmonic amplitudes and phases can be estimated via LS at this frequency [3].

Total harmonic distortion (THD)

Total harmonic distortion (THD) is a measure commonly used in electrical engineering, audio processing, and other fields to quantify how much the harmonic components of a signal differ from a pure sinusoid at the fundamental frequency. It is defined as the ratio of the root-sum-square of the harmonic amplitudes and the amplitude of the fundamental frequency,

\begin{equation} THD = \frac {\sqrt {\sum _{m=2}^M A_m^2}}{A_1}. \end{equation}

A lower THD value indicates that the signal is closer to a pure sinusoidal shape, whereas a higher THD signifies a stronger presence of higher-order harmonics.

19.5.2 Theoretical performance bounds (*)

• • The presented estimator, for both known and unknown \(\omega _0\), is an efficient ML estimator in the case of AWGN. It is a minimum variance unbiased estimator (MVUE) and attains the CRLB asymptotically, i.e. for a large number of samples \(L\) [3].

• • For a known \(\omega _0\) the model is linear, and the harmonic columns of \(\bX \) are (nearly) orthogonal to each other. Every harmonic therefore behaves as the isolated two-parameter problem of Sec. 19.3, with the bounds (19.38) and (19.39) applied per harmonic,

  \(\seteqnumber{0}{}{55}\)

  \begin{align} \Var [\hat {A}_m] &\gtrsim \frac {2\sigma ^2}{L}\quad [V^2] \label {eq-signals-harmonic_cr_A}\\ \Var [\hat {\theta }_m] &\gtrsim \frac {2\sigma ^2}{L A_m^2} = \frac {1}{L\cdot SNR_m}\quad [rad^2] \label {eq-signals-harmonic_cr_ph} \end{align} with the per-harmonic \(SNR_m = A_m^2/(2\sigma ^2)\). Weak harmonics have accurate amplitudes but poorly defined phases. The DC term is the single exception and is twice as accurate, \(\Var [\hat {A}_0]\gtrsim \sigma ^2/L\), because its column of ones has squared norm \(L\) instead of \(L/2\).

• • For an unknown \(\omega _0\), the whole set of harmonics contributes to the fundamental frequency estimate [3],

  \(\seteqnumber{0}{}{57}\)

  \begin{equation} \label {eq-signals-harmonic_cr} \Var [\hat {\omega }_0] \gtrsim \frac {24\sigma ^2}{L^3\sum _{m=1}^{M}m^2A_m^2} = \frac {12}{SNR_h\cdot L^3}\quad \left [\left (\frac {rad}{sample}\right )^2\right ], \end{equation}

  where

  \(\seteqnumber{0}{}{58}\)

  \begin{equation} \label {eq-signals-harmonic_snr} SNR_h = \frac {1}{2\sigma ^2}\sum _{m=1}^{M}m^2A_m^2 \end{equation}

  is the harmonic SNR, in which each harmonic is weighted by \(m^2\). This is exactly the single-tone bound of Sec. 19.4 with \(SNR\) replaced by \(SNR_h\), and it reduces to it for \(M=1\). The analog frequency bound follows as before, \(\Var [\hat {F}_0] = \Var [\hat {\omega }_0](F_s/2\pi )^2\).

• • The practical reading of (19.58) is that a harmonic model estimates the fundamental frequency far better than a single tone of the same power does. Every additional harmonic contributes \(m^2A_m^2\), so the high harmonics are the valuable ones even when their amplitudes are moderate. For \(M\) harmonics of comparable amplitude \(A\), \(\sum _m m^2A_m^2\approx A^2M^3/3\), and the frequency variance drops roughly as \(M^3\). Fitting the harmonics jointly is therefore much cheaper than lengthening the record, which only gives \(L^3\).

• • As above, for other noise distributions or short records, the estimation variance is readily approximated by Monte-Carlo simulations (Sec. 19.2.4).

19.6 Discrete Fourier Transform (DFT)

19.6.1 Definition

The discrete Fourier transform (DFT) is a special case of the harmonic model with:

• • \(N\ge L\) harmonics,

• • fundamental angular frequency \(\omega _0 = \dfrac {2\pi }{N}\).

Formally speaking, DFT is harmonic decomposition of a finite-length signal into a sum of \(N\) harmonically related sinusoids that are evenly spaced in frequency.

DFT representation states that:

• • any arbitrary, finite-time signal \(y[n]\) may be represented as a sum of sinusoidal signals,

  \(\seteqnumber{0}{}{59}\)

  \begin{equation} y[n] = \sum _{k=0}^{N-1}A_k\cos (k\frac {2\pi }{N}n+\theta _k),\quad n=0,\ldots ,L-1 \end{equation}

• • For \(N \geq L\), the DFT allows for perfect reconstruction of the signal using its harmonic representation parameters, \(\left \lbrace A_k,\theta _k\right \rbrace _{k=0}^{N-1}\).

The corresponding LS formulation follows the harmonic model above for known \(\omega _0\).

19.6.2 Properties

Zero-padding

The basic requirement is \(N=L\) harmonic signals. Sometimes, denser spectral representation with \(N>L\) is required. In this case, the corresponding time-domain signal can be interpreted as the original signal \(y[n]\) with \(N-L\) zeros padded at the end of this signal. This case it termed zero padding.

Symmetry

It is also worth noting the symmetry in the DFT,

\begin{equation} \begin{aligned} \cos ((N-k)\Delta \omega ) &= \cos (k\Delta \omega )\\ \sin ((N-k)\Delta \omega ) &= -\sin (k\Delta \omega ), \end {aligned} \end{equation}

resulting in redundant information for frequencies \(k\omega _0\) above and below \(\pi \),

\begin{equation} \begin{aligned} A_k &= A_{N-k}\\ \theta _k &= -\theta _{N-k} \end {aligned}. \end{equation}

As a result, only frequencies \(k\omega _0 \leq \pi \) need to be considered uniquely.

19.6.3 Power Spectral Density

The power of a signal of the form

\begin{equation} x_k[n] = A_k\cos (k\frac {2\pi }{N}n + \theta _k) \end{equation}

is

\begin{equation} S_\by [k]=P_{\by _k} = \frac {1}{L}\norm {\by _k}^2 = \frac {A_k^2}{2}. \end{equation}

This value is known as the power spectral density (PSD) at the frequency \(\omega = k\frac {2\pi }{N}\). The corresponding squared magnitude values \(A_k^2/2\) are known as the discrete-frequency periodogram (Eq. (19.44)), and this is the basic method for the PSD estimation of a signal. Plotting such a periodogram gives a frequency-domain representation of the signal’s power distribution, highlighting which frequencies carry the most power.

Parseval’s Theorem DFT is energy conservation transform with

\begin{equation} \label {eq:parseval} \sum _{k=0}^{N-1} A_k^2 = \frac {1}{L}\norm {\by }^2. \end{equation}

19.6.4 Advanced notes (*)

Fast Fourier transform (FFT) & Goertzel algorithm

The fast Fourier transform (FFT) algorithm efficiently computes \(Y[k]=\DFTtr {y[n]}\), providing \(A_k = |Y[k]|/N\) and \(\theta _k = \angle (Y[k])\) with significantly lower memory requirements and complexity than direct calculation in Eq. (19.76).

When only a single frequency value is of interest, Goertzel algorithm is more efficient method for the task. Moreover, it can be used for computationally effective peaking of the maximum in Eq. (19.43).

Parseval’s Theorem (Energy Preservation) Using FFT notation,

\begin{equation} A_k^2 = \abs {Y[k]}^2/L^2\Rightarrow \sum _{k=0}^{L-1} \abs {Y[k]}^2 = L\norm {\by }^2. \end{equation}

Sinusoidal signal power

For a more general case of an arbitrary \(\omega \) value, the signal of the form

\begin{equation} y[n] = A\cos (\omega _0 n) \end{equation}

has an \(\omega _0\)-dependent power. It is convenient to collect this dependence in a single quantity, the regressor energy

\begin{equation} \label {eq-cos-energy} \mathcal {E}(\omega _0,L) \triangleq \sum _{n=0}^{L-1}\cos [2](\omega _0 n) = \frac {L}{2}+\frac {\sin (L\omega _0)\cos ((L-1)\omega _0)}{2\sin (\omega _0)}, \end{equation}

which counts how much energy the cosine \(\cos (\omega _0n)\) accumulates over the \(L\) available samples. The same quantity is the denominator of the amplitude estimate (19.19), so it also determines the estimation accuracy: the more energy the regressor carries, the more accurate the estimate. With this notation the power is

\begin{equation} \begin{aligned}\label {eq-cos-power} P_\by &= \frac {A^2}{L}\mathcal {E}(\omega _0,L)\\ &=\frac {A^2}{4 L}\left (1+2 L-\frac {\sin (\omega _0-2 L \omega _0)}{\sin (\omega _0)}\right ), \end {aligned} \end{equation}

where the \(\omega _0\)-dependence results from the time-limited origin of the signal \(y[n]\). For almost every frequency the second term of (19.68) is negligible against \(L/2\), so

\begin{equation} \label {eq-cos-energy-approx} \mathcal {E}(\omega _0,L)\approx \frac {L}{2}\qquad \omega _0\ne 0,\pi , \end{equation}

and for the infinite number of samples the power converges to the continuous-time power expression,

\begin{equation} \lim \limits _{L\rightarrow \infty } P_\by \rightarrow \frac {A^2}{2}. \end{equation}

The two excluded frequencies are the DC level \(y[n]=A\) (\(\omega _0=0\)) and the alternating-sign signal \(y[n]=A(-1)^n\) (\(\omega _0=\pi \)). There every sample of \(y^2[n]\) equals \(A^2\) instead of oscillating between \(0\) and \(A^2\), so the energy is exactly \(\mathcal {E}=L\) and the power is \(A^2\), twice the value of any other frequency. The illustration of (19.69) over the whole band is presented in Fig. 19.7a, and the origin of this factor of two is shown in Fig. 19.7b.

Single frequency analysis

Consider a signal \(y_k[n]\) with a discrete frequency \(\omega _k = \dfrac {2\pi }{N}k\) is given,

\begin{equation} y_n[k] = A\cos (\omega _k n + \theta ) = A\cos (\frac {2\pi }{N}k n + \theta ) \end{equation}

To estimate amplitude and phase at this predefined frequency, we can form a matrix \(\bX \),

\[ \bX _k = \begin {bmatrix} \bx _c & \bx _s \end {bmatrix}, \]

where

\[ \bx _c = \begin {bmatrix} \cos (2\pi \frac {k}{N}\cdot 0) \\ \cos (2\pi \frac {k}{N}\cdot 1) \\ \vdots \\ \cos (2\pi \frac {k}{N}(L-1)) \end {bmatrix} \ \text {and} \ \bx _s = \begin {bmatrix} \sin (2\pi \frac {k}{N}\cdot 0) \\ \sin (2\pi \frac {k}{N}\cdot 1) \\ \vdots \\ \sin (2\pi \frac {k}{N}(L-1)) \end {bmatrix}. \]

By evaluating

\begin{equation} \bX _k^T \bX _k = \begin{bmatrix} \bx _c^T\bx _c & \bx _s^T\bx _c\\[5pt] \bx _s^T\bx _c & \bx _s^T\bx _s \end {bmatrix} \end{equation}

we find that the sine and cosine columns form an orthogonal basis for this single frequency, with for \(0<k<\frac {N}{2}\):

\begin{equation} \label {eq:orthogonality} \begin{aligned} \bx _{c}^T\bx _{c} &= \sum _{n=0}^{N-1}\cos [2](2\pi \frac {k}{N}n) = \frac {N}{2}\\ \bx _{c}^T\bx _{s} &= \sum _{n=0}^{N-1}\cos (2\pi \frac {k}{N}n)\sin (2\pi \frac {k}{N}n) = 0\\ &= \bx _s^T\bx _c\\ \bx _{s}^T\bx _{s} &= \sum _{n=0}^{N-1}\sin [2](2\pi \frac {k}{N}n) = \frac {N}{2} \end {aligned} \end{equation}

For the special cases \(k=0\) (DC) and, when \(N\) is even, \(k=N/2\) (Nyquist), the cosine column has squared norm \(N\) rather than \(N/2\) (it is all ones at DC and alternating \((-1)^n\) at Nyquist). At those same indices the sine column is identically zero.

Stacking these results for all \(k=0,\ldots ,N-1\) yields the complete DFT matrix forms a complete orthogonal basis for the \(L\)-sample signal space. The further discussion of \(\left (\bX ^T\bX \right )^{-1}\) matrix properties can be found in Examples 4.2 and 8.5 in [12].

Moreover, since \(\left (\bX _k^T\bX _k\right )^{-1}\) takes a particularly simple diagonal form,

\begin{equation} \left (\bX _k^T\bX _k\right )^{-1} = \begin{bmatrix} \frac {2}{N} & 0\\[5pt] 0 & \frac {2}{N} \end {bmatrix} \end{equation}

and the least squares solution \(\hat {\bw }\) for the parameters \(w_{c,k}\) and \(w_{s,k}\) (corresponding to amplitude and phase components at \(\omega _k\)) is

\begin{align} w_{c,k} &= \frac {2}{N}\sum _{n=0}^{L-1} y[n]\cos \left (2\pi \frac {k}{N}n\right ), \\ w_{s,k} &= \frac {2}{N}\sum _{n=0}^{L-1} y[n]\sin \left (2\pi \frac {k}{N}n\right ). \end{align} For the special cases \(k=0\) (DC) and, when \(N\) is even, \(k=N/2\) (Nyquist), the cosine column has squared norm \(N\) rather than \(N/2\), since it is all ones at DC and alternates as \((-1)^n\) at Nyquist,

\begin{equation} \label {eq:signals:dft_w_special} w_{c,0} = \frac {1}{N}\sum _{n=0}^{N-1} y[n], \qquad w_{c,N/2} = \frac {1}{N}\sum _{n=0}^{N-1} y[n](-1)^n, \qquad w_{s,0} = w_{s,N/2} = 0 \end{equation}

General orthogonality

The single-frequency result of Eq. (19.74) extends to different harmonic indices \(j\ne k\), which is what makes the full DFT matrix \(\bX ^T\bX \) orthogonal and its inverse trivial,

\begin{align} \sum _{n=0}^{N-1}\cos (2\pi \frac {j}{N}n)\cos (2\pi \frac {k}{N}n) &= \frac {N}{2}\delta \left [j-k\right ]\\ \sum _{n=0}^{N-1}\cos (2\pi \frac {j}{N}n)\sin (2\pi \frac {k}{N}n) &= 0\quad \forall j,k\\ \sum _{n=0}^{N-1}\sin (2\pi \frac {j}{N}n)\sin (2\pi \frac {k}{N}n) &= \frac {N}{2}\delta \left [j-k\right ]. \end{align}

19.6.5 Spectral Spreading and Leakage (*)

In an idealized setting, a pure cosine signal has a perfectly defined frequency representation. For instance, consider the discrete-time signal,

\begin{equation} \label {eq:signal:single_freq_cos} x[n] = A\cos (k_0\frac {2\pi }{L}n),\ k_0\in \left \lbrace 1,\cdots ,L-1 \right \rbrace \end{equation}

where \(k_0\) is the frequency index. The Fourier transform of this signal yields a single spectral component at frequency \(w_0=k_0\frac {2\pi }{L}\), such that the spectral amplitude \(A_k\) at each value of \(k\) is given by

\begin{equation} A_k = \begin{cases} \frac {A}{2} & k=k_0,N-k_0\\ 0 & \text {otherwise} \end {cases}. \end{equation}

Under these conditions, the signal’s spectral representation seems to be strictly localized at the specific frequency \(\omega _k\), with no energy distributed elsewhere in the spectrum.

However, practical scenarios deviate from this ideal case. In particular, if a denser frequency grid is employed (i.e. \(N>L\)) or the frequency varies continuously (as in Eq. (19.40)), the resulting spectral distribution can differ substantially from the discrete, single-peak ideal (Fig. 19.8). This difference arises because, in general, \(\bX (\omega )^T\bX (\omega )\) is not orthogonal as in Eq. (19.74). As a result, two effects are introduced:

• • The main frequency peak broadens, resulting in “spectral spreading”.

• • Additional frequency components emerge beyond the broadened main peak, termed “spectral leakage.”

19.7 Short-Time Fourier Transform (STFT)

All the analysis presented above, e.g. harmonic analysis and DFT, assume that the signal has time-constant harmonic components (termed stationary signal) over the entire observation interval of \(L\) samples. In many practical applications, however, the frequency content of a signal changes over time.

19.7.1 Definition

The short-time Fourier transform (STFT) addresses analysis of time-varying spectral behavior by applying the DFT to short, overlapping segments of the signal. A window function \(w[n]\) of length \(N_w\) is used to extract a local segment centered around a given time instant. The \(\ell \)-th windowed segment is defined as

\begin{equation} y_\ell [n] = w[n]\cdot y[\ell H + n],\quad n=0,\ldots ,N_w-1, \end{equation}

where \(H\) is the hop size (the number of samples between consecutive windows) and \(\ell = 0,1,2,\ldots \) is the segment (frame) index.

Each windowed segment \(y_\ell [n]\) is then analyzed using the DFT framework described above, i.e., the LS matrix \(\bX \) of dimension \(N_w\times 2N_w\) with cosine and sine columns at the DFT frequencies \(\omega _k = \frac {2\pi k}{N_w}\) is applied to the segment. The resulting DFT coefficients become time-dependent,

\begin{equation} A_k[\ell ],\quad \theta _k[\ell ],\quad k=0,\ldots ,N_w-1, \end{equation}

providing the amplitude and phase of each frequency component at each time frame \(\ell \).

The following example uses a chirp signal whose instantaneous frequency sweeps linearly from 50 Hz to 400 Hz over 1 s, sampled at \(F_s=1000\) Hz, with additive Gaussian noise (\(\sigma =0.5\)). The time-domain signal is shown in Fig. 19.9a and its spectrogram (\(N_w=128\), \(H=16\), Hann window) in Fig. 19.9b.

19.7.2 Time-frequency trade-off

The choice of the window length \(N_w\) involves a fundamental trade-off:

• • A longer window (large \(N_w\)) provides better frequency resolution (smaller \(\Delta \omega = \frac {2\pi }{N_w}\)) but poorer time localization, since each segment spans a longer time interval.

• • A shorter window (small \(N_w\)) provides better time localization but coarser frequency resolution.

This time-frequency trade-off is inherent to STFT and cannot be resolved by the choice of window function alone. Figure 19.9c illustrates this trade-off: a short window (\(N_w=32\)) provides sharp time localization but blurred frequency content, while a long window (\(N_w=256\)) yields precise frequency resolution at the cost of temporal smearing.

19.7.3 Spectrogram

The power spectral density of each segment defines a time-frequency representation known as the spectrogram,

\begin{equation} S[\ell , k] = \frac {A_k^2[\ell ]}{2}, \end{equation}

which displays the signal’s power distribution as a function of both time (frame index \(\ell \)) and frequency (index \(k\)). The spectrogram is one of the most widely used tools for visualizing non-stationary signals and serves as a basis for feature extraction in signal classification tasks.

Software reference

Common implementations of STFT and spectrogram:

• • Python — scipy.signal.stft, scipy.signal.spectrogram; librosa.stft.

• • MATLAB — stft, spectrogram, pspectrum.

Summary

The STFT extends the DFT to non-stationary signals by applying windowed DFT analysis to overlapping segments of the signal. The key parameters are:

• • Window length \(N_w\) — controls the time-frequency resolution trade-off.

• • Hop size \(H\) — determines the overlap between consecutive frames.

• • Window function \(w[n]\) — shapes the spectral characteristics of each segment (e.g., Hann, Hamming).

The resulting spectrogram \(S[\ell , k]\) provides a real-valued two-dimensional representation of the signal’s power as a function of time and frequency.

19.8 Quasi-periodic signals

A quasi-periodic signal has a clear repeating structure with an identifiable period, but no integer \(T\) satisfies \(x[n+T]=x[n]\) exactly. Successive cycles drift slowly in amplitude, phase, or local period.

Canonical forms

• • Slow amplitude/phase modulation,

  \(\seteqnumber{0}{}{85}\)

  \begin{equation} x[n] = A[n]\cos \left (\omega _0 n + \theta [n]\right ), \end{equation}

  where \(A[n]\) and \(\theta [n]\) vary on a much slower timescale than \(\omega _0\).

• • Sums of incommensurate sinusoids, e.g. \(\cos (n) + \cos (\sqrt {2}\,n)\), which look periodic locally but never exactly repeat.

• • Natural sources: voiced speech (wandering pitch), ECG (heart-rate variability), mechanical vibration under varying load.

Behavior of the chapter’s tools

• • Harmonic LS (Sec. 19.5) assumes constant \(A_m,\theta _m\) over the whole record. With slow modulation, LS returns time-averaged values and the residual systematically exceeds the noise floor.

• • DFT and periodogram (Sec. 19.4.2) show a broadened spectral line instead of a sharp delta. This spreading is of physical origin and is distinct from the windowing-induced leakage discussed in Sections. 19.6.5 and 19.7.

• • STFT (Sec. 19.7) recovers the slow variation as time-dependent coefficients \(A_k[\ell ],\theta _k[\ell ]\), trading frequency resolution for time localization.

19.9 Summary

The summary of the presented approach is shown in Table 19.1. The presented approach involves a design of matrix \(\bX \) and using LS to estimate unknown parameters. The key addressed task are as follows.

Amplitude Estimation With a known frequency \(\omega _0\), the amplitude \(A\) is found via LS. The resulting residuals provide noise variance and SNR estimates.

Amplitude and Phase Estimation For known \(\omega _0\), rewriting

\[ A\cos (\omega _0 n + \theta ) = w_c\cos (\omega _0 n) + w_s\sin (\omega _0 n) \]

transforms the problem into a two-parameter LS regression.

Frequency Estimation: If \(\omega _0\) is unknown, it is found by searching for the frequency that maximizes the fitted signal energy.

Harmonic Signal Analysis: Signals can be expressed as sums of multiple harmonics. Extending the LS approach to multiple harmonics allows estimation of each amplitude and phase. THD quantifies deviations from a pure tone.

Discrete Fourier Transform (DFT): The DFT is a special case of harmonic modeling, decomposing a signal into equally spaced frequency components. Efficiently computed by the FFT, the DFT is central to signal spectral analysis.

STFT: For non-stationary signals, the STFT applies the DFT to short, overlapping windowed segments, producing time-dependent spectral coefficients. The resulting spectrogram provides a two-dimensional time-frequency representation of the signal’s power distribution. The window length \(N_w\) controls the trade-off between time and frequency resolution.

Performance bounds For AWGN, all the estimators above are ML estimators and attain the Cramer-Rao lower bound (CRLB), so the achievable accuracy is known in closed form [12, 10]. The bounds collected along the chapter are summarized in Table 19.2. Three rules follow from it.

• • Amplitude accuracy is \(2\sigma ^2/L\) in every case. Estimating the phase, the frequency, or additional harmonics on top of it costs nothing in amplitude accuracy. The only exception is a constant (DC) or an alternating-sign signal, \(\omega _0=0,\pi \), whose amplitude is estimated twice as accurately, \(\sigma ^2/L\). The same holds for the DC term \(A_0\) of a harmonic model.

• • Phase accuracy scales as \(1/(L\cdot SNR)\) and frequency accuracy as \(1/(L^3\cdot SNR)\). Doubling the record length halves the phase error, but reduces the frequency error by a factor of about eight. Record length is thus far more valuable for frequency than for amplitude or phase.

• • For a harmonic signal only the frequency bound improves, through the \(m^2\)-weighted harmonic SNR (19.59): high harmonics sharpen the fundamental frequency estimate, roughly as \(M^3\).

These expressions hold for AWGN. For other noise distributions or short records, the bias and the variance of the estimators are readily approximated via Monte-Carlo simulations (Sec. 19.2.4) under any parameter setting of interest.

19.10 Exercise

1 (Amplitude estimation accuracy) Validate the amplitude estimation bound (19.25) by Monte-Carlo simulation, following the loop of Sec. 19.2.4. Use \(y[n] = A\cos (\omega _0n)+\epsilon [n]\) with AWGN, e.g. \(A=2\) and \(\sigma =2\) (\(SNR=-3\) dB), and estimate \(\hat {A}\) by (19.19). For every parameter point, average over at least \(100\) independent noise realizations.

• 1. Verify that the estimator is unbiased.

• 2. Fix \(L=128\) and sweep \(\omega _0\) over the whole band \((0,\pi ]\). Compare the sample standard deviation of \(\hat {A}\) with \(\sigma /\sqrt {\mathcal {E}(\omega _0,L)}\) and with the approximation \(\sigma \sqrt {2/L}\), where the regressor energy \(\mathcal {E}(\omega _0,L)\) is given by (19.68). An example of the expected result is shown in Fig. 19.10a.

• 3. Fix \(\omega _0\) and sweep \(L\) from \(8\) to \(2048\). Show that the result is a straight line of slope \(-1/2\) on log-log axes, as in Fig. 19.10b.

• 4. Explain what happens near \(\omega _0=0\) and \(\omega _0=\pi \)?

Hint: use the closed form of \(\mathcal {E}(\omega _0,L)\) given in (19.68), and keep in mind that it is \(0/0\) at \(\omega _0=0,\pi \), where the sum has to be evaluated directly.