fwhm.m

Hello,
It seems that if I want to propose a new function to the community, I must send an e-mail to this list.

As the title say, I have programmed a function that return the fwhm (full width at half maximum) of an input function. I programmed it because I didn't found it in the already existings functions.

So if you want to add this function to octave, I offer you my code. Every one can modify and improve my function if it seems that my own version is not the best possible.

I don't know how to proceed and if the list of octave-dev is important, I don't think that join my code in this mail as an attachment file is a good idea.

Autiwa
------------------------------------------------------------------------------

_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[Søren Hauberg]

Hi and welcome; it's always nice to see new contributors :-)

> As the title say, I have programmed a function that return the fwhm
> (full width at half maximum) of an input function. I programmed it
> because I didn't found it in the already existings functions.

Thanks. It seems that such a function is not part of Matlab core, so it
should probably go into one of the packages. Since I do not know what
this function does, I don't know where it should possibly go. Do you (or
anybody else) have any thoughts on which package might be suitable?

> I don't know how to proceed and if the list of octave-dev is
> important, I don't think that join my code in this mail as an
> attachment file is a good idea.

Just send the code as an attachment to this mail. We only do this for
new contributors and people looking for advice, so it doesn't happen
that often. Therefore it is quite acceptable to send code to this list
(but please remember to include license information in the file).

Søren


------------------------------------------------------------------------------
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[Søren Hauberg]

tor, 02 07 2009 kl. 22:36 +0200, skrev Autiwa:
> I think that the last 3 mail (this one count) were only in private
> conversation

Ohh, I didn't notice. I'm CC'ing the list.

> Here is the program with comments in english, maybe less complete, and
> with english error, because I talk easier in french...

Thanks

> I have also added the licence.

You have added:

%## This program is public domain. (GPL)

I hate to be focused on boring license stuff, but we need this to be
correct. The GPL is different from the public domain. So, either put it
in the public domain or release it under the GPL. You have to decide
which you prefer.

Autiwa has suggested that this function goes in the 'signal' package. I
don't really have an opinion on this. What does the rest of the people
on this list think?

Søren

[fwhm.m]

---

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function that return the full width at half maximum
%## This program is public domain. (GPL)
% x: define a list of abscisses. If this doesn't exist, the first test create a list that will contain a list of indices to allow the function to return a width that will be in fact a difference between two indices
% f: data of the function from wich we want the fwhm
% In addition, the function will return a warning if the fwhm doesn't exist, and return 0 as fwhm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function FWHM = fwhm(x,f)
if nargin < 2
        f = x;
        x=[1:length(f)];
end

[fmax,ifmax]=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 tha "half-mawimum" 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(2)+1))*(x(ind(2)+1)-x(ind(2)))/(f(ind(2)+1)-f(ind(2)))+x(ind(2)+1);
        FWHM=x2-x1;
else
        warning('FWHM is undefined, check if there is an impulsion. FWHM is set to 0')
        FWHM=0;%fwhm is set to 0 if it doesn't exist, to avoid errors in the rest of the programme. The warning is here to tell to the user that there is a problem.
end

---

------------------------------------------------------------------------------

_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

Ok, I think it's good now

Le jeudi 02 juillet 2009 à 22:44 +0200, Søren Hauberg a écrit :

tor, 02 07 2009 kl. 22:36 +0200, skrev Autiwa:
> I think that the last 3 mail (this one count) were only in private
> conversation

Ohh, I didn't notice. I'm CC'ing the list.

> Here is the program with comments in english, maybe less complete, and
> with english error, because I talk easier in french...

Thanks

> I have also added the licence.

You have added:

%## This program is public domain. (GPL)

I hate to be focused on boring license stuff, but we need this to be
correct. The GPL is different from the public domain. So, either put it
in the public domain or release it under the GPL. You have to decide
which you prefer.

Autiwa has suggested that this function goes in the 'signal' package. I
don't really have an opinion on this. What does the rest of the people
on this list think?

Søren

[fwhm.m]

---

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function that return the full width at half maximum
%## This program is public domain.
% x: define a list of abscisses. If this doesn't exist, the first test create a list that will contain a list of indices to allow the function to return a width that will be in fact a difference between two indices
% f: data of the function from wich we want the fwhm
% In addition, the function will return a warning if the fwhm doesn't exist, and return 0 as fwhm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function FWHM = fwhm(x,f)
if nargin < 2
        f = x;
        x=[1:length(f)];
