Get phase information from bursty time series with amplitude variability and sharp edges
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]&gt;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.