DIADEM script to determine rate of change of a signal based on local peaks and valley

[ebenellis]

Is it possible to determine the instantaneous rate of change of a signal between local peaks & valleys using DIADEM script. I have a table with below columns where signal and timestamp are known. The remaining columns need to be calculated as below:

i   --> index calculated

t   --> timestamp of valley or peak

x  --> peak /valley point

dx --> delta between peak and valley

dt --> delta between timestamp pf the corresponding peak and valley

Signal

timestamp

i

t

x

dx

dt

dx/dt

Below is the sample python code which implements the same logic by loading the test.csv file (attached).
How can we implement this logic when similar data is loaded in DIADEM which is currently implemented in python pandas using df.iloc function

import pandas

path = r"C:\Users\Downloads\test.csv"
filename = "post_process.csv"
path_write = r"C:\Users\Downloads\\" + filename
print(path_write)
fwrite = open(path_write, "w")
fwrite.write("i,t,x \n")
with open(path, 'r') as outfile:
    df = pandas.read_csv(outfile)
    local_peak = df.iloc[0, 0]
    local_valley = df.iloc[0, 0]
    i = 0
    fwrite.write("," + str(df.iloc[0, 1]) + "," + str(df.iloc[0, 0]) + "," + "\n")
    while i + 2 < len(df):
        if df.iloc[i, 0] < df.iloc[i + 1, 0] and df.iloc[i + 1, 0] > df.iloc[i + 2, 0]:
            fwrite.write(str(i + 1) + "," + str(df.iloc[i + 1, 1]) + "," + str(df.iloc[i + 1, 0]) + "," + "\n")
        if df.iloc[i, 0] > df.iloc[i + 1, 0] and df.iloc[i + 1, 0] < df.iloc[i + 2, 0]:
            fwrite.write(str(i + 1) + "," + str(df.iloc[i + 1, 1]) + "," + str(df.iloc[i + 1, 0]) + "," + "\n")
        i = i + 1
fwrite.close()
filename_dxdt = "post_process_dxdt.csv"
path_write_dxdt = r"C:\Users\Downloads\test_dxdt\\" + filename_dxdt
fwrite_dxdt = open(path_write_dxdt, "w")
fwrite_dxdt.write("dx,dt,dx/dt \n")
with open(path_write, 'r') as outfile:
    df = pandas.read_csv(outfile)
    print(df.iloc[0, 0], df.iloc[1, 0])
    if df.iloc[0, 1] < df.iloc[1, 1]:
        i = 0
    else:
        i = 1
    while i + 1 < len(df):
        fwrite_dxdt.write(
            str(df.iloc[i + 1, 1] - df.iloc[i, 1]) + "," + str(df.iloc[i + 1, 0] - df.iloc[i, 0]) + "," + str(
                (df.iloc[i + 1, 1] - df.iloc[i, 1]) / (df.iloc[i + 1, 0] - df.iloc[i, 0])) + "\n")
        i = i + 2
fwrite_dxdt.close()

 Attaching the sample output expected.

[Walter_Rick]

Hi joshilpa,

Please have a look at this script:

# prepare source groups and channels
oGroups       = dd.Data.Root.ChannelGroups
oSourceX      = oGroups.Item(1).Channels("TimeStamp")
oSourceY      = oGroups.Item(1).Channels("Signal")

# prepare result groups and channels
oResGroupChns = oGroups.Add("Results").Channels
oGroups.Item("Results").Activate()

oPeakX        = oResGroupChns.Add("t_peak",dd.DataTypeChnFloat64)
oPeakY        = oResGroupChns.Add("x_peak",dd.DataTypeChnFloat64)
oValX         = oResGroupChns.Add("t_valley",dd.DataTypeChnFloat64)
oValY         = oResGroupChns.Add("x_valley",dd.DataTypeChnFloat64)
odt           = oResGroupChns.Add("dt",dd.DataTypeChnFloat64)
odX           = oResGroupChns.Add("dx",dd.DataTypeChnFloat64)

