Data-Driven Time-Series Analysis

Chapter 3 Basic Signal Analysis

3.1 Signal Preliminaries

A general continuous-time cosine signal can be written as

\begin{equation} \begin{aligned} y(t) &= A\cos (2\pi F_0t + \theta )+ \epsilon (t),\\ &= A\cos (\Omega _0 t + \theta )+ \epsilon (t), \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,

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

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} y[n] &= y(nT)\\ &= A\cos \left (\omega _0 n+\theta \right ) + \epsilon [n]\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.

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.

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. 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\). Therefore, the sampling frequency requirement is \(0\leqslant \omega <\pi \). This relation holds for all the following discussions and derivations.

The energy of the 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} P_\bx = \frac {1}{L}E_\bx = \frac {1}{L}\norm {\bx }^2. \end{equation}

3.2 Amplitude estimation

Given a signal model with a known frequency \(\omega _0\),

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

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

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

where \(x[n]=\cos \left (\omega _0n\right )\). 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 linear LS regression formulation, we define the corresponding parameters \(\by ,\bX \) and \(\bw \). First, the required weight is \(\bw =A\). The matrix \(\bX \) is formed by samples of the signal \(\left \lbrace x[n]\right \rbrace _{n=0}^{L-1}\),

\begin{equation} \begin{aligned} \bX &= \begin{bmatrix} x[0] & x[1] & \cdots & x[L-1]\\ \end {bmatrix}^T\\ &= \begin{bmatrix} 1 & \cos (\omega _0) & \cos (2\omega _0) & \cdots & \cos ((L-1)\omega _0) \end {bmatrix}^T \end {aligned} \end{equation}

resulting in

\begin{equation} \hat \by = \bX \bw \end{equation}

Finally, \(\by \) is the vector of samples of \(\left \lbrace y[n]\right \rbrace _{n=0}^{L-1}\). 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 x[n]y[n]}{\sum _n x^2[n]}\\[2pt] &=\frac {\sum _n y[n]\cos (\omega _0n)}{\sum _n \cos [2](\omega _0n)} \end {aligned} \end{equation}

This expression be further approximated by

\begin{equation} \hat {A} \approx \frac {2}{L}\sum _n y[n]\cos (\omega _0n) \end{equation}

by noting that

\begin{equation} \sum _n \cos [2](\omega _0n)\approx \frac {L}{2}\qquad \omega _0\ne 0,\pi \end{equation}

If we substitute the model into the resulting solution,

\begin{equation} \begin{aligned} \hat {A} &=\frac {\sum _n x[n]\left (Ax[n]+\epsilon [n]\right )}{\sum _n x^2[n]}\\[4pt] &=A + \frac {\sum _n x[n]\epsilon [n]}{\sum _n x^2[n]}, \end {aligned} \end{equation}

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

The resulting residual error is given by

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

Since \(\mathbf {e}\perp \hat {\by }\) the power/energy terms can be decomposed as follows,

\begin{equation} \norm {\by }^2 = \norm {\hat {\by }}^2 + \norm {\mathbf {e}}^2,\quad E_\by = E_{\hat {\by }} + E_\be . \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 {\norm {\hat {\by }}^2}{\norm {\mathbf {e}}^2} \end{equation}

Moreover, due to zero-mean property of the noise, the estimated (biased) variance of the noise is

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

The following example (Fig. 3.1) uses a synthetic cosine signal of length \(L = 51\) samples, angular frequency \(\omega _0 = 0.1\pi \) and amplitude \(A = 1.5\). Gaussian noise with standard deviation \(\sigma = 5\) is then added to create a noisy observation. A least-squares regression is applied to estimate the amplitude, yielding \(\hat {\sigma }_\epsilon =4.43\).

3.3 Amplitude and phase estimation

The following analysis is provided for the more general model,

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

with two unknown parameters, the amplitude \(A\) and the phase \(\theta \).

The linear LS reformulation of the signal model

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

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]. The resulting LS formulation involves a two-valued vector of linear coefficients, \(\bw = \begin {bmatrix} w_c & w_s\end {bmatrix}^T\), the vector \(\by \) of samples of \(y[n]\), and the matrix \(\bX \) of dimensions \(L\times 2\) that is given by

\begin{equation} \label {eq-X-single-freq} \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}. \end{equation}

Once \(\hat {\bw }\) has been found, the amplitude and phase can be recovered from (3.24).

SNR and noise variance interpretations are similar to in the previous model in Eqs. (3.18) and (3.19).

The numerical example is presented in Fig. 3.2. The configuration is similar to the previous figure, expect the lower noise variance, \(\sigma =1.5\). Nevertheless, there is a decrease in performance, since two parameters are estimated simultaneously.

Implementation Tip

• • 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, 6].

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

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

  \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.

• • The approximated estimation variance can be easily evaluated by Monte-Carlo simulations for any set of parameters and any distribution of interest.

3.4 Frequency estimation

Since \(\omega _0\) is unknown, the matrix \(\bX \) may be parameterized as a frequency-dependent one, \(\bX (\omega )\). Here, 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. The corresponding frequency-dependent residual error is given by

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

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

\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), as mentioned above

\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.

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

Tip: Interpretation in Terms of the 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 measure of signal power that approximates the power spectral density (PSD) of the signal. By scanning over frequencies, the \(\omega \) that yields the maximum periodogram value is taken as the frequency estimate, \(\omega _0\).

The complexity of frequency estimation (without amplitude/phase) 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}

