Overcoming the negative instantaneous frequencies from Hilbert transform

Question

how to avoid negative frequencies that can be obtained from instantaneous frequency estimation using Hilbert transform?

Here is what I am doing:

1. compute analytic signal, X = hilbert(x);
2. from analytic signal, unwrap the instantaneous phase
3. calculate instantaneous frequency from derivation (np.diff) of instantaneous phase

The problem I have is that the instantaneous frequency can contains negative frequencies (e.g. chirp signal).

This issue is also well descriped her:

Negative instantaneous frequency with hilbert transform using scipy hilbert

Hilbert Huang Transform: Negative value in instantaneous frequency

The best solution seems to be descriped here:

Overcoming the negative frequencies - Instantaneous frequency and amplitude estimation using Osculating Circle method

and here:

Instantaneous frequency estimation using Osculating Circle Method

An other matlab code snipped is posted here, but it has no results:

InstantFrequencyOCM​ethod

The question is, how to calculate the velocity vector of the particle and the Osculating Circle method from the analytic signal (in matlab or python)?

Thanks, Tobias

Answer 1

I tried to reproduce the paper you attached to get the corrected Instantaneous Phase and Frequency. So as far as I understood, the velocity vector is the derivation of the analytical signal, so here is the step I tried:

1. Compute the analytical signal,

X=scipy.signal.hilbert(x)

2. Compute its derivation to get the velocity vector, here I used second order derivatite, as the first order still gives me negative frequency in the final result.

def second_derivative(amplitude, sampling_rate):

d_xt=[]
for i in range (len(amplitude)-2):
    d_xt.append((amplitude[i+2]-amplitude[i])/2)
return(d_xt)

dX= second_derivative(X, fs)

note:Based on the type of your data, maybe you want to use different type of derivation to maybe deal with the noisy gradient.

3. Compute the corrected angle, corrected_phase= np.unwrap(np.angle(dX)-0.5*np.pi)

4. Compute the instantaneous frequency freq= np.diff(corrected_phase)

Answer 2

Don’t unwrap phase as bounded, but unwrap it as always increasing (or delta(ph)>=0).