Get phase information from bursty time series with amplitude variability and sharp edges

Question

I have time series like this one, and I am trying to get the phase information from it. If it matters why, because I need to compute the precision of the cycles and amplitude variability is introducing a lot of "noise" to the signal that I want to get rid of.

Time series:

Purpose: At the end of the day I want to make an autocorrelation of the signal to compute the Quality Factor.

Problem:

• The problem is that since the signal has a very sharp edge at the beginning of every cycle, any approach I try to do give me a noisy phase information that I can't use.

• Computing the Fourier Transform of the original series also contains the amplitude variability, which I don't want to use.

I tired to approaches, using a Wavelet transform and a Hilbert transform. From any of these to approaches I get a Complex time series from which I can compute the phase, in principle with out amplitude variability.

The problem is that because of the sharp edges of the times series both approaches are failing to get a smooth phase profile of the signal.

Hilbert Transform:

I am using the Hilbert transform of Matlab, which uses FFT. The thing is that since the edges are sharp, FFT is failing the get the correct profile of the signal. I have the same problem if I make a FFT and then the inverse transform, I get this type of profiles, in which there are maximums where there are not supposed to be, and of course the phase profile, or the cosine, it's not smooth at all.

Wavelet:

Using a Complex Gaussian Kernel, "cgau1" that is the one that work best.

In the figure: Wavelet return and some cuts of the angles.

In the figure: making some cuts for different scales of the angle I get the followings cosine of the phase profile. Which again they are not smooth or they loose cycles, depending on the scale.

Question

How can I get a smooth phase profile, or the cosine of the phase, from a time series with high amplitude variability and with sharp edges?

python or matlab approaches are useful

Thank you very much.

Answer 1

I recently responded to a question elsewhere that shared some similarity with this one, so parts of that are re-used here but you can refer to it for the general idea anyway.

Since most of the information you are trying to capture exists in the phase spectra, you don't really "care" about the amplitude information as long as it doesn't produce any distortions.

Therefore, I would like to clarify that my response at this point is focusing on removing the amplitude modulation...carefully.

The main idea is rather simple: Detect the local maxima and reshape each pulse proportionally so that the amplitude modulation vanishes but the timing (and shape) of the pulses is preserved. The following code snippet contains a lot of inline comments on how this is done. It is in Python but it's not incredibly hard to port it to another language (e.g. MATLAB (?)).

from pandas import DataFrame
from scipy.interpolate import interp1d
import pylab
import numpy

def getEnvelopeModels(aTimeSeries):   
    '''Detects peaks and returns a suitable interpolation model
    Parameters:
        aTimeSeries: An Nx2 DataFrame that contains the timings and the actual samples in respective columns (col 0, col 1)
    '''        
    #Detect peaks and mark their location
    #Get the derivative of pulse train to get rid of "flat tops"
    df = numpy.diff(aTimeSeries[1])    
    #Append first data point
    res = [(aTimeSeries[0].iloc[0],aTimeSeries[1].iloc[0])]
    #Simple threshold detector on the derivative to recover peaks and their location
    for k in xrange(0,len(df)):
        if df[k]>10:
            res.append((Q[0].iloc[k+1],Q[1].iloc[k+1]))
    #Append last data point
    res.append((aTimeSeries[0].iloc[-1],aTimeSeries[1].iloc[-1]))
    #Convert to array
    res = numpy.array(res)               
    #Fit a model, here it is a 'zero' model which basically maintains
    #the y value until the next (x,y) interval. This is how the pulse width is preserved.
    u_p = interp1d(res[:,0],res[:,1], kind = 'zero',bounds_error = False, fill_value=0.0)
    #Return the model and the peak data
    return (u_p,res)


#This "if" is required here so that Python doesn't execute
#the whole script even one attempts to "borrow" the getEnvelopeModels
#from another script with a simple 'from myFile import getEnvelopeModels'
if __name__ == "__main__":
    #Load the dataset
    #This is the one that was uploaded earlier to pastebin
    Q = DataFrame.from_csv("file.csv", index_col=False, header=None)

    #Clamp min to zero
    #Q[1] = Q[1]-min(Q[1]);

    #Grab the peaks' model
    F = getEnvelopeModels(Q)

    #Evaluate model to get vector
    u = numpy.array(map(F[0],Q[0]))

    #Workout the scaling coefficient by finding the 
    # 'weakest' pulse and scaling everything to that one
    # The weakest pulse was chosen to avoid distorting the limited
    # range of the weak pulses by "stretching" them to higher amplitudes.
    coef = numpy.min(F[1][1:,1]) / numpy.max(F[1][1:,1])

    #Modify the evaluated model to workout a scaling vector
    u = coef/(u/numpy.max(u));

    #Plot everything
    #Original signal
    pylab.plot(Q[0],Q[1]);
    pylab.xlabel('x');
    pylab.xlabel('y');
    pylab.title('Original');    
    pylab.grid(True)    
    pylab.figure();
    #Scaled signal
    pylab.plot(Q[0],Q[1]*u);
    pylab.xlabel('x');
    pylab.xlabel('y');
    pylab.title('Scaled');        
    pylab.grid(True)
    pylab.figure();
    #Scaling coefficient
    pylab.plot(Q[0],u,'r');
    pylab.xlabel('x');
    pylab.xlabel('y');
    pylab.title('Scaling Vector');        
    pylab.grid(True)
    pylab.show()

    #Apply scaling to the original DataFrame
    Q[1]=Q[1]*u;
    #Save the DataFrame
    Q.to_csv('processed.csv',sep=' ',index=False,header=None)

The plots are as follows:

Original Signal

Scaling Vector

Processed Signal

By the way, processed.csv can be downloaded from here.

There is a tiny little problem here with the signal having a floating minimum. I tried to clamp it to zero by subtracting the min from the vector (it is commented in the code at the moment) but it requires more work than this.

So, although the local maxima envelope effect was diminished there is still the local minima which might still carry some effect.

I did not want to do anything more drastic to clamp the minima to zero before checking that this would be alright (?). It would basically take the same "scaling" idea but now applied to the "lower part" of the signal.

Hope this helps.