2
$\begingroup$

I’m wondering what is the best algorithm for unwrapping phase of a signal whose data type is single (float) that posess the least sensitivity to roundoff error. In my application I need to compute the instantaneous frequency from the unwrapped phase. For example, just for the sake of illustration, lets consider a simple sinusoid at 200 Hz that is sampled at a rate 2000 Hz. I used the code below to compute the instantaneous unwrapped phase

fs = 2000; % sample rate f0 = 200; % carrier frequency T = 32; % duration, the higher T, the more roundoff noise is acrrued t = (single(0):1/fs:T)'; x = cos(*pi*f0*t); xhilb = hilbert(x); xangl = angle(xhilb); puwrp = unwrap(xangl);

And then from puwrp, the instantaneous frequency is computed as below

finst = (fs/(2*pi))*diff(puwrp);

The $l_\infty$ of the error f0-finst is on the order of 40 Hz.

norm(f0-finst,inf) ans = single 40.8451

Based on a simple perturbation analysis, the culprit dwells in the unwrapping algorithm which has a perturbation gain of 4096!.

theta = xangl; puwrp = unwrap(theta); % Add perturbation (numerically theta and dtheta are the same) dtheta = theta + eps(theta)/2; dpuwrp = unwrap(dtheta); norm(puwrp-dpuwrp,inf)/norm(theta-dtheta,inf) ans = single 4096 ```
Improve this question
$\endgroup$
3
  • $\begingroup$ Have you tried $f[n] = \dfrac{F_s/2}{\pi} \mathrm{atan2}\left(z[n]z^*[n-1]\right)$? It is OK for small deviations. $\endgroup$ Commented Dec 11, 2024 at 1:53
  • $\begingroup$ Have you tried least squares phase unwrapping? Like this one mathworks.com/matlabcentral/fileexchange/… $\endgroup$ Commented Dec 11, 2024 at 2:16
  • $\begingroup$ @AndyWalls: I'm not sure how this approach would improve accuracy in single precision. $\endgroup$ Commented Dec 12, 2024 at 19:30

1 Answer 1

Reset to default
3
$\begingroup$

I think that phase unwrapping is not the best solution for computing the instantaneous frequency. You can't avoid (approximative) differentiation, but it's more stable to directly use the signal and its derivatives instead of a phase function which may contain artifacts due to unwrapping.

The instantaneous frequency is defined as

$$f(t)=\frac{1}{2\pi}\phi'(t)\tag{1}$$

where $\phi(t)$ is the instantaneous phase of an analytic signal

$$w(t)=u(t)+jv(t)=r(t)e^{j\phi(t)}\tag{2}$$

Taking the derivative of $(2)$ gives

\begin{align*} w'(t) &= r'(t)e^{j\phi(t)}+jr(t)e^{j\phi(t)}\phi'(t)\\ & = w(t)\left(\frac{r'(t)}{r(t)}+j\phi'(t)\right) \tag{3} \end{align*}

from which it follows that

\begin{align*} \phi'(t) &= \textrm{Im}\left\{\frac{w'(t)}{w(t)}\right\} \\ &= \textrm{Im}\left\{\frac{u'(t)+jv'(t)}{u(t)+jv(t)}\right\} \\ &= \frac{u(t)v'(t)-u'(t)v(t)}{u^2(t)+v^2(t)} \tag{4} \end{align*}

Improve this answer
$\endgroup$
3
  • $\begingroup$ Thanks @Matt L. for your prompt reply. This is exactly what I tried before using the unwrap approach. In single precision arithimetic, the result of this approach is worse than the unwrap approach, which motivated me to take the unwrap approach. Is there any other ways to find the instantaneous frequency from the equivalent analytic form of a signal? $\endgroup$ Commented Dec 10, 2024 at 17:46
  • $\begingroup$ @AHT: I'm surprised that your results are worse. How did you compute the derivatives? There's a paper comparing several practical ways of computing the expression $(4)$ in my answer. $\endgroup$ Commented Dec 10, 2024 at 17:58
  • $\begingroup$ @ Matt L.: I dont have access to that paper. I used both forward difference and central difference approximations which lead to similar results. $\endgroup$ Commented Dec 11, 2024 at 19:58

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.