end

[fmax,ifmax]=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 tha "half-mawimum" 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(2)+1))*(x(ind(2)+1)-x(ind(2)))/(f(ind(2)+1)-f(ind(2)))+x(ind(2)+1);
        FWHM=x2-x1;
else
        warning('FWHM is undefined, check if there is an impulsion. FWHM is set to 0')
        FWHM=0;%fwhm is set to 0 if it doesn't exist, to avoid errors in the rest of the programme. The warning is here to tell to the user that there is a problem.
end

---

------------------------------------------------------------------------------

_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[Søren Hauberg]

tor, 02 07 2009 kl. 22:48 +0200, skrev Autiwa:
> Ok, I think it's good now :)

Thanks. Since nobody has commented, I assume people think your stuff is
good :-)

I've (hopefully) made the documentation somewhat more clear, and changed
the coding-style a bit. Could you check that I haven't broken anything?
I'm attaching the code.

By the way, what is your full name? I think you should be given credit
for the code.

Thanks
Søren

[fwhm.m]

---

## -*- 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, an error is raised.
## @seealso{std}
## @end deftypefn

## This program is public domain.

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);
  endif

  if (numel (x) != numel (f))
    error ("fwhm: number of elements does not match");
  endif

  ## 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
    error ("fwhm: full width at half maximum is not defined.");
  endif
endfunction

---

------------------------------------------------------------------------------

_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[Søren Hauberg]

søn, 05 07 2009 kl. 14:10 +0000, skrev Autiwa:
> I modified the .m because I prefer that if fwhm is not defined, that
> fwhm is set to 0. En fact, it happens sometimes that it's a normal
> behaviour for a signal not to have a fwhm. It's disturbing to break
> the run of the program only because of that. I think it's sufficient
> to prompt a warning and set the value to 0. It avoid to have errors in
> the program.

The reason why I changed it is that people tend to complain when they
get a warning, because they cannot easily get rid of it. Errors are
different because you can catch these (using try-catch-blocks). So, I
would say that the function should either return 0 (or something
negative) in which case the behaviour should be documented or raise an
error.

> In addition, I tried to have a syntax that would have allowed someone
> to run it in matlab. But you changed the syntax and now, you can't run
> it under matlab environment.

Oh, I didn't realise you needed Matlab compatibility. Try the attached
version then.

> By the way, you have modified a lot of things so I don't think that
> all credit remain only to me.

I've mostly just inserted spaces to make the code more readable.

> My full name : Christophe Cossou

I've listed you as the author.

Søren

> Le dimanche 05 juillet 2009 à 15:48 +0200, Søren Hauberg a écrit :
> > tor, 02 07 2009 kl. 22:48 +0200, skrev Autiwa:
> > > Ok, I think it's good now :)
> >
> > Thanks. Since nobody has commented, I assume people think your stuff is
> > good :-)
> >
> > I've (hopefully) made the documentation somewhat more clear, and changed
> > the coding-style a bit. Could you check that I haven't broken anything?
> > I'm attaching the code.
> >
> > By the way, what is your full name? I think you should be given credit
> > for the code.
> >
> > Thanks
> > Søren

[fwhm.m]

---

%% -*- 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

---

------------------------------------------------------------------------------

_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[Søren Hauberg]

Hi

Since nobody had any comments, and you didn't raise any further concerns
about Matlab compatibility, I've added the 'fwhm' function to SVN.

Søren


------------------------------------------------------------------------------
Enter the BlackBerry Developer Challenge  
This is your chance to win up to $100,000 in prizes! For a limited time,
vendors submitting new applications to BlackBerry App World(TM) will have
the opportunity to enter the BlackBerry Developer Challenge. See full prize  
details at: http://p.sf.net/sfu/Challenge
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[Petr Mikulik]

> > Ok, I think it's good now :)
>
> Thanks. Since nobody has commented, I assume people think your stuff is
> good :-)
>
> I've (hopefully) made the documentation somewhat more clear, and changed
> the coding-style a bit. Could you check that I haven't broken anything?
> I'm attaching the code.

I've tried it now, and I think the code is buggy:

A.
        x=-pi:0.01:pi; y=cos(x);
        plot(x,y)
        fwhm(x,y)
