NUMLIB

Numerical Routines Library

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

Back to Top


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
This page last modified on : Mar. 23, 1999