Pulse duration detection

Question

I have a signal sequence where multiple pulses with different amplitude and duration may occur at random times. White noise is present as well. I was searching for a Matlab function which could detect the pulses successfully so I can calculate their duration but couldn't find anything appropriate. I am attaching a snapshot of one pulse. Does anyone know of a function or a way to detect such pulses?

Thank you in advance! ps: One of the functions I tried was findchangepts with findchangepts(data,'MaxNumChanges',10,'MinDistance',1e3) but the function is very slow for big data and the result was not what I need. MinDistance is set to 1e3 because I expect the pusles to be at least 1e3 samples long. The one shown in the figure is 1e4 samples long.

Answer 1

That looks a lot like an exponentially decaying sinusoid. What you are primarily interested in is the decay rate. Where it ends would then have to be defined as when it reached some threshold level.

Create a subsequence consisting of the peak values and the negative values of the troughs. This should give you a nice exponential decay function.

Then follow my answer in this question: least square fitting to inverse exponential function

You may be able to ignore the $U$ component to simplify it. If you can, then there is alternative approach to find the decay rate where you simply take the log of your signal and do a standard linear regression.

Hope this helps.

Ced

---

Followup

I don't use MATLAB, but here is some Python code.

Since you have so many points, and a lot of noise, I think it would be a good idea to smooth your data some to make the peaks easier to find. The values in my program were taken from your pic and are in units of pixels.

An exponentially decaying sinusoidal is a solution of a linear differential equation in the form $ y" + a y' + b y = 0 $ where the determinant of the characteristic equation is negative.

 
import numpy as np

#=========================================================================
def main():

#---- The curve graph in pixels

        theTopCurve = [ [ 283 ,  34 ], \
                        [ 320 ,  151 ], \
                        [ 358 ,  320 ], \
                        [ 395 ,  469 ], \
                        [ 430 ,  452 ], \
                        [ 470 ,  482 ], \
                        [ 504 ,  489 ], \
                        [ 546 ,  497 ] ]

        theBotCurve = [ [ 302 ,  938 ], \
                        [ 339 ,  851 ], \
                        [ 379 ,  619 ], \
                        [ 412 ,  591 ], \
                        [ 450 ,  564 ], \
                        [ 487 ,  559 ], \
                        [ 526 ,  548 ] ]

        theMidHeight = 524

#---- Combine into one

        theCurve = []

        for i in range( 0, len( theTopCurve ) ):
            theTuple = theTopCurve[i]
            theTuple[1] = theMidHeight - theTuple[1]
            theCurve.append( theTuple )               

        for i in range( 0, len( theBotCurve ) ):
            theTuple = theBotCurve[i]
            theTuple[1] -= theMidHeight
            theCurve.append( theTuple )               

        theCurve.sort()

        print theCurve

#---- Find fit and toss outliers until good

        theTossedCount = 0
        theTossedLimit = len( theCurve ) / 4

        for a in range( 0, 10 ):
            print
            c, m, b = Iterate( theCurve )
            print c, m, b
            theTossedCount += c
            if( theTossedCount >= theTossedLimit ): break

#---- Print the Results

        print "The decay rate is ", ( "%7.5f" % m )
        print " The amplitude is ", ( "%7.5f" % np.exp( b ) )


#=========================================================================
def Iterate( argCurve ):

#---- Set relative start time

        theStartTime = argCurve[0][0]

#---- Convert into Vectors

        thePointCount = len( argCurve )

        theT = np.zeros( thePointCount )
        theY = np.zeros( thePointCount )
        theL = np.zeros( thePointCount )

        for p in range( 0, thePointCount ):
            theT[p] = argCurve[p][0] - theStartTime
            theY[p] = argCurve[p][1]
            theL[p] = np.log( theY[p] )

#---- Find the Decay Rate Parameters

        m, b = np.polyfit( theT, theL, 1 )

#---- Measure the Fit

        theR = m * theT + b

        theError  = theL - theR
        theError2 = theError * theError
        theRms    = np.sqrt( np.mean( theError2 ) )

        print theError
        print "RMS=", theRms

#---- Toss the Outliers

        theTossThreshold = theRms * 1.5
        theTossCount = 0

        for p in range( thePointCount - 1, -1, -1 ):
            if( abs( theError[p] ) > theTossThreshold ):
                print "Toss:", p, argCurve[p]
                theTossCount += 1
                argCurve.pop( p )

        return theTossCount, m, b

#=========================================================================
main()

# ln( y ) = m * t + b

# y = e^b * e^{ m * t }