=> FWHM is pi, but the script returns 2.094

B.
        x=-pi:0.01:pi; y=cos(x) + 5;
        plot(x,y)
        fwhm(x,y)
=> FWHM is pi, but the script says "full width at half maximum is not
=> defined.".

It looks like the script works not between min and max values of the curve
(x,y), but between max(y) and zero. That's wrong and contrary to its
documentation.

I think that "#!test" case should be added.

Further, could it produce vector of fwhm's for matricial f in fwhm(x,f)?

Greetings,
  Petr

------------------------------------------------------------------------------
Enter the BlackBerry Developer Challenge  
This is your chance to win up to $100,000 in prizes! For a limited time,
vendors submitting new applications to BlackBerry App World(TM) will have
the opportunity to enter the BlackBerry Developer Challenge. See full prize  
details at: http://p.sf.net/sfu/Challenge
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

The function will not work if there is many impulsion (maximum) in the vector imput.

To be honest, I didn't know that cosinus had a fwhm. I just code that function for my own use and share it. I will not improve this to work on every mathematical function it already exist. I haven't enought knowledge in math or coding. I only work with impulsion defined from 0 to maximum, it's only a single maximum, not even another impulsion with a lower maximum. I know that with many others shapes it will not work. That's why I wanted the code to return 0 and print a warning when the function would not be able to find a fwhm.

But I think you're right about one thing, the fwhm is defined by half the value between the min and the max. I can't modify it because I haven't the last version of my own function. I have a version for my own use wich is quite different from the version you use (because you have added many test and feature that I don't care for my use and in addition, I want for my use, to have a warning and not an error, so I deleted the public version in my folder )

If you can send me the last version or modify by your own, thank you. I think you juste have to replace 0 by min(f) but I don't have the code in mind.

Regards,
christophe

PS : I don't think it would work for matrix, and however, I don't know how make it works.

Le lundi 13 juillet 2009 à 12:01 +0200, Petr Mikulik a écrit :

> > Ok, I think it's good now :)
> 
> Thanks. Since nobody has commented, I assume people think your stuff is
> good :-)
> 
> I've (hopefully) made the documentation somewhat more clear, and changed
> the coding-style a bit. Could you check that I haven't broken anything?
> I'm attaching the code.

I've tried it now, and I think the code is buggy:

A.
	x=-pi:0.01:pi; y=cos(x);
	plot(x,y)
	fwhm(x,y)
=> FWHM is pi, but the script returns 2.094

B.
	x=-pi:0.01:pi; y=cos(x) + 5;
	plot(x,y)
	fwhm(x,y)
=> FWHM is pi, but the script says "full width at half maximum is not 
=> defined.".

It looks like the script works not between min and max values of the curve 
(x,y), but between max(y) and zero. That's wrong and contrary to its 
documentation.

I think that "#!test" case should be added.

Further, could it produce vector of fwhm's for matricial f in fwhm(x,f)?

Greetings,
  Petr

------------------------------------------------------------------------------
Enter the BlackBerry Developer Challenge  
This is your chance to win up to $100,000 in prizes! For a limited time,
vendors submitting new applications to BlackBerry App World(TM) will have
the opportunity to enter the BlackBerry Developer Challenge. See full prize  
details at: http://p.sf.net/sfu/Challenge
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[Petr Mikulik]

> The function will not work if there is many impulsion (maximum) in the
> vector imput.
>
> To be honest, I didn't know that cosinus had a fwhm. I just code that
> function for my own use and share it. I will not improve this to work on
> every mathematical function it already exist.

Well, it is necessary to test a (contributed) code into all possible inputs
users can supply. I've noticed that your code does neither work for
continuous and monotonous functions.

The enclosed script  fwhm.m  is rewritten from scratch (I was also using
such a function earlier).
- It works on both vectorial and matricial input.
- It returns -1 if FWHM does not exist.
- It was tested against various inputs.
- It works on both Octave and Matlab.
- It includes test case for Octave.

Please test it and commit to Octave if it is found OK (with whatever
comments in the file header).

---
Petr Mikulik

------------------------------------------------------------------------------
Enter the BlackBerry Developer Challenge  
This is your chance to win up to $100,000 in prizes! For a limited time,
vendors submitting new applications to BlackBerry App World(TM) will have
the opportunity to enter the BlackBerry Developer Challenge. See full prize  
details at: http://p.sf.net/sfu/Challenge
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[Søren Hauberg]

