1 %% Author: Petr Mikulik (2009)
2 %% This program is granted to the public domain.
4 %% Compute peak full-width at half maximum (FWHM) or at another level of peak
5 %% maximum for vector or matrix data y, optionally sampled as y(x). If y is
6 %% a matrix, return FWHM for each column as a row vector.
8 %% f = fwhm({x, } y {, 'zero'|'min' {, 'rlevel', rlevel}})
9 %% f = fwhm({x, } y {, 'alevel', alevel})
13 %% f = fwhm(x, y, 'zero')
14 %% f = fwhm(x, y, 'min')
15 %% f = fwhm(x, y, 'alevel', 15.3)
16 %% f = fwhm(x, y, 'zero', 'rlevel', 0.5)
17 %% f = fwhm(x, y, 'min', 'rlevel', 0.1)
19 %% The default option 'zero' computes fwhm at half maximum, i.e. 0.5*max(y).
20 %% The option 'min' computes fwhm at the middle curve, i.e. 0.5*(min(y)+max(y)).
22 %% The option 'rlevel' computes full-width at the given relative level of peak
23 %% profile, i.e. at rlevel*max(y) or rlevel*(min(y)+max(y)), respectively.
24 %% For example, fwhm(..., 'rlevel', 0.1) computes full width at 10 % of peak
25 %% maximum with respect to zero or minimum; FWHM is equivalent to
26 %% fwhm(..., 'rlevel', 0.5).
28 %% The option 'alevel' computes full-width at the given absolute level of y.
30 %% Return 0 if FWHM does not exist (e.g. monotonous function or the function
31 %% does not cut horizontal line at rlevel*max(y) or rlevel*(max(y)+min(y)) or
32 %% alevel, respectively).
34 %% Compatibility: Octave 3.x, Matlab
36 function myfwhm = fwhm (y, varargin)
38 if nargin < 1 || nargin > 5
47 if ischar(varargin{1})
55 while k <= length(varargin)
56 if strcmp(varargin{k}, 'alevel')
59 if k > length(varargin)
60 error('option "alevel" requires an argument');
63 if ~isreal(level) || length(level) > 1
64 error('argument of "alevel" must be real number');
69 if any(strcmp(varargin{k}, {'zero', 'min'}))
73 if k > length(varargin) break; end
74 if strcmp(varargin{k}, 'rlevel')
76 if k > length(varargin)
77 error('option "rlevel" requires an argument');
80 if ~isreal(level) || length(level) > 1 || level(1) < 0 || level(:) > 1
81 error('argument of "rlevel" must be real number from 0 to 1 (it is 0.5 for fwhm)');
88 if k ~= length(varargin)+1
89 error('fwhm: extraneous option(s)');
95 if (nr == 1 && nc > 1)
96 y = y'; nr = nc; nc = 1;
100 error('dimension of input arguments do not match');
103 % Shift matrix columns so that y(+-xfwhm) = 0:
105 % case: full-width at the given absolute position
108 if strcmp(opt, 'zero')
109 % case: full-width at half maximum
110 y = y - level * repmat(max(y), nr, 1);
112 % case: full-width above background
113 y = y - level * repmat((max(y) + min(y)), nr, 1);
117 % Trial for a "vectorizing" calculation of fwhm (i.e. all
118 % columns in one shot):
119 % myfwhm = zeros(1,nc); % default: 0 for fwhm undefined
120 % ind = find (y(1:end-1, :) .* y(2:end, :) <= 0);
121 % [r1,c1] = ind2sub(size(y), ind);
122 % ... difficult to proceed further.
123 % Thus calculate fwhm for each column independently:
124 myfwhm = zeros(1,nc); % default: 0 for fwhm undefined
127 ind = find((yy(1:end-1) .* yy(2:end)) <= 0);
128 if length(ind) >= 2 && yy(ind(1)) > 0 % must start ascending
131 [mx, imax] = max(yy); % protection against constant or (almost) monotonous functions
132 if length(ind) >= 2 && imax >= ind(1) && imax <= ind(end)
135 xx1 = x(ind1) - yy(ind1) * (x(ind2) - x(ind1)) / (yy(ind2) - yy(ind1));
138 xx2 = x(ind1) - yy(ind1) * (x(ind2) - x(ind1)) / (yy(ind2) - yy(ind1));
139 myfwhm(n) = xx2 - xx1;
145 %! x=-pi:0.001:pi; y=cos(x);
146 %! assert( abs(fwhm(x, y) - 2*pi/3) < 0.01 );
149 %! assert( fwhm(-10:10) == 0 && fwhm(ones(1,50)) == 0 );
153 %! y1=-4+zeros(size(x)); y1(4:10)=8;
154 %! y2=-2+zeros(size(x)); y2(4:11)=2;
155 %! y3= 2+zeros(size(x)); y3(5:13)=10;
156 %! assert( max(abs(fwhm(x, [y1;y2;y3]') - [20.0/3,7.5,9.25])) < 0.01 );
159 %! x=1:3; y=[-1,3,-1]; assert(abs(fwhm(x,y)-0.75)<0.001 && abs(fwhm(x,y,'zero')-0.75)<0.001 && abs(fwhm(x,y,'min')-1.0)<0.001);
162 %! x=1:3; y=[-1,3,-1]; assert(abs(fwhm(x,y, 'rlevel', 0.1)-1.35)<0.001 && abs(fwhm(x,y,'zero', 'rlevel', 0.1)-1.35)<0.001 && abs(fwhm(x,y,'min', 'rlevel', 0.1)-1.40)<0.001);
165 %! x=1:3; y=[-1,3,-1]; assert(abs(fwhm(x,y, 'alevel', 2.5)-0.25)<0.001 && abs(fwhm(x,y,'alevel', -0.5)-1.75)<0.001);
168 %! x=-10:10; assert( fwhm(x.*x) == 0 );
171 %! x=-5:5; y=18-x.*x; assert( abs(fwhm(y)-6.0) < 0.001 && abs(fwhm(x,y,'zero')-6.0) < 0.001 && abs(fwhm(x,y,'min')-7.0 ) < 0.001);