how to get the envelope of oscillating signal

8 Ansichten (letzte 30 Tage)
Jorge
Jorge am 8 Okt. 2024
Kommentiert: Jorge am 24 Okt. 2024

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

Akzeptierte Antwort

Star Strider
Star Strider am 8 Okt. 2024
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)
.
  3 Kommentare
Jorge
Jorge am 9 Okt. 2024
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
Star Strider am 9 Okt. 2024
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.
.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (1)

Image Analyst
Image Analyst am 9 Okt. 2024
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
  3 Kommentare
Image Analyst
Image Analyst am 9 Okt. 2024
Bearbeitet: Image Analyst am 10 Okt. 2024
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
Jorge am 24 Okt. 2024
thanks!

Melden Sie sich an, um zu kommentieren.

Kategorien

Mehr zu Descriptive Statistics finden Sie in Help Center und File Exchange

Tags

Produkte


Version

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by