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. Thusthe 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).

If you have questions or suggestions about this document, please send e-mail to kochhar@physiology.wisc.edu

Return to Computing Page

Return to Basement Page