man, 13 07 2009 kl. 16:29 +0200, skrev Petr Mikulik:
> The enclosed script  fwhm.m  is rewritten from scratch (I was also using
> such a function earlier).

It seems you forgot to attach the file.

Søren


------------------------------------------------------------------------------
Enter the BlackBerry Developer Challenge  
This is your chance to win up to $100,000 in prizes! For a limited time,
vendors submitting new applications to BlackBerry App World(TM) will have
the opportunity to enter the BlackBerry Developer Challenge. See full prize  
details at: http://p.sf.net/sfu/Challenge
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[Petr Mikulik]

> The function will not work if there is many impulsion (maximum) in the
> vector imput.
>
> To be honest, I didn't know that cosinus had a fwhm. I just code that
> function for my own use and share it. I will not improve this to work on
> every mathematical function it already exist.

Well, it is necessary to test a (contributed) code into all possible inputs
users can supply. I've noticed that your code does neither work for
continuous and monotonous functions.

The enclosed script  fwhm.m  is rewritten from scratch (I was also using
such a function earlier).
- It works on both vectorial and matricial input.
- It returns -1 if FWHM does not exist.
- It was tested against various inputs.
- It works on both Octave and Matlab.
- It includes test case for Octave.

Please test it and commit to Octave if it is found OK (with whatever
comments in the file header).

---
Petr Mikulik
%% Calculation of full-width at half maximum (FWHM) for y or y(x), i.e. full
%% width of the curve between 0 and max(y). Return -1 if FWHM does not exist
%% (e.g. monotonous function).
%%
%% Usage:
%% f = fwhm(y)
%% f = fwhm(x,y)
%% where x is vector and y is vector or matrix.
%% For a matrix argument, return fwhm for each column as a row vector.
%%
%% Compatibility: Octave 3.x, Matlab
%% Version: 13. 7. 2009
%% This program is public domain.


function myfwhm = fwhm (x, y)

    if nargin < 1 || nargin > 2
        error('1 or 2 input arguments required');
    elseif nargin==1
        y = x;
        x = 1:length(y);
    end

    [nr, nc] = size(y);
    if (nr == 1 && nc > 1)
        y = y'; nr = nc; nc = 1;
    end

    if length(x) ~= nr
        error('dimension of input arguments do not match');
    end

    % Shift matrix columns so that y(+-xfwhm) = 0:
    y = y - 0.5 * repmat(max(y), nr, 1); % full-width at half maximum
    % y = y - 0.5 * repmat((max(y) + min(y)), nr, 1); % full-width above background case

    if 0 % Try to do it via a "vectorizing" approach
        ind = find (y(1:end-1, :) .* y(2:end, :) <= 0);
        [r1,c1] = ind2sub(size(y), ind);
        % ... how to take min and max value in each column?
        % => need to do it column-by-column
        myfwhm = -1 + zeros(1,nc);
    else % Analyze the matrix column by column
        myfwhm = -1 + zeros(1,nc); % default: -1 for fwhm undefined
        for n=1:nc
            yy = y(:, n);
            ind = find((yy(1:end-1) .* yy(2:end)) <= 0);
            if length(ind) >= 2 && yy(ind(1)) > 0 % must begin below zero
                ind = ind(2:end);
            end
            [mx, imax] = max(yy); % protection against constant or (almost) monotonous functions
            if length(ind) >= 2 && imax > ind(1) && imax < ind(end)
                ind1 = ind(1);
                ind2 = ind1 + 1;
                xx1 = x(ind1) - yy(ind1) * (x(ind2) - x(ind1)) / (yy(ind2) - yy(ind1));
                ind1 = ind(end);
                ind2 = ind1 + 1;
                xx2 = x(ind1) - yy(ind1) * (x(ind2) - x(ind1)) / (yy(ind2) - yy(ind1));
                myfwhm(n) = xx2 - xx1;
            end
        end
    end

