Accurate phase calculation in sinusoidal linear regression?

Question

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?

Thanks in advance.

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

---

Answer 1

What other "modern" methods exist for accurate phase extraction?

Unless the frequencies in the signal are phase locked to your sampling clock an FFT is not a great way to determine either frequency or phase of a sinusoidal component.

In many cases a Phase Locked Loop (PLL) or Delay Lock Loop works much better. You can use an FFT to quickly determine number of rough frequencies of the components and than use a PLL to track each one. A PPL allows you to trade off accuracy vs speed to convergence and even dynamically adjust it based on changing signal conditions.

Answer 2

Do an FFTShift (rotate the data halfway) before doing an FFT for phase analysis. That will re-reference the measured phase to the center of your data, not to potential circular discontinuities at the ends. Any discontinuities near the phase reference point (due to any non-integer-periodic-in-aperture frequencies) will corrupt phase interpolation.

Do not "correct" the magnitude and phase results using linear interpolation. Instead first interpolate the complex components of the FFT result using Sinc reconstruction. Then compute the magnitude and phase of the interpolated complex IQ. Successive approximation can be used to find a magnitude maxima (for frequency estimation in a low enough adjacent noise floor).

Interpolate the I and Q results of the FFT using Sinc (or windowed Sinc) interpolation to find your magnitude maxima, instead of linear interpolation (which has a much noisier frequency response).