Technical Report no. 8 Ravi Kochhar Mar 22 1988 Department of Physiology Rev. 2.03, Mar. 23 1999 University of Wisconsin (rep008) Madison, Wi. 53706
INTRODUCTION
The library NUMLIB contains a series of general purpose routines which can be called by any program on VAX/VMS or RSX systems. The routines are mostly for numeric computations, and include FFT, Regression, Random number, Matrix inversion and other routines.
On the VAX, the library resides in directory L: and may be specified at LINK time as in the following example :
$ LINK program,L:NUMLIB/L
Individual routines are described below in more detail.
PROGRAMMING NOTES
(1) As with all numerical computations, accuracy depends on many factors, including the numbers supplied to the routines in a particular application. Thus no claim is made as to the accuracy or reliability of these routines and they are provided for general use on a strictly as-is basis. Users are advised to test the routines in their particular application before relying on them.
(2) Unless otherwise specified, the typing of variables follows the FORTRAN convention, i.e. names that begin with letters I to N are type INTEGER and the rest are single precision floating point (REAL). Exceptions to this rule are noted with individual routines.
FOURIER TRANSFORM ROUTINES
SUBROUTINE FFA(M,X)
Compute the Fast Fourier Transform of the specified REAL series. M : Power of 2 (N=2**M, where N=no. of points) X : Real array containing N values
Both M and X must be supplied by the calling program. This routine will replace the values in X with the corresponding Fourier Transform. Thus the contents of X will be lost. The number of points in X must be an exact power of 2.
The (N/2+1) Fourier coefficients are returned in array X, arranged as follows :
F(0) = X(1)
F(N/2) = X(2)
F(k) = X(2k+1) + i.X(2k+2)
for k=1,2,3.....(N/2-1)
Thus the real and imaginary parts alternate after the first two values, which are the DC term and the amplitude of the N/2 th harmonic, respectively.
To get amplitude and phase from the the Fourier coefficients returned above, use the following :
If F(k) = X(2k+1) + i.X(2k+2)
then :
Amplitude(rms) = sqrt(2).sqrt(X(2k+1)**2+X(2k+2)**2)
Amplitude(peak)= 2.sqrt(X(2k+1)**2+X(2k+2)**2)
phase(radians) = arctan(X(2k+2)/X(2k+1))
The phase returned by the routine FFA is expressed relative to a Cosine, and a phase lag is a negative angle. Thus, the phase for a Cosine wave will be returned as zero, and for a sine wave will be returned as minus PI/2.
3.2 FFS
SUBROUTINE FFS(M,X)
Compute the Inverse Fast Fourier Transform of specified
series.
M : Power of 2 (N=2**M, where N=no. of points)
X : Real array containing the Fourier Coefficients
Both M and X must be supplied by the calling program. This
routine will replace the coefficients in X with the corresponding
real series. Thus the contents of X will be lost. The number of
points in X must be an exact power of 2.
Note that the arrangement of the coefficients in X must be
exactly as output by the routine FFA (see above).
3.3 FOURI
SUBROUTINE FOURI(N,X,Y)
Compute the straight Fourier Transform of the specified REAL
series.
N : Number of points
X : Array containing real series
Y : Array containing Fourier coefficients
N and X must be supplied, while Y will be returned by this
routine.
The arrangement of the coefficients in Y is exactly as output
by the routine FFA (see above).
Unlike FFA, this routine does not destroy the contents of X.
Also, N need not be a power of 2.
This routine runs considerably slower than FFA for all values
of N greater than 256.
3.4 FOUR3
SUBROUTINE FOUR3(B,R)
Compute the straight Fourier Transform of the specified 3
point REAL series.
B : Input time domain array
R : Output freq. domain array
B must be supplied, while R will be returned by this routine.
B must contain exactly 3 points in the time domain.
The arrangement of the coefficients in R is exactly as output
by the routine FFA (see above).
Unlike FFA, this routine does not destroy the contents of B.
3.5 DFFA
SUBROUTINE DFFA(M,X)
Double precision version of FFA.
This routine is identical to FFA in every respect except that
here X is a double precision array, and all computations within
DFFA are done in double precision arithmetic.
3.6 DFFS
SUBROUTINE DFFS(M,X)
Double precision version of FFS.
This routine is identical to FFS in every respect except that
here X is a double precision array, and all computations within
DFFS are done in double precision arithmetic.
3.7 DFOURI
SUBROUTINE DFOURI(N,X,Y)
Double precision version of FOURI.
This routine is identical to FOURI in every respect except
that here B and R are double precision arrays, and all
computations within DFOURI are done using double precision
arithmetic.
4 HISTOGRAM ROUTINES
4.1 HSTMSK
SUBROUTINE HSTMSK(HIST,NBIN,FVAL,BW,RNS,RM,STDV,SKEW,RKURT,NER)
Compute statistics for specified histogram.
HIST : Real array containing histogram
NBIN : Number of bins
FVAL : Starting value of histogram
BW : Histogram bin width
RNS : Total number of events in histogram
RM : Mean value
STDV : Standard deviation
SKEW : Skewness
RKURT : Kurtosis
NER : Error code, normally zero
HIST,NBIN,FVAL and BW must be supplied by the calling
program, while RNS,RM,STDV,SKEW,RKURT and NER will be returned by
this routine.
Note that RM and STDV will be in units of FVAL and BW.
The array HIST must be dimensioned to at least NBIN.
4.2 HSTPH
SUBROUTINE HSTPH(NHIST,NBIN,NS,SYNC,THETA,NER)
Compute synchronization coefficient and phase for specified
histogram.
NHIST : Integer array containing histogram
NBIN : Number of bins
NS : Total number of events in histogram
SYNC : Sync. coefficient (vector strength) (0 to 1)
THETA : Phase (0 to 1)
NER : Error code, normally zero
NHIST and NBIN must be supplied by the calling program, while
NS,SYNC,THETA and NER will be returned by this routine.
The histogram is supplied as an INTEGER array. If your array
is REAL, use the routine HSTPHR instead (described below).
Page 9
The possible error codes (NER) are as follows :
320 : Invalid number of bins
4.3 HSTPHR
SUBROUTINE HSTPHR(HIST,NBIN,NS,SYNC,THETA,NER)
This routine is identical to HSTPH (above) in every way
except that here the histogram is supplied in the REAL array HIST.
5 RANDOM NUMBER ROUTINES
5.1 GAUSS
SUBROUTINE GAUSS(NSEED,N,X)
Compute a series of normally distributed random numbers
NSEED : Seed for random number generator (INTEGER*4)
N : Number of random numbers to be generated
X : Real array in which random numbers are returned
NSEED and N must be supplied, while X will be returned.
The Mean and Variance of the series have expected values of
zero and one respectively
The value of NSEED will be modified by this routine. A
typical way to generate NSEED might be to call the SECNDS routine
available in VMS/RSX FORTRAN and get the time of day in seconds
since midnight.
This routine uses the routine RAN available in VMS/RSX to
generate uniformly distributed random numbers and then converts
them to a gaussian (normal) distribution using an algorithm
supplied by Bob Stanny.
5.2 POISN
SUBROUTINE POISN(NSEED,N,RLMB,IX)
Compute Poisson distributed random numbers
NSEED : Seed for random number generator (INTEGER*4)
N : Number of values to be generated
RLMB : Mean and variance of random series
IX : Integer array to hold Poisson distributed series
NSEED,N and RLMB must be supplied, while IX will be returned
by this routine.
Note that the numbers are returned in the INTEGER array IX.
The value of NSEED will be modified by this routine. A
typical way to generate NSEED might be to call the SECNDS routine
available in VMS/RSX FORTRAN and get the time of day in seconds
since midnight.
This routine uses the routine RAN available in VMS/RSX to
generate uniformly distributed random numbers and then converts
them to a Poisson distribution.
6 MATRIX ROUTINES
6.1 MINV
SUBROUTINE MINV(A,N,D,L,M)
Invert a matrix
A : Matrix to be inverted (N x N) (real)
N : Order of the matrix
D : Resultant determinant
L : Work vector of length N (integer)
M : Work vector of length N (integer)
A and N must be supplied by the calling program, while A and
D will be returned by this routine.
This routine destroys the contents of A and replaces them by
the inverse.
The standard Gauss-Jordan method is used. The determinant is
also calculated. A value of zero for the determinant means that
the matrix is singular and the inverse cannot be computed.
The routine expects the matrix A to be in "vector storage
format", i.e. each column of the matrix is immediately followed
in storage by the next column. If the number of rows and columns
in the matrix are the same as in your dimension statement, then
this is always true. If N does not equal the value in your
dimension statement then you must first arrange your 2-dimensional
matrix in a single dimension vector so that each column is
followed by the next without any gaps.
7 TIME SERIES ROUTINES
7.1 AUTOC
SUBROUTINE AUTOC(A,N,L,R)
Compute auto-correlation for series "A" for lags of 0 to L-1
A : Single dimensioned real array containing points
N : No. of values in A
L : Auto-correlation computed for lags of 0,1,2...,L-1
R : Array of length L to hold results
A,N and L must be supplied, while R will be returned. A and
R are single precision floating point arrays. N must be greater
than L. If not, R(1) is set to zero and the routine returns
immediately.
7.2 CROSSC
SUBROUTINE CROSSC(A,B,N,L,R,S)
Compute the cross-correlation of series A with series B
A : Real array containing first time series
B : Real array containing second time series
N : No. of points in A and B
L : Compute cross-correlation for lags and leads of 0 to L-1
R : Real array to hold crosscorrelation where B lags A
S : Real array to hold crosscorrelation where B leads A
A,B,N and L must be supplied, while R and S will be returned.
N must be greater than L, else R(1) and S(1) will be set to zero
and the routine will return immediately. A,B,R and S are all
single dimensioned single precision arrays.
7.3 SMTHM
SUBROUTINE SMTHM(X,NP,N,XMP,NER)
Smooth the values in array X using N-point smoothing
X : Real single dimensioned array containing points
NP : No. of points in X
N : Smoothing factor (odd number from 1 to 31)
XMP : Missing Value code. Any values in X equal to XMP
Page 14
are ignored for purposes of smoothing
NER : Error code, normally zero
X,NP,N and XMP must be supplied. NER will be returned. The
smoothed points will replace the original contents of array X. X
must be a single precision floating point, single dimensioned
array.
The routine uses a rectangular window for smoothing.
7.4 RSAMPL
SUBROUTINE RSAMPL(X,MAX,NP,NEWP,NER)
Re-sample points in array X
X : Real single dimensioned array containing points
MAX : Max no. of points X can hold (i.e. dimension of X)
NP : No. of points currently in X
NEWP : No. of points in X after re-sampling
NER : Error code, normally zero
X,MAX,NP and NEWP must be supplied. The re-sampled points
will replace the contents of X.
The routine uses linear interpolation to re-sample the points
in X so there are NEWP points instead of NP. NEWP can be less
than or greater than NP, but it cannot be greater than MAX.
8 OTHER ROUTINES
8.1 SORTA
SUBROUTINE SORTA(X,N)
Sort an array in ascending order
X : Single dimensioned array of real numbers to be sorted
N : No. of values in X
Both X and N must be supplied. X is a single precision real
array. The sorted values replace the original contents of X.
This routine uses the bubble sort technique.
8.2 LININC
SUBROUTINE LININC(IFL,A,B,X1,Y1,RD,X2,Y2,XC,YC,DL,NER)
Determine line coordinates at a particular distance
IFL : If IFL=0 then eqn. is Y=AX+B, if=1 then X=A
A : Slope of line
B : Y-axis intercept of line
X1 : X-coordinate of origin point
Y1 : Y-coordinate of origin point
RD : Distance from (X1,Y1)
X2 : X-coordinate of destination point
Y2 : Y-coordinate of destination point
XC : X-coordinate of point on line
YC : Y-coordinate of point on line
DL : Distance of (XC,YC) from (X2,Y2)
NER : Error code, normally zero
IFL,A,B,X1,Y1,RD,X2,Y2 must be supplied by the calling
program, while XC,YC,DL and NER will be returned.
Given a straight line with the equation "Y=A.X+B" (or, if
IFL=1 then given a line with the equation "X=A"), compute the
coordinates of the point on the line which is a linear distance
"RD" from the point (X1,Y1). Return the point which is closer to
(X2,Y2). The results are returned in (XC,YC). Also return "DL",
the distance from (XC,YC) to (X2,Y2)
8.3 INT2D
SUBROUTINE INT2D(A,B,X1,Y1,X2,Y2,XC,YC,KODE)
Determine if the given straight line intersects another.
A : Slope of first line (whose eqn. is Y=A.X+B)
B : Y-intercept of first line (eqn: Y=A.X+B)
X1 : X-coordinate of first point
Y1 : Y-coordinate of first point
X2 : X-coordinate of second point
Y2 : Y-coordinate of second point
XC : X-coordinate of intersection point
YC : Y-coordinate of intersection point
KODE : returned as KODE=-1 if computation not possible
: KODE=0 if no intersection (lines are parallel)
: KODE=1 if intersection point is between (X1,Y1),(X2,Y2)
: KODE=2 if intersection point is beyond (X1,Y1),(X2,Y2)
A,B,X1,Y1,X2 and Y2 must be supplied, while XC,YC and KODE
will be returned.
This routine determines whether the straight line whose
equation is Y=A.X+B intersects the straight line whose endpoints
are (X1,Y1) and (X2,Y2). If a point of intersection is found then
it is returned in (XC,YC). Further information is returned in
KODE.
8.4 FORDR
SUBROUTINE FORDR(X,Y,N,A,B,NER)
Fit a straight line using least squares
X : Array containing X-coordinates of points to be fitted
Y : Array containing Y-coordinates of points to be fitted
N : No. of points in X,Y
A : Slope of fitted line (Y=A.X+B)
B : Y-intercept of fitted line (Y=A.X+B)
NER : Error code, normally zero
X,Y,and N must be supplied. A,B and NER will be returned. X
and Y are both single dimensioned, single precision, floating
point arrays.
8.5 DFORDR
SUBROUTINE DFORDR(X,Y,N,A,B,NER)
Double precision version of FORDR.
This routine is identical to FORDR in every way except that
here X and Y are double precision arrays, and A and B are also
double precision. All computations within DFORDR are done using
double precision arithmetic.
8.6 RYCOFF
SUBROUTINE RYCOFF(NS,SYNC,RY,NER)
Compute the Rayleigh Coefficient, given the no. of spikes
and sync. coeff.
NS : No. of spikes
SYNC : Sync. coefficient (0.0 to 1.0)
RY : Rayleigh coefficient
NER : Error code, normally zero
NS and SYNC must be supplied. RY and NER are returned. The
routine uses a table lookup get the rayleigh coefficient.
8.7 GCF
SUBROUTINE GCF(N1,N2,NGCF,NER)
Compute the greatest common factor of two integers
N1 : First number
N2 : Second number
NGCF : GCF of (N1,N2)
NER : Error code, normally zero
N1 and N2 must be supplied, while NGCF and NER will be
returned.
The Greatest Common Factor (GCF) is the largest integer
number that divides both N1 and N2. For example, the GCF of 1200
and 1300 is 100.
SUBROUTINE MNSDV(X,N,RMN,RSDV,NER)
Compute mean and standard deviation of a series of real numbers.
X : Array of Real numbers
N : Number of values in X
RMN : Mean of the values in X
RSDV : Standard deviation of the values in X
NER : Error code, normally zero
X and N must be supplied, while RMN, RSDV and NER will be returned by the routine.
SUBROUTINE MNSDVI(IX,N,RMN,RSDV,NER)
Compute mean and standard deviation of a series of integer numbers.
IX : Array of Integer numbers
N : Number of values in IX
RMN : Mean of the values in IX
RSDV : Standard deviation of the values in IX
NER : Error code, normally zero
IX and N must be supplied, while RMN, RSDV and NER will be returned by the routine.
ACKNOWLEDGEMENTS
NUMLIB was written under the direction of Dr. W.S.Rhode.
The routines in NUMLIB have been adapted from a variety of sources.
The routines FFA and FFS have been adapted from the following article :
"A Radix-eight Fast Fourier Transform Subroutine for Real Valued Series" by G.D.Bergland, IEEE trans. Audio and Electroacoustics, June 1969, pp 138-144.
The multiple regression routines are from the book :
"IBM System/360 Scientific Subroutine Package" - Version III
The routine GAUSS was provided by Bob Stanny.
For time-series techniques, the reference is :
"Digital Time Series Analysis" by R.K.Otnes and L.Enochson, Wiley 1972.
Supported in part by a grant from the National Institutes of Health (NIH).