% Script that compares my "halfwidth.m" function FWHM.m to ChatGPT's
% version computeFWHM, to measure the halfwidth of individual peaks.

xo=2; % Approximate x-value of the peak whose FWHM is desired
deltax=.01; % simulated signal sampling interval
x=-4:deltax:5; % Create x axis
y=sin(x).^2; % Could cos or gaussian or any other peak function
plot(x,y,'.');
grid;

FWHM=halfwidth(x,y,xo);
disp('FWHM result = ' )
disp (FWHM)

FWHMchat=computeFWHM(x,y,xo);
disp('FWHMChatGPT result = ')
disp(FWHMchat(1))

function [FWHM,slope1,slope2,hwhm1,hwhm2]=halfwidth(x,y,xo)
% function [FWHM,slope1,slope2,x1,x2]=halfwidth(x,y,xo) computes the full
% width at half maximum of the maximum of any shape peak that peaks near xo
% and has a zero baseline. If xo is omitted, it computes the halfwidth from
% the maximum y value. Not highly accurate if the function is too sparsely
% sampled. Also optionally returns the leading and trailing edge slopes
% (slope1 and slope 2, respectively) and the leading and trailing edge
% half-width-at-half-maximum, if those output arguments are incuded.
% Tom O'Haver (toh@umd.edu) Version 2, July 2016
%
% Example 1:
% x=-5:.01:5;
% y=sinc(x).^2;
% FWHM=halfwidth(x,y,0) % FWHM of central peak
% FWHM=halfwidth(x,y,1.5) % FWHM of positive side peak
%
% Example 2:
% x=[0:.1:10];
% W=3; % W is the true half-width
% y=gaussian(x,5,W);
% FWHM=halfwidth(x,y,5)
% FWHM=halfwidth(x,y) % gives the same result since there is only one peak
%
% Example 3: % Returns NaN because signal does not contain the halfwidth
% x=[0:.1:1];
% y=exp(-(x).^2)
% halfwidth(x,y)
%
% Example 4: Also returns the leading and trailing edge slopes (+1 and -1)
% and the leading and trailing edge half-width-at-half-maximum
% x=[0:.1:10];
% y=triangle(x,5,1);
% [FWHM,slope1,slope2,x1,x2]=halfwidth(x,y)
%
% Example 5. Comparison of Gaussian (g) and BiGaussian peaks (bg)
% x=[1:100];
% g=gaussian(x,50,20);
% bg=bigaussian(x,50,20,3);
% plot(x,g,x,bg)
% [GaussianFWHM,slope1,slope2,hwhm1,hwhm2]=halfwidth(x,g)
% [BiGaussianFWHM,slope1,slope2,hwhm1,hwhm2]=halfwidth(x,bg)

if nargin==2,xo=x(val2ind(y,max(y)));end

try
    indmax=val2ind(x,xo);
    maxy=y(indmax);
    oy=y-maxy/2;
    
    n=indmax;
    while oy(n)>0
        n=n-1;
    end
    x1=interp1([y(n) y(n+1)],[x(n) x(n+1)],maxy/2);
    slope1=(y(n+1)-y(n))./(x(n+1)-x(n));
    
    n=indmax;
    while oy(n)>0
        n=n+1;
    end
    x2= interp1([y(n-1) y(n)],[x(n-1) x(n)],maxy/2);
    slope2=(y(n-1)-y(n))./(x(n-1)-x(n));
    
    FWHM=x2-x1;
    hwhm1=xo-x1;
    hwhm2=x2-xo;
catch me
    FWHM=NaN;
    slope1=NaN;
    slope2=NaN;
    hwhm1=NaN;
    hwhm2=NaN;
    me=me
end
end

function [index,closestval]=val2ind(x,val)
% Returns the index and the value of the element of vector x that is closest to val
% If more than one element is equally close, returns vectors of indicies and values
% Tom O'Haver (toh@umd.edu) October 2006
% Examples: If x=[1 2 4 3 5 9 6 4 5 3 1], then val2ind(x,6)=7 and val2ind(x,5.1)=[5 9]
% [indices values]=val2ind(x,3.3) returns indices = [4 10] and values = [3 3]
dif=abs(x-val);
index=find((dif-min(dif))==0);
closestval=x(index);
end

function fwhm = computeFWHM(x, y, xo)
% fwhm = computeFWHM(x, y, xo)
% ChatGPT request: "write a matlab function that accepts a time series x,y,
% and computes the full width at half maximum of the of the peak in y
% nearest xo If xo is omitted, it computes the halfwidth from the maximum y
% value." December, 2022.

  % Find the nearest x value to xo
  [~, i] = min(abs(x - xo));
  
  % Find the maximum y value within a certain range around the nearest x value
  ymax = max(y(max(1, i - 10):min(length(y), i + 10)));
  imax = find(y == ymax);
  
  % If xo is not provided, use the x value of the maximum y as the reference point
  if nargin < 3
    xo = x(imax);
  end
  
  % Find the half maximum value
  yhalf = ymax / 2;
  
  % Find the indices of the points where y crosses the half maximum value
  i1 = find(y(1:imax) < yhalf, 1, 'last');
  i2 = find(y(imax:end) < yhalf, 1, 'first') + imax - 1;
  
  % Calculate the full width at half maximum
  fwhm = abs(x(i2) - x(i1));
end