# calculation
dd.ChnPeakFind(oSourceX, oSourceY, oPeakX, oPeakY, 65535, "Max.Peaks", "Time")
dd.ChnPeakFind(oSourceX, oSourceY, oValX, oValY, 65535, "Min.Peaks", "Time")
dd.ChnSub(oPeakX, oValX, odt)
dd.ChnSub(oPeakY, oValY, odX)

Greetings

Walter

[ebenellis]

Hi Walter,

This method will not work out because the no of peaks here is not fixed. The data here is huge captured at 4000000 samples per second.
The below script implements the expected logic but it takes very long time to execute because of large data size. Could you please suggest a way to optimize the below script?

# --------------------------------------------------------------------
# -- Python Script File
# -- Created on 07/18/2022 09:33:24
# -- Author: joshilpa
# -- Comment: Script to determine slew rate
# --------------------------------------------------------------------
import sys

if 'DIAdem' in sys.modules:
    from DIAdem import Application as dd

    if dd.AppEnableScriptDebugger:
        import debugpy

        debugpy.configure(python=sys.prefix + '\\python.exe')
        if not debugpy.is_client_connected():
            try:
                debugpy.listen(5678)
            except:
                pass
            debugpy.wait_for_client()

# --------------------------------------------------------------------
# -- Beginning of user code --
# prepare source groups and channels
oGroups = dd.Data.Root.ChannelGroups
oSourceX = oGroups.Item(1).Channels("timeStamp")
oSourceY = oGroups.Item(1).Channels("Signal")
# prepare result groups and channels
oResGroupChns = oGroups.Add("Results").Channels
oGroups.Item("Results").Activate()
# gives the sample number
oCount = oResGroupChns.Add("count", dd.DataTypeChnFloat64)
oPeakX = oResGroupChns.Add("t_peak", dd.DataTypeChnFloat64)
oPeakY = oResGroupChns.Add("signal_peak", dd.DataTypeChnFloat64)
odY = oResGroupChns.Add("dY", dd.DataTypeChnFloat64)
odX = oResGroupChns.Add("dX", dd.DataTypeChnFloat64)
odY_dX = oResGroupChns.Add("dY_dX", dd.DataTypeChnFloat64)
# count Signal samples
chn_len = oSourceY.Size
# value in row 1 cf column x
# print(oSourceY.Values[1])

local_peak = oSourceY.Values[1]
local_valley = oSourceY.Values[1]
i = 1
j = 1
while i + 2 < chn_len:
    # print(oSourceY.Values[i + 2])
    if oSourceY.Values[i + 2] is not None:
        if oSourceY.Values[i] < oSourceY.Values[i + 1] and oSourceY.Values[i + 1] > oSourceY.Values[i + 2]:
            oCount.Values[j] = i + 1
            oPeakX.Values[j] = oSourceX.Values[i + 1]
            oPeakY.Values[j] = oSourceY.Values[i + 1]
            j = j + 1
        if oSourceY.Values[i] > oSourceY.Values[i + 1] and oSourceY.Values[i + 1] < oSourceY.Values[i + 2]:
            oCount.Values[j] = i + 1
            oPeakX.Values[j] = oSourceX.Values[i + 1]
            oPeakY.Values[j] = oSourceY.Values[i + 1]
            j = j + 1
        if oSourceY.Values[i] == oSourceY.Values[i + 1]:
            if oSourceY.Values[i - 1] < oSourceY.Values[i + 1] and oSourceY.Values[i + 1] > oSourceY.Values[i + 2]:
                oCount.Values[j] = i + 1
                oPeakX.Values[j] = oSourceX.Values[i + 1]
                oPeakY.Values[j] = oSourceY.Values[i + 1]
                j = j + 1
            if oSourceY.Values[i - 1] > oSourceY.Values[i + 1] and oSourceY.Values[i + 1] < oSourceY.Values[i + 2]:
                oCount.Values[j] = i + 1
                oPeakX.Values[j] = oSourceX.Values[i + 1]
                oPeakY.Values[j] = oSourceY.Values[i + 1]
                j = j + 1
            # i = i + 1

    i = i + 1

# count peak signal samples
peak_chn_len = oPeakY.Size
k = 1
if oPeakY.Values[1] < oPeakY.Values[2]:
    a = 1
