how to get the envelope of oscillating signal

13 次查看（过去 30 天）
Jorge
Jorge 2024-10-8
评论： Jorge 2024-10-24，16:09

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

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2024-10-8
Ran in:
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 个评论
显示 1更早的评论
Jorge
Jorge 2024-10-9
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 2024-10-9
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.
.

请先登录，再进行评论。

更多回答（1 个）

Image Analyst
Image Analyst 2024-10-9
Ran in:
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 个评论
显示 1更早的评论
Image Analyst
Image Analyst 2024-10-9
编辑：Image Analyst 2024-10-10
Ran in:
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Descriptive Statistics 的更多信息

标签

产品


版本

R2023b

Community Treasure Hunt

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

Start Hunting!

Translated by