How to use findpeaks to find the width of the peak at 1/e^2.
3 ビュー (過去 30 日間)
古いコメントを表示
I have a data and I want to find the peak value at 1/e^2. What I use is the following code.
x=field5(:,1);
y=field5(:,2);
e = exp(1)
[pks,locs,fwhm,proms] = findpeaks(y, x);
[~,i] = sort(pks,'descend');
peak_1 = pks(i);
esquared = max(peak_1)*1/e^2;
z = abs(y-esquared);
[p,loc] = findpeaks(-z);
plot(x,y,'b',x(loc),y(loc),'+r');
angle = diff(x(loc))/2
What I want to get is the width of the peak at the value of 1/e^2 (like FWHM but lower), but I get wrong results. The two points are not at the correct locations expecting them to be infront of each other.
1 件のコメント
dpb
2021 年 11 月 26 日
I use findpeaks only to locate the actual peak location(s) of interest and then interp1 to locate the intersection with the desired level if want precise intersection(*).
It would be a useful enhancement if the ' WidthReference' named parameter would take an absolute level or percentage of the peak height besides the two named references. I've always wondered why TMW didn't do that from the git-go.
(*) I know I had a utility function coded in prior occupation/when still had a day job, but I didn't seem to get it off the corporate machine along with a fair amount of other things I've occasionally missed since. But, it's not difficult to implement.
採用された回答
Mathieu NOE
2021 年 11 月 26 日
hello
this would be my suggestion
this code will find the peak value then from there you can tell at which y threshold value you want to search for crossing points
then x values difference will give you the peak width
clc
clearvars
% dummy data
n=100;
x=(0:n-1)/n;
y = 0.75*sin(5*x -0.5);
[max_val,locs] = max(y); % find the peak and its x axis location
threshold = max_val/2; % your value here
[t0_pos1,s0_pos1,t0_neg1,s0_neg1]= crossing_V7(y,x,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points
% ind => time index (samples)
% t0 => corresponding time (x) values
% s0 => corresponding function (y) values , obviously they must be equal to "threshold"
figure(1)
plot(x,y,'b',x(locs),max_val,'dr',t0_pos1,s0_pos1,'db',t0_neg1,s0_neg1,'dg','linewidth',2,'markersize',12);grid on
legend('signal 1','signal 1 peak value','signal 1 positive slope crossing points','signal 1 negative slope crossing points' );
% time difference between crossing points
time_width = - t0_pos1 + t0_neg1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t0_pos,s0_pos,t0_neg,s0_neg] = crossing_V7(S,t,level,imeth)
% [ind,t0,s0,t0close,s0close] = crossing_V6(S,t,level,imeth,slope_sign) % older format
% CROSSING find the crossings of a given level of a signal
% ind = CROSSING(S) returns an index vector ind, the signal
% S crosses zero at ind or at between ind and ind+1
% [ind,t0] = CROSSING(S,t) additionally returns a time
% vector t0 of the zero crossings of the signal S. The crossing
% times are linearly interpolated between the given times t
% [ind,t0] = CROSSING(S,t,level) returns the crossings of the
% given level instead of the zero crossings
% ind = CROSSING(S,[],level) as above but without time interpolation
% [ind,t0] = CROSSING(S,t,level,par) allows additional parameters
% par = {'none'|'linear'}.
% With interpolation turned off (par = 'none') this function always
% returns the value left of the zero (the data point thats nearest
% to the zero AND smaller than the zero crossing).
%
% check the number of input arguments
error(nargchk(1,4,nargin));
% check the time vector input for consistency
if nargin < 2 | isempty(t)
% if no time vector is given, use the index vector as time
t = 1:length(S);
elseif length(t) ~= length(S)
% if S and t are not of the same length, throw an error
error('t and S must be of identical length!');
end
% check the level input
if nargin < 3
% set standard value 0, if level is not given
level = 0;
end
% check interpolation method input
if nargin < 4
imeth = 'linear';
end
% make row vectors
t = t(:)';
S = S(:)';
% always search for zeros. So if we want the crossing of
% any other threshold value "level", we subtract it from
% the values and search for zeros.
S = S - level;
% first look for exact zeros
ind0 = find( S == 0 );
% then look for zero crossings between data points
S1 = S(1:end-1) .* S(2:end);
ind1 = find( S1 < 0 );
% bring exact zeros and "in-between" zeros together
ind = sort([ind0 ind1]);
% and pick the associated time values
t0 = t(ind);
s0 = S(ind);
if ~isempty(ind)
if strcmp(imeth,'linear')
% linear interpolation of crossing
for ii=1:length(t0)
%if abs(S(ind(ii))) >= eps(S(ind(ii))) % MATLAB V7 et +
if abs(S(ind(ii))) >= eps*abs(S(ind(ii))) % MATLAB V6 et - EPS * ABS(X)
% interpolate only when data point is not already zero
NUM = (t(ind(ii)+1) - t(ind(ii)));
DEN = (S(ind(ii)+1) - S(ind(ii)));
slope = NUM / DEN;
slope_sign(ii) = sign(slope);
t0(ii) = t0(ii) - S(ind(ii)) * slope;
s0(ii) = level;
end
end
end
% extract the positive slope crossing points
ind_pos = find(sign(slope_sign)>0);
t0_pos = t0(ind_pos);
s0_pos = s0(ind_pos);
% extract the negative slope crossing points
ind_neg = find(sign(slope_sign)<0);
t0_neg = t0(ind_neg);
s0_neg = s0(ind_neg);
else
% empty output
ind_pos = [];
t0_pos = [];
s0_pos = [];
% extract the negative slope crossing points
ind_neg = [];
t0_neg = [];
s0_neg = [];
end
end
2 件のコメント
その他の回答 (0 件)
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!