else:
    a = 2
print("a ", a)
while a + 1 < peak_chn_len:
    odY.Values[k] = oPeakY.Values[a + 1] - oPeakY.Values[a]
    odX.Values[k] = oPeakX.Values[a + 1] - oPeakX.Values[a]
    odY_dX.Values[k] = odY.Values[k] / odX.Values[k]
    a = a + 2
    k = k + 1

[Walter_Rick]

Hi joshilpa,

I did a performance test with the data you have provided. Here are the results:

Your script above took: 51.5s

Your script optimized: 13.4 s

My script above: 0.05s

The optimization was to remove all “.Values” because this is the default and not necessary. Example:

Old:

oPeakX.Values[j] = oSourceX.Values[i + 1]

New:

oPeakX[j] = oSourceX[i + 1]

Greetings

Walter

[ebenellis]

Hi Walter,

The optimization you suggested in my script helps.

When I try your script I don't get the same results as my script, could you please suggest how can I modify your script to get similar results to that of my script?

[Walter_Rick]

Hi joshhipa,

The difference is if two identical values follow each other. In general, I find the first value and you the second one. But in some cases, you don’t find such a value.

Part of the signal

Part of the signal with your results:

Part of the signal with my results:

What you can take into consideration is to use my script to calculate your channel section by section and attach the results.

Greetings

Walter

[ebenellis]

Hi Walter,

There is an issue with your script where you specify the 65535 in below script:
I originally have 143700000 samples in my data how do I determine the peaks and valley in that case?

dd.ChnPeakFind(oSourceX, oSourceY, oPeakX, oPeakY, 65535, "Max.Peaks", "Time")

[Walter_Rick]

Hi joshilpa,

Yes, that is a current limitation of this function. I already informed R&D about it, and they plan to remove the limit for one of the next versions.

For this reason, I mentioned to analyze your data section by section. Copy the first 60000 values into new channels, calculate it, take the next 60000 values, calculate it and append the results, etc.

I have not tested this but assume it will be faster than your script.

Greetings

Walter

[ebenellis]

Hi Walter,

Could you please suggest a way we can split 143700000 samples into smaller parts using script and move it into a sperate channel to calculate the peaks and valley followed by appending the results to an existing channel using script.
A sample script with smaller sample set will also help resolve this issue

[Walter_Rick]

Hi joshilpa,

With the following script it took 70s for a channel pair with 10000000 values:

# --------------------------------------------------------------------
def CalcPeakValley(oSourceX, oSourceY, oPeakX, oPeakY, oValX, oValY, odt, odX):
    # calculation
    dd.ChnPeakFind(oSourceX, oSourceY, oPeakX, oPeakY, 65535, "Max.Peaks", "Time")
    dd.ChnPeakFind(oSourceX, oSourceY, oValX, oValY, 65535, "Min.Peaks", "Time")
    dd.ChnSub(oPeakX, oValX, odt)
    dd.ChnSub(oPeakY, oValY, odX)
    
# --------------------------------------------------------------------
# --------------------------------------------------------------------
# --------------------------------------------------------------------

dd.StopWatchReset(1)

# prepare source groups and channels
oGroups           = dd.Data.Root.ChannelGroups
oSourceX          = oGroups.Item(1).Channels("TimeStamp")
oSourceY          = oGroups.Item(1).Channels("Signal")
# prepare result groups and channels
oResGroupChns     = oGroups.Add("Results").Channels
# prepare temporary result groups and channels
oResTempGroupChns = oGroups.Add("ResTemp").Channels
oTempSourceX      = oResTempGroupChns.Add("TimeTemp", dd.DataTypeChnFloat64)
oTempSourceY      = oResTempGroupChns.Add("SignalTemp", dd.DataTypeChnFloat64)
oPeakX            = oResGroupChns.Add("t_peak",dd.DataTypeChnFloat64)
oPeakY            = oResGroupChns.Add("x_peak",dd.DataTypeChnFloat64)
oValX             = oResGroupChns.Add("t_valley",dd.DataTypeChnFloat64)
oValY             = oResGroupChns.Add("x_valley",dd.DataTypeChnFloat64)
odt               = oResGroupChns.Add("dt",dd.DataTypeChnFloat64)
odX               = oResGroupChns.Add("dx",dd.DataTypeChnFloat64)
oPeakX_temp       = oResTempGroupChns.Add("t_peak",dd.DataTypeChnFloat64)
oPeakY_temp       = oResTempGroupChns.Add("x_peak",dd.DataTypeChnFloat64)
oValX_temp        = oResTempGroupChns.Add("t_valley",dd.DataTypeChnFloat64)
oValY_temp        = oResTempGroupChns.Add("x_valley",dd.DataTypeChnFloat64)
odt_temp          = oResTempGroupChns.Add("dt",dd.DataTypeChnFloat64)
odX_temp          = oResTempGroupChns.Add("dx",dd.DataTypeChnFloat64)

