> 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@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/octave-dev
