how to get the envelope of oscillating signal

I have a decaying oscllatory signal for which I want to get the envelope to estimate the time constant of the system.
load data.mat
plot(t,il5)

[pks,locs]=findpeaks(il5);
hold
plot(t(locs),pks,'or')
I tried using the 'Minpeakdistance' constraint, set to 40e-6,but that still did not pick the right peaks ---in fact, I didn't see much difference.
[pks2,locs2]=findpeaks(il5,'MinPeakDistance',40e-6);
figure;
plot(t,il5)
hold
plot(t(locs2),pks2,'oc')
Anyway, how do I get the [positive] "envelope" of this signal -- the black line below?

Thanks!
Jorge

Answer 1

[Star Strider]

Perhaps this —

LD = load('data.mat')

LD = struct with fields:
    il5: [100002x1 double]
      t: [100002x1 double]

t = LD.t;
il5 = LD.il5;
Data = rmmissing([t il5]);
t = Data(:,1);
il5 = Data(:,2);

nc = find(il5 > 0, 1, 'first')-1

nc = 25001

figure
plot(t, il5)
hold on
plot(t(nc), il5(nc), '+r')
hold off
grid

tv = t(nc:end);
il5v = il5(nc:end);

[pks,locs] = findpeaks(il5v, MinPeakDistance=150);

figure
plot(tv, il5v)
hold on
plot(tv(locs), il5v(locs), '+r')
hold off
grid

% xlim([0.0049 0.006])

objfcn = @(b,x) b(1).*exp(b(2).*(x + b(3))) + b(4);
[B,rn] = fminsearch(@(b) norm(il5v(locs) - objfcn(b,tv(locs))), [2; -1E3; tv(1); 0])

B = 4×1

1.0e+03 * 0.0241 -1.2234 -0.0000 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

rn = 0.1793

figure
plot(tv, il5v)
hold on
plot(tv, objfcn(B,tv), '-r', 'LineWidth',1.5)
hold off
xlabel('Time')
ylabel('Amplitude')
grid
text(0.0075, 1.5, sprintf('$il5(t) = %.2f \\cdot e^{%.2f \\cdot (t %+.4f)} %+.4f$', B), 'Interpreter','LaTeX', 'FontSize',14)

.

[Star Strider]

Error using envelope (line 63)

Expected input to be finite.

The first element in each vector is NaN. (I checked for non-finite values initially.) That throws the envelope error.

My Answer provides the information necessary to calculate τ from the data —

EDIT — (9 Oct 2024 at 16:57)

I did not separately do a regression on the lower envelope beecause it appeears to be symmetrical with the uppeer envelope, and regreessing it separately would likely not add any useful information. The approach would be similar, with two distinct diifferences. It would be neeceessary to deeteeeeereminee the locations of the ‘valleys’ using findpeaks:

[vys,vlocs] = findpeaks(-il5v, MinPeakDistance=150);

and the objective function would be the negative of the original objective function, although that could be negated in the fminsearch call:

[B,rn] = fminsearch(@(b) norm(il5v(locs) - objfcn(b,tv(locs))), [2; -1E3; tv(1); 0])

The entire result would be —

LD = load('data.mat')

LD = struct with fields:
    il5: [100002x1 double]
      t: [100002x1 double]

t = LD.t;
il5 = LD.il5;
Data = rmmissing([t il5]);
t = Data(:,1);
il5 = Data(:,2);

nc = find(il5 > 0, 1, 'first')-1

nc = 25001

figure
plot(t, il5)
hold on
plot(t(nc), il5(nc), '+r')
hold off
grid

tv = t(nc:end);
il5v = il5(nc:end);

[vys,vlocs] = findpeaks(-il5v, MinPeakDistance=150);

figure
plot(tv, il5v)
hold on
plot(tv(vlocs), il5v(vlocs), '+r')
hold off
grid

% xlim([0.0049 0.006])

objfcn = @(b,x) b(1).*exp(b(2).*(x + b(3))) + b(4);
[BL,rn] = fminsearch(@(b) norm(il5v(vlocs) + objfcn(b,tv(vlocs))), [2; -1E3; tv(1); 0])

BL = 4×1

1.0e+03 * 0.0084 -1.1832 -0.0000 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

rn = 0.2999

figure
plot(tv, il5v)
hold on
plot(tv, -objfcn(BL,tv), '-r', 'LineWidth',1.5)
hold off
xlabel('Time')
ylabel('Amplitude')
grid
text(0.0075, 1.5, sprintf('$il5(t) = %.2f \\cdot e^{%.2f \\cdot (t %+.4f)} %+.4f$', BL), 'Interpreter','LaTeX', 'FontSize',14)

mdl = fitnlm(tv(vlocs), il5v(vlocs), @(b,t)-objfcn(b,t), BL)

Warning: The Jacobian at the solution is ill-conditioned, and some model parameters may not be estimated well (they are not identifiable). Use caution in making predictions.

mdl = 
Nonlinear regression model:
    y ~ F(b,t)

Estimated Coefficients:
           Estimate         SE          tStat       pValue 
          __________    __________    __________    _______

    b1        8.4215    5.5677e-10    1.5126e+10          0
    b2       -1183.2        6.2557       -189.15          0
    b3    -0.0038286    8.3495e-06       -458.55          0
    b4    0.00083852    0.00098477       0.85149    0.39509


