# Pulse duration detection

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.

• One possibility is to calculate the signal's correlation to a bank of pre-calculated pulses of different durations. Whenever you get a large peak value out of one of the correlators, that pulse has been detected.
– MBaz
Mar 21 '18 at 17:18
• Finding the onset might be as easy as using a threshold detector, possibly with a filter to suppress noise variance a bit. I, however, claim that you can't ever put an end to pulses like yours: your end-of-pulse marker seems totally arbitrary, and in fact, I can still "see" oscillations after that end marker. Mar 21 '18 at 17:30
• Try googling power-law detector
– user28715
Mar 21 '18 at 17:34
• Thank you all for the comments. I will follow your recommendations and check these options. @ Marcus Müller: you are right. The end-of-pulse marker is arbitrary. I just wanted to visualize what I am looking for. And indeed it is hard to define the end of a pulse. Mar 24 '18 at 13:53

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 }

• Nice answer. "When estimating parameters in some noisy observation, be sure to estimate the right parameters" :) Mar 21 '18 at 17:34
• AKA "all models are wrong, some models are useful". :D
– Peter K.
Mar 22 '18 at 15:24
• Hi Ced, thank you very much for your detailed answer. It looks very promising. I will definitely try it out and let you know how it works. Thank you. Mar 24 '18 at 13:59
• @DbergOber, You're welcome. This approach is basically a rough cut. If you need higher precision values send me an email (address on my profile page) and I will share some different techniques with you. Mar 24 '18 at 14:21