0
$\begingroup$

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! enter image description here 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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – MBaz Mar 21 '18 at 17:18
  • $\begingroup$ 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. $\endgroup$ – Marcus Müller Mar 21 '18 at 17:30
  • $\begingroup$ Try googling power-law detector $\endgroup$ – Stanley Pawlukiewicz Mar 21 '18 at 17:34
  • $\begingroup$ 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. $\endgroup$ – DbergOber Mar 24 '18 at 13:53
4
$\begingroup$

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 }
$\endgroup$
  • 1
    $\begingroup$ Nice answer. "When estimating parameters in some noisy observation, be sure to estimate the right parameters" :) $\endgroup$ – Marcus Müller Mar 21 '18 at 17:34
  • 1
    $\begingroup$ AKA "all models are wrong, some models are useful". :D $\endgroup$ – Peter K. Mar 22 '18 at 15:24
  • $\begingroup$ 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. $\endgroup$ – DbergOber Mar 24 '18 at 13:59
  • $\begingroup$ @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. $\endgroup$ – Cedron Dawg Mar 24 '18 at 14:21

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.