Number of observations: 345, Error degrees of freedom: 342
Root Mean Squared Error: 0.0162
R-Squared: 0.998,  Adjusted R-Squared 0.998
F-statistic vs. constant model: 7.31e+04, p-value = 0

The value for the time constant from thisi regression is

I added the statistics here out of interest.

.

[Jorge]

this works great, thank you.

What are the units of the 150 in 'MinPeakDistance=150'?

How is

objfcn = @(b,x) b(1).*exp(b(2).*(x + b(3))) + b(4);
[B,rn] = fminsearch(@(b) norm(il5v(locs) - objfcn(b,tv(locs))), [2; -1E3; tv(1); 0])

doing the regression?

I had created a function to do exponential regression using a linear model:

function [m,alpha] = fitexp(xvec,yvec)
%% Curve Fit -- use exponential regression using linear model
%  xvec, yvec are matrices or row vecs
%  y= alpha*exp(B x) 
%  ln y = ln a  + Bx    map to log space to use linear regression...
%  b0=intercept, b1=slope.
% 
% y= alpha*exp(beta*x)
% ln y = ln alpha + beta*x --> use linear regression...

if(size(xvec,2)==1)
    xvec = xvec';
end
if(size(yvec,2)==1)
    yvec = yvec';
end

% organize data matrices
ymeanvec = mean(yvec,1); 
logymeanvec = log(ymeanvec)';  % take log of mean y vector
xvec_aug = [ones(length(xvec),1) xvec'];

B = mldivide(xvec_aug,logymeanvec);
intercept = B(1);
m = B(2);
alpha = exp(intercept);

but I'd like to understand how those 2 statements work....starting with the @....

excellent, thank you.

[Star Strider]

As always, my pleasure!

‘What are the units of the 150 in 'MinPeakDistance=150'?’

Those are the index values. I prefer to use index values with findpeaks (and similar functions) because I can then apply them to all necessary variables. (I left out some intermediate calculations and plots where I plotted ‘il5’ against its indices rather than against the time values to understand how the index values worked with it.)

The objective function:

objfcn = @(b,x) b(1).*exp(b(2).*(x + b(3))) + b(4);

is a generic one that generally fits exponential data. This one includes an offset (‘b(3)’) for the inidependent variable because the data do not begin at t=0. I coded it as an anonymous function. See the documentatiion on Anoymous Functions for details.

The fminsearch function is a derivative-free (that is, it does not use gradient-decent) optimisation algorithm that finds the minimum of its function argument. Here, that means subtracting the function value at every value of the independent variable from the dependent variable at every corresponding value of them. When it finds that minimum, then it reports back the parameter vector (‘B’) that corresponds to the minimum. I chose fminsearch because it is part of core MATLAB, so everyone has it.

[B,rn] = fminsearch(@(b) norm(il5v(locs) - objfcn(b,tv(locs))), [2; -1E3; tv(1); 0])

The ‘rn’ (residual norm) is the norm of the reesiduals at convergence (when the function argument value is at its minimum). The function I use calculates the norm of the residuals (essentially the square root of the sum of the squares of its argument, that being the difference between the dependent variable and the function value for both vectors) at each step.

I will gladly enlarge on this discussion of you would like more detail.

.

Answer 2

[Image Analyst]

Did youi look at envelope?

help envelope

 envelope - Signal envelope
    This MATLAB function returns the upper and lower envelopes of the input
    sequence, x, as the magnitude of its analytic signal.

    Syntax
      [yupper,ylower] = envelope(x)

      [yupper,ylower] = envelope(x,fl,'analytic')
      [yupper,ylower] = envelope(x,wl,'rms')
      [yupper,ylower] = envelope(x,np,'peak')

      envelope(___)

    Input Arguments
      x - Input sequence
        vector | matrix
      fl - Hilbert filter length
        positive integer scalar
      wl - Window length
        positive integer scalar
      np - Peak separation
        positive integer scalar

    Output Arguments
      yupper,ylower - Upper and lower signal envelopes
        vectors | matrices

    Examples
      openExample('signal/AnalyticEnvelopesOfDecayingSinusoidExample')
      openExample('signal/AnalyticEnvelopesOfDecayingSinusoidUsingFilterExample')
      openExample('signal/MovingRMSEnvelopesOfAudioRecordingExample')
      openExample('signal/PeakEnvelopesOfSpeechSignalExample')
      openExample('signal/EnvelopeOfAsymmetricSequenceExample')

    See also findpeaks, hilbert, movmax, movmean, movmin, rms

    Introduced in Signal Processing Toolbox in R2015b
    Documentation for envelope
       doc envelope

[Jorge]

Thanks, but I get the error shown below.....Not sure why the vector of length =100002 is regarded as infinite...

LD = load('data.mat')

LD = struct with fields:
    il5: [100002x1 double]
      t: [100002x1 double]

t= LD.t;
il5 = LD.il5;

[yupper,ylower] = envelope(il5)

Error using envelope (line 63)
Expected input to be finite.

[Image Analyst]

Did you look at your data to see if all the values make sense? You have nan's in there. Remove them.

load('jorge data.mat')
plot(t, il5, 'b.');
grid on;
nanIndexes = isnan(il5);
t(nanIndexes) = [];
il5(nanIndexes) = [];
[yupper,ylower] = envelope(il5);
hold on;
plot(t, ylower, 'r-');
plot(t, yupper, 'r-');

Zooming way in you see:

Looks like an envelope to me.

[Jorge]

thanks!