1
$\begingroup$

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:

original 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.

enter image description here

$\endgroup$
  • 2
    $\begingroup$ 1: Is your data properly sampled (Shannon-Nyquist)? 2: Are you doing phase unwrapping? $\endgroup$ – Emanuel Landeholm Dec 12 '15 at 16:54
  • $\begingroup$ By "precision of cycles" do you mean whether or not the time interval between peaks is constant? If yes, there is a much simpler way to quantify it. Is it possible to share a little bit more information about the application? What does the signal represent? Why do you need the phase spectrum? And also, 'phase' with respect to what? (A reference square wave perhaps? ) $\endgroup$ – A_A Dec 13 '15 at 11:47
  • $\begingroup$ @EmanuelLandeholm Thanks for the answer. 1: Yes, it is satisfyng Nyquist theorem. 2: No, I am not doing unwrapping. $\endgroup$ – Iván Dec 13 '15 at 18:45
  • $\begingroup$ @A_A Sorry if my question wasn't precise enough. With precision I want to see how much a cycle differs from the other ones, and not only the interpeak interval. In particular what I am interested in is, once I got the phase of the cycle I want to compute the autocorrelation of the cosine of the phase. If all the cycles are equal, I would get a cosine, if the cycles are not equal but are similar, I will get a dumped cosine, from which I get a lot of information. $\endgroup$ – Iván Dec 13 '15 at 18:51
  • $\begingroup$ @A_A The signal is the stochastic simulation of the expression of a gene. I don't need particular the phase spectrum. What I ideally want is from obtain the amplitude and the phase as a function of time. From the phase(t) I can compute everything I need. $\endgroup$ – Iván Dec 13 '15 at 18:52
1
$\begingroup$

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 Original

Scaling Vector Scaling Coeficient

Processed Signal Without Amplitude

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.

$\endgroup$
  • $\begingroup$ Yes, it does help indeed. It's not what I am looking for absolutely, still need to compute the phase, but at least I get rid of the amplitude fluctuation. Thank you very much. Python it's ok for me also ;) $\endgroup$ – Iván Dec 14 '15 at 20:43
  • $\begingroup$ @Iván 1:If it was helpful you might want to consider upvoting it later on...Just saying :) 2: Phase plot with unwraping as suggested earlier by Emmanuel (imgbox.com/DzZRWTbD). Q[1] there is the final one, after mul with (u). I am not 100% sure "where you are going" wtih the phase though 3: I am intrigued, where can I find out more about this? (tinyurl.com/hghe59n) $\endgroup$ – A_A Dec 14 '15 at 22:26
  • $\begingroup$ I would love to but I don't have enough reputation in this community to upvote :P. Heres is the result I was looking for es.scribd.com/doc/293522011/Autocorrelation-Test $\endgroup$ – Iván Dec 17 '15 at 13:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.