How to use findpeaks to find the width of the peak at 1/e^2.
조회 수: 4 (최근 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!