A numerical example of the signal with additive white Gaussian noise (AWGN), and with the parameters \(A=1.5,\omega _0=0.1\pi ,\theta =-\pi /4\) and \(\sigma _\epsilon ^2=1\), is presented in Fig. 3.3. First, periodogram peak is found (Fig. 3.3a). Than, the subsequent amplitude/phase estimation result is presented (Fig. 3.3b).

Tip: 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 [6]

\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} 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}

In practice, for short data lengths or non-Gaussian noise, these bounds provide only approximate guides to achievable performance.

3.5 Harmonic Signal Analysis

A particularly important class of signals encountered in many practical applications is the harmonic or periodic signal. 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), \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.

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. 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. The resulting SNR and noise variance estimates are similar to those described in the previous sections.

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 \).

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) 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.

The example is the sampled current of a switch-mode power supply in a 50Hz network sampled at a 50kHz frequency [2]. Figure 3.4a shows a reconstruction of the signal with \(M=250\) harmonics. The estimated amplitudes \(\hat {A}_m\) are shown (Fig. 3.4b) as a function of the harmonic index \(m\), including the DC term at \(m = 0\). A larger magnitude indicates a more prominent harmonic component. The first non-DC harmonic amplitude \(m=1\) corresponds to the fundamental frequency, \(\omega _0\), while higher indices capture additional harmonics in the signal. The estimated fundamental frequency is 50.104Hz with the corresponding THD of about 1.6 ratio (\(\approx 160\%\)). Figure 3.4c shows estimated SNR (top) and the noise standard deviation (bottom) vary as the number of harmonics \(M\) in the model increases.

Tip: The frequency estimator is an effective ML estimator with known analytical CRLB [3].

3.6 Discrete Fourier Transform (DFT)

The discrete Fourier transform (DFT) can be viewed as a systematic way of decomposing a finite-length signal into a sum of harmonically related sinusoids. In fact, it is a special case of the harmonic signal representation discussed earlier. Specifically, setting the fundamental angular frequency to \(\omega _0 = \frac {2\pi }{N}\) and using \(N\ge L\) harmonics, the harmonic model reduces exactly to a DFT decomposition that provides a natural harmonic decomposition of the signal into \(N\) harmonics that are evenly spaced in frequency.

For simplicity, the \(L=N\) relation is used in the discussion below by zero padding the signal \(y[n]\) up to length \(N\).

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

\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}

When \(N \geq L\), the DFT allows for perfect reconstruction of the signal using its harmonic representation:

\[ \by = \bX \hat {\bw }. \]

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.

3.6.1 Single frequency analysis

Consider a signal \(y[n]\) assume a discrete frequency \(\omega _0 = \frac {2\pi }{N}k\) is given. To estimate amplitude and phase at this predefined frequency, we can form a matrix \(\bX \),

\[ \bX _1 = \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 \(\bX _1^T \bX _1\), we find that the sine and cosine columns form an orthogonal basis for this single frequency, with

\begin{align} \bx _c^T \bx _c &= \frac {N}{2}, \\ \bx _s^T \bx _s &= \frac {N}{2}, \\ \bx _c^T \bx _s &= 0. \end{align} 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 may be found in Examples 4.2 and 8.5 in [6].

Moreover, since \(\left (\bX ^T\bX \right )^{-1}\) takes a particularly simple diagonal form, 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}

Tip 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. (3.45). 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. (3.30).

3.6.2 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. (3.31)), 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.

DFT is energy conservation transform (Parseval’s Theorem) that states the relation

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

3.6.3 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. (3.27)), the resulting spectral distribution can differ substantially from the discrete, single-peak ideal (Fig. 3.5). This difference arises because, in general, \(\bX (\omega )^T\bX (\omega )\) is not orthogonal as in Eq. (3.44). 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.”

3.7 Summary

The summary of the presented approach is shown in Table 3.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.

Although the estimators presented above have been extensively analyzed for the specific case of additive white Gaussian noise (AWGN) in the statistical signal processing literature [6, 5], conducting such an analysis requires a significantly more extensive mathematical framework. Furthermore, it is worth noting that any bias and variance in these estimators can be readily approximated via Monte Carlo simulations under various parameter settings and noise distributions.

Appendices

3.A Single frequency analysis

3.A.1 Theory

• • \(\bX ^T\bX \) analysis:

  \(\seteqnumber{0}{3.}{50}\)

  \begin{equation} \bX ^T\bX = \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}

  with the following values

  \(\seteqnumber{0}{3.}{51}\)

  \begin{equation} \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}

• • \(\left (\bX ^T\bX \right )^{-1}\) analysis: The resulting matrix

  \(\seteqnumber{0}{3.}{52}\)

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

• • Finally, \(\left (\bX ^T\bX \right )^{-1}\bX ^T\by \) is

  \(\seteqnumber{1}{3.54}{0}\)

  \begin{align} w_{c,k} &=\frac {2}{N}\sum _{n=0}^{L-1} y[n]\cos (2\pi \frac {k}{N}n)\\ w_{s,k} &=\frac {2}{N}\sum _{n=0}^{L-1} y[n]\sin (2\pi \frac {k}{N}n) \end{align}

The orthogonality in more general form is given by

\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 {k}{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}

3.1.2 Power

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

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

has the \(\omega _0\)-dependent power,

\begin{equation} \begin{aligned} P_\by &= \frac {A^2}{L}\sum _{n=0}^{L-1}\cos [2](\omega _0 n)\\ &=\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}

that results from the time-limited origin of the signal \(y[n]\). For the infinite number of samples, the resulting power converges to a continuous-time power expression,

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