Separating/ Filtering data of multiple waveforms

[wiseac]

Hello,

I have a waveform(s) that i’m trying to separate out into individual max-fall-min for each set. This is what it looks like. (Without time added)

The data comes in a one dimensional array (Voltage) then i add a time range to create a dataframe. I also add a “Row column” that is just the range 1:length(df) so that i can consistently filter data based on that row value instead of the changing df length. (If this is a bad idea please tell me so.)

The start and stops of each rise and fall fluctuate so i cant just say this one starts at 1000 then the other starts at 2000.

I want to get them separated out into a matrix where A is the first max-fall-min, then B is the next and so on. If possible id also like to exclude the rise of the waveform and just have the max-fall-min.

Thanks,
W

[rafael.guerra]

Could you annotate your picture or explain what do you mean by “max-fall-min”?

[wiseac]

Yes, Sorry im bad with lingo

I hope that helps with some understanding.

Best,
W

[wiseac]

This is also what my DataFrame looks like.

(Dont mind the corrected Voltage, im not fully calculating it in this example)

[rafael.guerra]

So you would like to filter out the high-frequency noise and estimate the times when the waveforms start to fall and rise, so that you can split the waveform in individual segments corresponding to each “period”?

[Jake]

I am still not clear on what you want.

If you look here, you will find some idea of how to find the tach crossing times. You may need to change it to find both positive and negative crossing locations.

If that is what you want then you are done.

If you want the rise time and the fall time, then once you have these locations, you can just go backwards one point at a time until the slope changes and then forwards one point at a time until the slope changes. Then you will have the rise and fall times.

With the signal oscillating at the maximum and minimum values, it is not clear what you want for these values, but once you have the tach crossing locations you can figure that out.

[wiseac]

Id like to separate the cycles basically. So each Max-Fall-Min of those lets say 5 waveform i want to separate out and leave out the rise. It would be basically the same as separating them out by their periods but I dont want the rise of each period.

If it helps this is being used for Cavity Ring Down.

(Not my image https://www.researchgate.net/figure/Typical-cavity-ring-down-signals-The-fitting-residuals-are-shown-in-the-lower-panel_fig1_335399450)

[Jake]

The approach I suggested should work with some modifications. Looking somewhat closer there is some difficulty in finding the locations automatically because of signal jitter. I don’t think you want to filter the data before your exponential fit because of potential phase changes. If you do use the DSP.jl library and the filtfilt command.

I would find the fall location by using the tach approach with a threshold set to about 0.95*(max-min)+min and specifying a negative slope. Then you can go backwards until the slope changes.

Obtain the end of signal with a similar process with a positive slope and maybe a threshold of 0.1*(max-min)+min.

Once you are happy that you have isolated the signal of interest then do your exponential fit.

In Matlab I tend to use ginput() and pick these points off using a mouse if there are only a few. I have seen a similar approach discussed on discourse for Julia. By recollection it was using Makie but I don’t have the reference. However I think there is a reasonable chance your data should work with an automated approach as I described.

[wiseac]

Hello Jake,

I will surely try it (im just not familiar with using what you describe as the tach approach but it doesnt seem to be hard). Because im not familiar with it could it also be used to separate the “pulses” (Another lingo term). We want them separated into matrices so we can perform an exponential fit on each one. (this data is coming from an oscilloscope so its not a static signal, there is small differences in each “fall” which we care about. We then will average those exponential fits together at the end but thats the easy part.

[sairus7]

Why not just a threshold detector, like vec .> 0.25 ?

Also you can use theshold detection on preprocessed data, like derivative over neighbour points y[i] - y[i-1], over a window y[i] - y[i-w], or make you own differentiator filter, for example I like this one.

Or even make you own sliding window algorithm that maximizes some criteria between left and right parts of the window: y1 = y[i-w:i], y2 = y[i:i+w], f(y1, y2) \rightarrow max

[rafael.guerra]

If the pulses are periodic, another idea would be to fit a model to the noisy data and invert for the exponential time constants in one go.

The example below fits 6 parameters: period, time-shift, max and min voltages, and the exponentials raise and fall time constants:

LsqFit code

# 1 - GENERATE INPUT SIGNAL WITH NOISE

function expsignal(t, T, t0, V1, V2, a, b)
    raise(t,t0,V1,V2) = (t<t0 || t>=t0+T/2) ? 0.0 : V1 + (V2-V1)*(1 - exp(-a*(t-t0)))
    fall(t,t0,V1,V2)  = (t<t0 || t>=t0+T/2) ? 0.0 : V1 + (V2-V1)*exp(-b*(t-t0))
    n = round(Int, length(t)/T)
    y = zero(t)
    for i in 0:n
        y += raise(t, t0 + 2*(i-1)*T/2, V1, V2) + fall(t, t0 + (2*i-1)*T/2, V1, V2)
    end
    return y
end

T  = 3e-3                   # Period [s]
V1 = 0.10                   # Min signal amplitude [V]
V2 = 0.45                   # Max signal amplitude [V]
σN = 0.02                   # Noise level
α = 2.0e4                   # 1/raise time constant [1/s]
β = 5.0e4                   # 1/fall time constant [1/s]
dt = 5e-7                   # sampling rate [s]
t = 0:dt:(3*T - dt)         # time interval sampled
t0 = T*rand()               # random start time [s]
N = length(t)

y = expsignal.(t, T, t0, V1, V2, α, β) + σN*(rand(N) .- 0.5)


# 2 - ANALYSE DATA AND FIT PERIODIC EXPONENTIAL SIGNAL

using FFTW, Statistics, LsqFit

ix0 = findmax(abs.(fft(y .- mean(y))[1:N÷2]))[2]
T0 = (N*dt)/(ix0-1)         # approximate main cycle period from FFT
V1 = maximum(y)
V2 = minimum(y)
α₀ = 1.0e5                    # [1/s] need educated initial guess (chose high value)
β₀ = 1.0e5                    # [1/s] need educated initial guess (chose high value)

P0 = [T0, t0, V1, V2, α₀, β₀]
@. model(t, p) = expsignal(t, p[1], p[2], p[3], p[4], p[5], p[6])
fit = curve_fit(model, t, y, P0)


# 3 - DISPLAY RESULTS

println("Period [s] = ", round(coef(fit)[1]; sigdigits=4))
println("Time-shift [s] = ", round(coef(fit)[2]; sigdigits=4))
println("Minimum Voltage [V] = ", round(coef(fit)[3]; sigdigits=3))
println("Maximum Voltage [V] = ", round(coef(fit)[4]; sigdigits=3))
println("1/Fall time constant [1/s] = ", round(coef(fit)[5]; sigdigits=3))
println("Fall time constant [s] = ", round(1/coef(fit)[5]; sigdigits=3))
β₁ = round(coef(fit)[6]; sigdigits=3)
println("1/Raise time constant [1/s] = ", β₁)
τᵦ = round(1/coef(fit)[6]; sigdigits=3)
println("Raise time constant [s] = ", τᵦ)

using Plots
yt = model(t,fit.param)
plot(t, y, c=:blues)
τᵦ = round(1/β; sigdigits=5)
plot!(t, yt, c=:darkblue, lw=0.5, title = "LsqFit: fall time constant, τ = $τᵦ [s], inverse = $β₁ [1/s]", titlefont=8)

[wiseac]

Thank you for that. I’m trying it right now. I appreciate the time it took you to contrive this.