Mercurial > forge
changeset 5815:eb2fe511a6e5 octave-forge
New function: fwhm
author | hauberg |
---|---|
date | Sat, 11 Jul 2009 19:20:27 +0000 |
parents | d0e55a481460 |
children | 9bd92218b882 |
files | main/signal/INDEX main/signal/inst/fwhm.m |
diffstat | 2 files changed, 54 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- a/main/signal/INDEX Sat Jul 11 07:30:31 2009 +0000 +++ b/main/signal/INDEX Sat Jul 11 19:20:27 2009 +0000 @@ -34,6 +34,7 @@ grpdelay impz zplane + fwhm Filter conversion ac2poly poly2ac
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/signal/inst/fwhm.m Sat Jul 11 19:20:27 2009 +0000 @@ -0,0 +1,53 @@ +%% -*- texinfo -*- +%% @deftypefn {Function File} {} fwhm (@var{f}) +%% @deftypefnx{Function File} {} fwhm (@var{x}, @var{f}) +%% Return the full width at half maximum if a given signal. +%% +%% The full width at half maximum is an expression of the extent of a signal, +%% defined as the difference between the two extremal points of the signal where +%% it equals half of its maximum value. +%% +%% The signal is given by the vector @var{f} with optional @math{x}-values +%% @var{x}. If the latter is not given, it defaults to the indexes of @var{f}. +%% +%% If the full width at half maximum is not defined, the function returns 0. +%% @seealso{std} +%% @end deftypefn + +%% This program is public domain. +%% Author: Christophe Cossou + +function retval = fwhm (x, f) + %% Check input + if (nargin == 0) + error ('fwhm: not enough input arguments'); + elseif (nargin == 1) + f = x; + x = 1:length (f); + end + + if (numel (x) != numel (f)) + error ('fwhm: number of elements does not match'); + end + + %% Locate maximum + fmax = max (f); + f_renorm = f - 0.5 * fmax; + + ind = find (f_renorm(1:end-1) .* f_renorm (2:end) <= 0); + + %% If the product is negative, this means that the 'half-maximum' is between + %% theses two values + if (length (ind) == 2) + %% We make a linear regression between the two values to get a more precise + %% estimation of the fwhm. + x1 = (0.5 * fmax - f (ind (1) + 1)) * (x (ind (1) + 1) - x (ind (1))) ... + / (f (ind (1) + 1) - f (ind (1))) + x (ind (1) + 1); + x2 = (0.5 * fmax - f (ind (end) + 1)) * (x (ind (end) + 1) - x (ind (end))) ... + / (f (ind (end) + 1) - f (ind (end))) + x (ind (end) + 1); + retval = x2 - x1; + else + % warning ('fwhm: full width at half maximum is not defined. FWHM is set to 0'); + retval = 0; + end +