# Accurate phase calculation in sinusoidal linear regression?

I've been trying to work out a way to minimize the error in phase calculation. The underlying model is the following

$$s(t) = \sum_{i=1}^{M} A_i\sin\left(\frac{2\pi t}{T_i} + \phi_i\right) + \epsilon_t,$$

where $$\epsilon_t\sim N(0,\sigma^2)$$. I'm assuming that the nature of an arbitrary signal (without linear trend) may be described in terms of the most relevant harmonic components; these components are to be extracted from the FFT.

Until now, I've noticed that the main hindrance in accurate phase determination is due to spectral leakage. In this regard, I've tried using the following correction algorithm (as described in [1]) assuming no windowing is performed:

\begin{align} \nu_k &= (K + \Delta K)\dfrac{\nu_s}{N} &\qquad \text{Frequency correction} \\ Y'_k &= \dfrac{\pi\Delta K Y_k}{\sin\left(\pi\Delta K\right)} &\qquad \text{Amplitude correction}\\ \phi_k &= \arctan\left(\frac{I_k}{R_k}\right) + \Delta K\pi &\qquad \text{Phase correction} \end{align}

where

$$\Delta K = \begin{cases} \frac{Y_{k+1}}{Y_{k+1}+Y_{k}} & Y_{k+1} \geq Y_{k-1}\\ \frac{-Y_{k}}{Y_{k}+Y_{k-1}} & Y_{k+1} < Y_{k-1} \end{cases},$$

and $$I_k$$ and $$R_k$$ are the imaginary and real components of the $$k$$-th element in the FFT. Additionally, I've tried using a similar correction algorithm (also described in [1]) which employs the Hann window in order to deal with spectral leakage.

The questions:

1. Aside from spectral leakage, what other factors affect phase calculation?
2. Assuming a linear regression model as $$\hat{y}(t) = \sum_{k=1}^{N/2 - 1} \left[a_k\sin\left(\frac{2\pi k}{N}t\right) + b_k\cos\left(\frac{2\pi k}{N}t\right)\right],$$ what's the difference between phase calculation in this model and using the FFT?
3. What other "modern" methods exist for accurate phase extraction?

References

1. Ming, X., & Kang, D. (1996). Corrections for frequency, amplitude and phase in a fast Fourier transform of a harmonic signal. Mechanical Systems and Signal Processing, 10(2), 211-221.

EDIT 1:

My answer to question two would be: if linear regression is solved via minimum mean square error (MMSE), then the DFT and the LR are equivalent.

Given an N-sequence $$x[n] = \{x_0,x_1,\dots,x_{N-1} \}$$, the DFT of $$x$$ is defined as

\begin{align*} X[k] & = \frac{1}{N}\sum_{j=0}^{N-1}x[j]e^{2\pi ij \frac{k}{N}}, \qquad k=0,1,\dots,N-1\\ & = \frac{1}{N}\sum_{j=0}^{N-1}\left[\cos\left(\frac{2\pi k}{N}j\right) + i\sin\left(\frac{2\pi k}{N}j\right)\right]x[j] \\ & = \frac{1}{N}\sum_{j=0}^{N-1}\left\{x[j]\cos\left(\frac{2\pi k}{N}j\right) + ix[j]\sin\left(\frac{2\pi k}{N}j\right) \right\} \\ & = \frac{1}{N}\left(a_k + ib_k\right) \end{align*}

The radius $$r_k$$ and the angle $$\theta_k$$ of $$X[k]$$ are given by

\begin{align} r_k & = \frac{1}{N}\sqrt{a_k^2 + b_k^2} \nonumber\\ & = \frac{1}{N}\sqrt{\left[\sum_{j=0}^{N-1} x[j]\cos\left(\frac{2\pi k}{N}j\right)\right]^2 + \left[\sum_{j=0}^{N-1} x[j]\sin\left(\frac{2\pi k}{N}j\right)\right]^2} \label{mangitude_DFT} \\ \theta_k & = \arctan\left(\frac{b_k}{a_k}\right) \nonumber\\ & = \arctan\left[\frac{\sum_{j=0}^{N-1}x[j]\sin\left(\frac{2\pi k}{N}j\right)}{\sum_{j=0}^{N-1}x[j]\cos\left(\frac{2\pi k}{N}j\right)}\right]\label{angle_DFT} \end{align}

On the other hand, the coefficients of LR are given by

\begin{align} a_m & = \frac{2}{N}\sum\limits_{i=1}^{N}y_i\sin\left(\frac{2\pi m}{N}i\right) \\[5pt] b_m & = \frac{2}{N}\sum\limits_{i=1}^{N}y_i\cos\left(\frac{2\pi m}{N}i\right). \end{align}

and the LR model may be equivalently expressed as

$$\widehat{y}(t) = A_k\sum_{k=1}^{N/2}\cos\left(2\pi\nu_k t - \phi_k\right)$$

where the parameters are given by

\begin{align} A_k & = \sqrt{a_k^2 + b_k^2} = \frac{2}{N}\sqrt{\left[\sum\limits_{i=1}^{N}y_i\sin\left(\frac{2\pi k}{N}i\right)\right]^2 + \left[\sum\limits_{i=1}^{N}y_i\cos\left(\frac{2\pi k}{N}i\right)\right]^2} \label{mangitude_LR}\\ \sin\phi_k & = \frac{a_k}{A_k} \nonumber\\ \cos\phi_k & = \frac{b_k}{A_k} \nonumber\\ \phi_k & = \arctan\left( \frac{a_k}{b_k} \right) = \arctan\left[ \frac{\sum\limits_{i=1}^{N}y_i\sin\left(\frac{2\pi k}{N}i\right)}{\sum\limits_{i=1}^{N}y_i\cos\left(\frac{2\pi k}{N}i\right)} \right]. \end{align}

A quick inspection reveals that both the DFT and LR are equivalent in the frequency interval $$\nu_k = \left\{1/N,\dots,\frac{N/2-1}{N}\right\}$$.

• Could you say something about the range of frequencies, sampling rate, number of samples?
– Royi
Apr 28 at 14:48
• AFAIK with option 2, you will reach the Cramer-Rao Lower Bound
– Ben
Apr 28 at 14:58
• @Royi Since I'm working with arbitrary signals I don't think I could say something specfic about the sampling characteristics. Nevertheless, I have an example of phase extraction in a dual tone model implemented in MATLAB, maybe that could be useful? Apr 28 at 15:01
• Eric Jacobsen has a few frequency estimation algorithms on his page that might be relevant.
– Peter K.
Apr 28 at 15:04
• The difficulty of estimating the parameters you're after is governed by the ration between the 2 closest tones and the observation interval (Given you're above Nyquist rate). So having those defined is important to assist you solving the problem.
– Royi
Apr 28 at 17:45