octave_packages/signal-1.1.3/fwhm.m

   1 %% Author: Petr Mikulik (2009)
   2 %% This program is granted to the public domain.
   3
   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.
   7 %%   Syntax:
   8 %%      f = fwhm({x, } y {, 'zero'|'min' {, 'rlevel', rlevel}})
   9 %%      f = fwhm({x, } y {, 'alevel', alevel})
  10 %%   Examples:
  11 %%      f = fwhm(y)
  12 %%      f = fwhm(x, y)
  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)
  18 %%
  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)).
  21 %%
  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).
  27 %%
  28 %% The option 'alevel' computes full-width at the given absolute level of y.
  29 %%
  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).
  33 %%
  34 %% Compatibility: Octave 3.x, Matlab
  35
  36 function myfwhm = fwhm (y, varargin)
  37
  38     if nargin < 1 || nargin > 5
  39         print_usage;
  40     end
  41     opt = 'zero';
  42     is_alevel = 0;
  43     level = 0.5;
  44     if nargin==1
  45         x = 1:length(y);
  46     else
  47         if ischar(varargin{1})
  48             x = 1:length(y);
  49             k = 1;
  50         else
  51             x = y;
  52             y = varargin{1};
  53             k = 2;
  54         end
  55         while k <= length(varargin)
  56             if strcmp(varargin{k}, 'alevel')
  57                 is_alevel = 1;
  58                 k = k+1;
  59                 if k > length(varargin)
  60                     error('option "alevel" requires an argument');
  61                 end
  62                 level = varargin{k};
  63                 if ~isreal(level) || length(level) > 1
  64                     error('argument of "alevel" must be real number');
  65                 end
  66                 k = k+1;
  67                 break
  68             end
  69             if any(strcmp(varargin{k}, {'zero', 'min'}))
  70                 opt = varargin{k};
  71                 k = k+1;
  72             end
  73             if k > length(varargin) break; end
  74             if strcmp(varargin{k}, 'rlevel')
  75                 k = k+1;
  76                 if k > length(varargin)
  77                     error('option "rlevel" requires an argument');
  78                 end
  79                 level = varargin{k};
  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)');
  82                 end
  83                 k = k+1;
  84                 break
  85             end
  86             break
  87         end
  88         if k ~= length(varargin)+1
  89             error('fwhm: extraneous option(s)');
  90         end
  91     end
  92
  93     % test the y matrix
  94     [nr, nc] = size(y);
  95     if (nr == 1 && nc > 1)
  96         y = y'; nr = nc; nc = 1;
  97     end
  98
  99     if length(x) ~= nr
 100         error('dimension of input arguments do not match');
 101     end
 102
 103     % Shift matrix columns so that y(+-xfwhm) = 0:
 104     if is_alevel
 105             % case: full-width at the given absolute position
 106             y = y - level;
 107     else
 108         if strcmp(opt, 'zero')
 109             % case: full-width at half maximum
 110             y = y - level * repmat(max(y), nr, 1);
 111         else
 112             % case: full-width above background
 113             y = y - level * repmat((max(y) + min(y)), nr, 1);
 114         end
 115     end
 116
 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
 125     for n=1:nc
 126         yy = y(:, n);
 127         ind = find((yy(1:end-1) .* yy(2:end)) <= 0);
 128         if length(ind) >= 2 && yy(ind(1)) > 0 % must start ascending
 129             ind = ind(2:end);
 130         end
 131         [mx, imax] = max(yy); % protection against constant or (almost) monotonous functions
 132         if length(ind) >= 2 && imax >= ind(1) && imax <= ind(end)
 133             ind1 = ind(1);
 134             ind2 = ind1 + 1;
 135             xx1 = x(ind1) - yy(ind1) * (x(ind2) - x(ind1)) / (yy(ind2) - yy(ind1));
 136             ind1 = ind(end);
 137             ind2 = ind1 + 1;
 138             xx2 = x(ind1) - yy(ind1) * (x(ind2) - x(ind1)) / (yy(ind2) - yy(ind1));
 139             myfwhm(n) = xx2 - xx1;
 140         end
 141     end
 142 end
 143
 144 %!test
 145 %! x=-pi:0.001:pi; y=cos(x);
 146 %! assert( abs(fwhm(x, y) - 2*pi/3) < 0.01 );
 147 %!
 148 %!test
 149 %! assert( fwhm(-10:10) == 0 && fwhm(ones(1,50)) == 0 );
 150 %!
 151 %!test
 152 %! x=-20:1:20;
 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 );
 157 %!
 158 %!test
 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);
 160 %!
 161 %!test
 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);
 163 %!
 164 %!test
 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);
 166 %!
 167 %!test
 168 %! x=-10:10; assert( fwhm(x.*x) == 0 );
 169 %!
 170 %!test
 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);