From sources-request at octave dot org Thu Jun 16 13:03:53 2005
Subject: find_cross.m, version 2.0 (typo fix)
From: "Phillip M. Feldman" <pfeldman at verizon dot net>
To: sources at octave dot org
Date: Thu, 16 Jun 2005 12:53:09 -0500

This is a multi-part message in MIME format.
--------------030807070305080105000909
Content-Type: text/plain; charset=us-ascii; format=flowed
Content-Transfer-Encoding: 7bit

Following a suggestion from Paul Kienzle (National Institute of Standards and Technology),
I've switched to a vectorized algorithm that is much faster when working with large
datasets.  The underlying assumption is that, even for a large dataset, the number of
level crossings tends to be small, so that the initial step of finding array indices in
the vicinity of level crossings is the critical one from the standpoint of performance.
The subsequent step of estimating level crossing locations via interpolation can be done
within a loop.

Phillip



--------------030807070305080105000909
Content-Type: text/plain;
 name="find_cross.m"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="find_cross.m"

%                             find_cross.m
%
%          Phillip M. Feldman, Northrop Grumman Space Technology
%
%                      Version 2.0  June 16, 2005
%
% Overview
%
% find_cross() uses a combination of search and interpolation to
% determine the point or points at which a given level is crossed.  If
% there are multiple level crossings, find_cross finds all of the level
% crossings.  If the input sequence "kisses" the specified level but
% does not actually cross it, this is not counted as a crossing.  If
% there are two or more consecutive zeroes, each of them is counted as a
% crossing.
%
% Calling syntax:
%
% X_cross= find_cross(X, Y, Y_target)
%
% where:
%
% X is a column vector containing a sequence of monotone increasing
% values of the independent variable.
%
% Y is a column vector containing values of the dependent variable.
%
% Y_target is the desired level crossing.
%
% X_cross (the returned value) is the estimated value of X at which the
% dependent variable crosses the level Y_target.
%
%
% Acknowledgment
%
% This code uses a vectorized algorithm that was suggested to me by Paul
% Kienzle of the National Institute of Standards and Technology,
% Gaithersburg, MD.


% Revision History

% 2.0 Vectorized the code for faster execution.  The underlying
% assumption is that, even for a large dataset, the number of level
% crossings tends to be small, so that the initial step of finding array
% indices in the vicinity of level crossings is the critical one from
% the standpoint of performance.  The subsequent step of estimating
% level crossing locations via interpolation can be done within a loop.

% 1.1 Fixed code to handle situation where Y_target exactly matches
% Y(i), with Y(i-1) and Y(i+1) bracketing Y_target.  Previously, this
% was not being recognized as a crossing.

function X_cross= find_cross(X, Y, Y_target)


% Section 1: Preliminaries.

% Check dimensions of input arrays:

[m n]= size(X);
[r p]= size(Y);

if (n ~= 1 | p ~= 1)
   error('X and Y must be column vectors.');
end

if (m ~= r)
   error('Lengths of X and Y vectors must be the same.');
end

% Initialize output array to guarantee that it is defined:

X_cross= [];

% To prevent a subscript-out-of-bounds error in Section 2, augment Y by
% duplicating the last element of the array.  (We are operating on a
% copy of the Y array, so the Y array in the calling program is not
% affected).

Y(m+1)= Y(m);

% Subtract Y_target to simplify subsequent code:

Y= Y - Y_target;


% Section 2: Find array indices in crossing vicinities.

% Find potential crossings:

idx= find ( Y(1:end-2) .* Y(2:end-1) <= 0 )

% At this point, the largest possible element of idx is m-1, and the
% Y array has m+1 elements.

% Eliminate "kisses", i.e., cases where a zero is surrounded by two
% same-sign values:

idx(Y(idx+1)==0 & Y(idx).*Y(idx+2) > 0) = [];


% Section 3: Examine each of the crossing vicinities to determine
% whether the crossing in question is an exact crossing.  If it is not
% an exact crossing, use linear interpolation to estimate crossing
% point.  (idx should be a very small array, so we can afford to use a
% loop).

for i= 1 : size(idx,1)

   j= idx(i);


   % There are now two cases to consider:

   % Case 1: (Non-kiss) zero ==> Exact crossing.

   if (Y(j) == 0)
      X_cross(i)= X(j);


   % Case 2: Level crossing occurs between X(j) and X(j+1).

   else

      % Estimate crossing point via linear interpolation:

      X_cross(i)= ...
        X(j) + (X(j+1) - X(j)) * Y(j) / (Y(j) - Y(j+1));

   end

end

--------------030807070305080105000909--