# calculate base parameter
iChnLenForCalc    = 65500
iNoOfIterations   = int(oSourceX.Size / iChnLenForCalc)
iIterationRest    = oSourceX.Size - (iChnLenForCalc * iNoOfIterations)

# loop for chunk calculation
for iLoop in range(1, iNoOfIterations+1):
    iStart = (iLoop - 1) * iChnLenForCalc + 1
    
    if iLoop == 1:
        iNoOfFormerChnk = 0
    else:
        iNoOfFormerChnk = 3
    
    oTempSourceX.RemoveValues()
    oTempSourceY.RemoveValues()
    oPeakY_temp.RemoveValues()
    oValX_temp.RemoveValues()
    oValY_temp.RemoveValues()
    odt_temp.RemoveValues()
    odX_temp.RemoveValues()
    
    oTempSourceX.SetValuesBlock(oSourceX.GetValuesBlock(iStart - iNoOfFormerChnk, iChnLenForCalc+ iNoOfFormerChnk))
    oTempSourceY.SetValuesBlock(oSourceY.GetValuesBlock(iStart - iNoOfFormerChnk, iChnLenForCalc + iNoOfFormerChnk))
    CalcPeakValley(oTempSourceX, oTempSourceY, oPeakX_temp, oPeakY_temp, oValX_temp, oValY_temp, odt_temp, odX_temp)
    dd.ChnConcatenateChannels(oPeakX, oPeakX_temp, True, False, False)
    dd.ChnConcatenateChannels(oPeakY, oPeakY_temp, True, False, False)
    dd.ChnConcatenateChannels(oValX, oValX_temp, True, False, False)
    dd.ChnConcatenateChannels(oValY, oValY_temp, True, False, False)
    dd.ChnConcatenateChannels(odt, odt_temp, True, False, False)
    dd.ChnConcatenateChannels(odX, odX_temp, True, False, False)
    
# calculate rest of ieteration
iStart = iNoOfIterations * iChnLenForCalc
oTempSourceX.RemoveValues()
oTempSourceY.RemoveValues()
oPeakY_temp.RemoveValues()
oValX_temp.RemoveValues()
oValY_temp.RemoveValues()
odt_temp.RemoveValues()
odX_temp.RemoveValues()
oTempSourceX.SetValuesBlock(oSourceX.GetValuesBlock(iStart - iNoOfFormerChnk, iIterationRest + iNoOfFormerChnk))
oTempSourceY.SetValuesBlock(oSourceY.GetValuesBlock(iStart - iNoOfFormerChnk, iIterationRest + iNoOfFormerChnk))
CalcPeakValley(oTempSourceX, oTempSourceY, oPeakX_temp, oPeakY_temp, oValX_temp, oValY_temp, odt_temp, odX_temp)
dd.ChnConcatenateChannels(oPeakX, oPeakX_temp, True, False, False)
dd.ChnConcatenateChannels(oPeakY, oPeakY_temp, True, False, False)
dd.ChnConcatenateChannels(oValX, oValX_temp, True, False, False)
dd.ChnConcatenateChannels(oValY, oValY_temp, True, False, False)
dd.ChnConcatenateChannels(odt, odt_temp, True, False, False)
dd.ChnConcatenateChannels(odX, odX_temp, True, False, False)

print(dd.StopWatch(1))

Greetings

Walter