%!test
%! x=-pi:0.001:pi; y=cos(x);
%! assert( abs(fwhm(x, y) - 2*pi/3) < 0.01 );
%!
%!test
%! assert( fwhm(-10:10) == -1 );
%!
%!test
%! assert( fwhm(ones(1,50)) == -1 );
%!
%!test
%! x=-20:1:20;
%! y1=-4+zeros(size(x)); y1(4:10)=8;
%! y2=-2+zeros(size(x)); y2(4:11)=2;
%! y3= 2+zeros(size(x)); y3(5:13)=10;
%! assert( max(abs(fwhm(x, [y1;y2;y3]') - [20.0/3,7.5,9.25])) < 0.01 );

------------------------------------------------------------------------------
Enter the BlackBerry Developer Challenge  
This is your chance to win up to $100,000 in prizes! For a limited time,
vendors submitting new applications to BlackBerry App World(TM) will have
the opportunity to enter the BlackBerry Developer Challenge. See full prize  
details at: http://p.sf.net/sfu/Challenge
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

The code doesn't work for me. In fact, there are two cases where fwhm is set to 0. With the previou code, I had the right fwhm. That's why I don't like a very complex code, because it has strange behavious for ome cases.

My function wa very simple, did only one thing, but did it well.

My program is too complicated and in fact I must keep secret what it contains but the fact is the function you have coded don't reach to find as many fwhm as the previous code.

Your wish is good, but I don't thing it's important to include as many features as this curent version. I think matricial input are not needed. I really think that the most important is to provide a simple function that return the fwhm of a function and that's all. No need to return the fwhm in 5 dimensional space with relativity corrections

By the way, do as you wish, the code is public domain, but I will keep my simple version for my own use

Now, your version is too complicated for me to help you find the bug. I tried to find why it return 0 in place of the right fwhm. Hope you'll find the bug. good luck

Le lundi 13 juillet 2009 à 16:37 +0200, Petr Mikulik a écrit :

> The function will not work if there is many impulsion (maximum) in the
> vector imput. 
> 
> To be honest, I didn't know that cosinus had a fwhm. I just code that
> function for my own use and share it. I will not improve this to work on
> every mathematical function it already exist.

Well, it is necessary to test a (contributed) code into all possible inputs 
users can supply. I've noticed that your code does neither work for 
continuous and monotonous functions.

The enclosed script  fwhm.m  is rewritten from scratch (I was also using 
such a function earlier).
- It works on both vectorial and matricial input. 
- It returns -1 if FWHM does not exist.
- It was tested against various inputs.
- It works on both Octave and Matlab.
- It includes test case for Octave.

Please test it and commit to Octave if it is found OK (with whatever 
comments in the file header).

---
Petr Mikulik

------------------------------------------------------------------------------
Enter the BlackBerry Developer Challenge  
This is your chance to win up to $100,000 in prizes! For a limited time,
vendors submitting new applications to BlackBerry App World(TM) will have
the opportunity to enter the BlackBerry Developer Challenge. See full prize  
details at: http://p.sf.net/sfu/Challenge
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

here is the function I use. On it, I corrected the way to calculate fwhm, and with cos, it works now, try this line :

x=pi:0.01:3*pi; plot(x,cos(x)); fwhm(x,cos(x))

Hope it will help you

Regards,
Christophe

Le lundi 13 juillet 2009 à 16:37 +0200, Petr Mikulik a écrit :

> The function will not work if there is many impulsion (maximum) in the
> vector imput. 
> 
> To be honest, I didn't know that cosinus had a fwhm. I just code that
> function for my own use and share it. I will not improve this to work on
> every mathematical function it already exist.

Well, it is necessary to test a (contributed) code into all possible inputs 
users can supply. I've noticed that your code does neither work for 
continuous and monotonous functions.

The enclosed script  fwhm.m  is rewritten from scratch (I was also using 
such a function earlier).
- It works on both vectorial and matricial input. 
- It returns -1 if FWHM does not exist.
- It was tested against various inputs.
- It works on both Octave and Matlab.
- It includes test case for Octave.

Please test it and commit to Octave if it is found OK (with whatever 
comments in the file header).

---
Petr Mikulik

[fwhm.m]

---

%% -*- 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);
  fmin = min(f);
  renorm = 0.5 * (fmax-fmin) + fmin;
  f_renorm = f - renorm;

  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 = (renorm - f (ind (1) + 1)) * (x (ind (1) + 1) - x (ind (1))) ...
       / (f (ind (1) + 1) - f (ind (1))) + x (ind (1) + 1);
    x2 = (renorm - 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

---

------------------------------------------------------------------------------
Enter the BlackBerry Developer Challenge  
This is your chance to win up to $100,000 in prizes! For a limited time,
vendors submitting new applications to BlackBerry App World(TM) will have
the opportunity to enter the BlackBerry Developer Challenge. See full prize  
details at: http://p.sf.net/sfu/Challenge
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[Petr Mikulik]

> The code doesn't work for me. In fact, there are two cases where fwhm is set
> to 0. With the previou code, I had the right fwhm. That's why I don't like a
> very complex code, because it has strange behavious for ome cases.

I have tried many testing data and found no fail. If you think you have
found data the routine does not work on, please show them.

BTW, your function fails for parabola-like functions, e.g.
        x=-20:1:20;
        y = x.*x;
        fwhm(x,y)
My code tests for this (maybe that's why you find it complicated)?

> here is the function I use. On it, I corrected the way to calculate
> fwhm, and with cos, it works now, try this line

I think we should stay with the original definition of FWHM, i.e. values
taken at fmax*0.5, not at 0.5*(fmax+fmin). The "background" can be
subtracted outside. Or even better, an option can be added. So you should
use
  renorm = 0.5 * fmax;
instead of
  renorm = 0.5 * (fmax-fmin) + fmin;
and fix documentation in the header.

BTW: it is easily
 renorm = 0.5 * (fmax+fmin);

> I think matricial input are not needed. I really think that the most
> important is to provide a simple function that return the fwhm of a
> function and that's all. No need to return the fwhm in 5 dimensional space
> with relativity corrections.

Well, the routine works on data, not on functions. Consider an experiment
with measured 100 data sets put columnwise into a matrix y. Then max(y),
min(y), mean(y), std(y), etc. give statistics for each data set at one call.
Therefore fwhm(y) should behave the same.

I'm enclosing an update for my previous routine fwhm.m that does work on
matrix data, works on both potential definitions of fwhm, and passes all of
the accompanying tests. I consider this being general and powerful enough
for being enclosed into an octave's package.

---
PM
%% Compute full-width at half maximum (FWHM) for vector or matrix data y,
%% optionally sampled as y(x). If y is a matrix, return fwhm for each column
%% as a row vector.
%% f = fwhm(y)
%% f = fwhm(x, y)
%% f = fwhm(x, y, 'zero')
%% f = fwhm(x, y, 'min')
%%
%% The default option 'zero' computes fwhm at half maximum, i.e. 0.5*max(y).
%% The option 'min' computes fwhm at the middle curve, i.e. 0.5*(min(y)+max(y)).
%%
%% Return 0 if FWHM does not exist (e.g. monotonous function or the function
%% does not cut horizontal line at 0.5*max(y) or 0.5(max(y)+min(y)),
%% respectively).
%%
%% Compatibility: Octave 3.x, Matlab
%% Author: Petr Mikulik
%% Version: 15. 7. 2009
%% This program is public domain.


function myfwhm = fwhm (x, y, opt)

    if nargin < 1 || nargin > 3
        error('1, 2 or 3 input arguments required');
    end
    if nargin==1
        y = x;
        x = 1:length(y);
        opt = 'zero';
    elseif nargin==2
        opt = 'zero';
    end

    if ~isstr(opt)
        error('opt must be "zero" or "min"');
    end
    if ~(strcmp(opt, 'zero') || strcmp(opt, 'min'))
        error('opt must be "zero" or "min"');
    end

    [nr, nc] = size(y);
    if (nr == 1 && nc > 1)
        y = y'; nr = nc; nc = 1;
    end

    if length(x) ~= nr
        error('dimension of input arguments do not match');
    end

    % Shift matrix columns so that y(+-xfwhm) = 0:
    if strcmp(opt, 'zero')
        % case: full-width at half maximum
        y = y - 0.5 * repmat(max(y), nr, 1);
    else
        % case: full-width above background
        y = y - 0.5 * repmat((max(y) + min(y)), nr, 1);
    end

    if 0 % Trying to do it via a "vectorizing" approach failed:
        myfwhm = zeros(1,nc); % default: 0 for fwhm undefined
        ind = find (y(1:end-1, :) .* y(2:end, :) <= 0);
        [r1,c1] = ind2sub(size(y), ind);
        % ... how to take min and max value in each column?
        % => need to do it column-by-column
    else % Analyze the matrix column by column
        myfwhm = zeros(1,nc); % default: 0 for fwhm undefined
        for n=1:nc
            yy = y(:, n);
            ind = find((yy(1:end-1) .* yy(2:end)) <= 0);
            if length(ind) >= 2 && yy(ind(1)) > 0 % must start ascending
                ind = ind(2:end);
            end
            [mx, imax] = max(yy); % protection against constant or (almost) monotonous functions
            if length(ind) >= 2 && imax > ind(1) && imax < ind(end)
                ind1 = ind(1);
                ind2 = ind1 + 1;
                xx1 = x(ind1) - yy(ind1) * (x(ind2) - x(ind1)) / (yy(ind2) - yy(ind1));
                ind1 = ind(end);
                ind2 = ind1 + 1;
                xx2 = x(ind1) - yy(ind1) * (x(ind2) - x(ind1)) / (yy(ind2) - yy(ind1));
                myfwhm(n) = xx2 - xx1;
            end
        end
    end

%!test
%! x=-pi:0.001:pi; y=cos(x);
%! assert( abs(fwhm(x, y) - 2*pi/3) < 0.01 );
%!
%!test
%! assert( fwhm(-10:10) == 0 );
%!
%!test
%! assert( fwhm(ones(1,50)) == 0 );
%!
%!test
%! x=-20:1:20;
%! y1=-4+zeros(size(x)); y1(4:10)=8;
%! y2=-2+zeros(size(x)); y2(4:11)=2;
%! y3= 2+zeros(size(x)); y3(5:13)=10;
%! assert( max(abs(fwhm(x, [y1;y2;y3]') - [20.0/3,7.5,9.25])) < 0.01 );
%!
%!test
%! x=-10:10; assert( fwhm(x.*x) == 0 );
%!
%!test
%! x=-5:5; assert( abs( fwhm(x, 18-x.*x, 'zero') - 6.0 ) < 0.01);
%!
%!test
%! x=-5:5; assert( abs( fwhm(x, 18-x.*x, 'min') - 7.0 ) < 0.01);

------------------------------------------------------------------------------
Enter the BlackBerry Developer Challenge  
This is your chance to win up to $100,000 in prizes! For a limited time,
vendors submitting new applications to BlackBerry App World(TM) will have
the opportunity to enter the BlackBerry Developer Challenge. See full prize  
details at: http://p.sf.net/sfu/Challenge
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev

[Søren Hauberg]

ons, 15 07 2009 kl. 14:18 +0200, skrev Petr Mikulik:
> Well, the routine works on data, not on functions. Consider an experiment
> with measured 100 data sets put columnwise into a matrix y. Then max(y),
> min(y), mean(y), std(y), etc. give statistics for each data set at one call.
> Therefore fwhm(y) should behave the same.

I agree with this. These kind of API's are very helpful for writing
vectorised code.

> I'm enclosing an update for my previous routine fwhm.m that does work on
> matrix data, works on both potential definitions of fwhm, and passes all of
> the accompanying tests. I consider this being general and powerful enough
> for being enclosed into an octave's package.

I am not knowledgeable enough to be able to decide which of the two
competing implementations are better. Looking at the code, I get the
impression that Petr's version is more stable, but I would appreciate it
if somebody else could make the decision.

Petr, I have a few comments on your code:

    if ~isstr(opt)
        error('opt must be "zero" or "min"');
    end
    if ~(strcmp(opt, 'zero') || strcmp(opt, 'min'))
        error('opt must be "zero" or "min"');
    end

Shouldn't that just be

    if ~isstr(opt) || ~(strcmp(opt, 'zero') || strcmp(opt, 'min'))
        error('opt must be "zero" or "min"');
    end

? Also, in Octave you should be able to use 'strcmp (opt, {'zero',
'min'})', but I don't know if that works in Matlab.

You have some code that is removed using

    if 0
        ...

This should probably be changed into a comment.

Søren


------------------------------------------------------------------------------
Enter the BlackBerry Developer Challenge  
This is your chance to win up to $100,000 in prizes! For a limited time,
vendors submitting new applications to BlackBerry App World(TM) will have
the opportunity to enter the BlackBerry Developer Challenge. See full prize  
details at: http://p.sf.net/sfu/Challenge
_______________________________________________
Octave-dev mailing list
Octave-dev@...
https://lists.sourceforge.net/lists/listinfo/octave-dev