Octave Miscellaneous Functions



Miscellaneous Functions (Not yet properly filed/documented)
===========================================================

 - Function File:  axis2dlim (AXDATA)
     determine axis limits for 2-d data(column vectors); leaves a 10%
     margin around the plots.  puts in margins of +/- 0.1 if data is
     one dimensional (or a single point)

     *Inputs* AXDATA nx2 matrix of data [x,y]

     *Outputs* AXVEC vector of axis limits appropriate for call to
     axis() function

 - Function File:  moddemo (INPUTS)
     Octave Controls toolbox demo: Model Manipulations demo

 - Function File:  prompt (INPUTS)
     function prompt([str])
     Prompt user to continue
     str: input string. Default value: "\n ---- Press a key to continue ---"

 - Function File:  rldemo (INPUTS)
     Octave Controls toolbox demo: Root Locus demo

 - Function File:  rlocus (INPUTS)
     [rldata, k] = rlocus(sys[,increment,min_k,max_k])
     Displays root locus plot of the specified SISO system.
     
            -----   ---     --------
        --->| + |---|k|---->| SISO |----------->
            -----   ---     --------        |
            - ^                             |
              |_____________________________|
     
     inputs: sys = system data structure
     min_k, max_k,increment: minimum, maximum values of k and
     the increment used in computing gain values
     Outputs: plots the root locus to the screen.
     rldata: Data points plotted column 1: real values, column 2: imaginary
     values)
     k: gains for real axis break points.

 - Function File:  sortcom (INPUTS)
     [yy,idx] = sortcom(xx[,opt]): sort a complex vector
     xx: complex vector
     opt: sorting option:
      "re": real part (default)
      "mag": by magnitude
      "im": by imaginary part
     
     if opt != "im" then complex conjugate pairs are grouped together,
     a - jb followed by a + jb.
     yy: sorted values
     idx: permutation vector: yy = xx(idx)

 - Function File:  ss2tf (INPUTS)
     [num,den] = ss2tf(a,b,c,d)
     Conversion from tranfer function to state-space.
     The state space system
           .
           x = Ax + Bu
           y = Cx + Du
     
     is converted to a transfer function
     
                     num(s)
               G(s)=-------
                     den(s)
     
     used internally in system data structure format manipulations

 - Function File:  ss2zp (INPUTS)
     Converts a state space representation to a set of poles and zeros.
     
     [pol,zer,k] = ss2zp(a,b,c,d) returns the poles and zeros of the state space
     system (a,b,c,d).  K is a gain associated with the zeros.
     
     used internally in system data structure format manipulations

 - Function File:  starp (INPUTS)

     sys = starp(P, K, ny, nu)
     
     Redheffer star product or upper/lower LFT, respectively.
     
     
                    +-------+
          --------->|       |--------->
                    |   P   |
               +--->|       |---+  ny
               |    +-------+   |
               +-------------------+
                                |  |
               +----------------+  |
               |                   |
               |    +-------+      |
               +--->|       |------+ nu
                    |   K   |
          --------->|       |--------->
                    +-------+
     
     If ny and nu "consume" all inputs and outputs of K then the result
     is a lower fractional transformation. If ny and nu "consume" all
     inputs and outputs of P then the result is an upper fractional
     transformation.
     
     ny and/or nu may be negative (= negative feedback)

 - Function File:  tf2ss (INPUTS)
     Conversion from tranfer function to state-space.
     The state space system
           .
           x = Ax + Bu
           y = Cx + Du
     
     is obtained from a transfer function
     
                     num(s)
               G(s)=-------
                     den(s)
     
     via the function call [a,b,c,d] = tf2ss(num,den).
     The vector 'den' must contain only one row, whereas the vector 'num'
     may contain as many rows as there are outputs of the system 'y'.
     The state space system matrices obtained from this function will be
     in controllable canonical form as described in "Modern Control Theory",
     [Brogan, 1991].

 - Function File:  tf2zp (INPUTS)
     Converts transfer functions to poles / zeros.

     [zer,pol,k] = tf2zp(num,den) returns the zeros and poles of the
     SISO system defined by num/den.  K is a gain associated with the
     system zeros.

 - Function File: [A, B, C, D] = zp2ss (ZER, POL, K)
     Conversion from zero / pole to state space.  *Inputs*
    ZER
    POL
          vectors of (possibly) complex poles and zeros of a transfer
          function.  Complex values must come in conjugate pairs (i.e.,
          x+jy in zer means that x-jy is also in zer)

    K
          real scalar (leading coefficient) *Outputs* A, B, C, D The
     state space system
          .
          x = Ax + Bu
          y = Cx + Du
     is obtained from a vector of zeros and a vector of poles via the
     function call `[a,b,c,d] = zp2ss(zer,pol,k)'.  The vectors `zer'
     and `pol' may either be row or column vectors.  Each zero and pole
     that has an imaginary part must have a conjugate in the list.  The
     number of zeros must not exceed the number of poles.  `k' is
     `zp'-form leading coefficient.

 - Function File: [NUM, DEN] = zp2tf (ZER, POL, K)
     Converts zeros / poles to a transfer function.  *Inputs*
    ZER
    POL
          vectors of (possibly complex) poles and zeros of a transfer
          function.  Complex values should appear in conjugate pairs

    K
          real scalar (leading coefficient) `[num,den] =
     zp2tf(zer,pol,k)' forms the transfer function `num/den' from the
     vectors of poles and zeros.


Signal Processing
*****************

   I hope that someday Octave will include more signal processing
functions.  If you would like to help improve Octave in this area,
please contact .

 - Function File:  detrend (X, P)
     If X is a vector, `detrend (X, P)' removes the best fit of a
     polynomial of order P from the data X.

     If X is a matrix, `detrend (X, P)' does the same for each column
     in X.

     The second argument is optional.  If it is not specified, a value
     of 1 is assumed.  This corresponds to removing a linear trend.

 - Loadable Function:  fft (A, N)
     Compute the FFT of A using subroutines from FFTPACK.  If A is a
     matrix, `fft' computes the FFT for each column of A.

     If called with two arguments, N is expected to be an integer
     specifying the number of elements of A to use.  If A is a matrix,
     N specifies the number of rows of A to use.  If N is larger than
     the size of A, A is resized and padded with zeros.

 - Loadable Function:  ifft (A, N)
     Compute the inverse FFT of A using subroutines from FFTPACK.  If A
     is a matrix, `fft' computes the inverse FFT for each column of A.

     If called with two arguments, N is expected to be an integer
     specifying the number of elements of A to use.  If A is a matrix,
     N specifies the number of rows of A to use.  If N is larger than
     the size of A, A is resized and padded with zeros.

 - Loadable Function:  fft2 (A, N, M)
     Compute the two dimensional FFT of A.

     The optional arguments N and M may be used specify the number of
     rows and columns of A to use.  If either of these is larger than
     the size of A, A is resized and padded with zeros.

 - Loadable Function:  ifft2 (A, N, M)
     Compute the two dimensional inverse FFT of A.

     The optional arguments N and M may be used specify the number of
     rows and columns of A to use.  If either of these is larger than
     the size of A, A is resized and padded with zeros.

 - Function File:  fftconv (A, B, N)
     Return the convolution of the vectors A and B, as a vector with
     length equal to the `length (a) + length (b) - 1'.  If A and B are
     the coefficient vectors of two polynomials, the returned value is
     the coefficient vector of the product polynomial.

     The computation uses the FFT by calling the function `fftfilt'.  If
     the optional argument N is specified, an N-point FFT is used.

 - Function File:  fftfilt (B, X, N)
     With two arguments, `fftfilt' filters X with the FIR filter B
     using the FFT.

     Given the optional third argument, N, `fftfilt' uses the
     overlap-add method to filter X with B using an N-point FFT.

 - Loadable Function: y = filter (B, A, X)
 - Loadable Function: [Y, SF] = filter (B, A, X, SI)
     Return the solution to the following linear, time-invariant
     difference equation:

             N                   M
            SUM a(k+1) y(n-k) = SUM b(k+1) x(n-k)      for 1<=n<=length(x)
            k=0                 k=0

     where  N=length(a)-1 and M=length(b)-1.  An equivalent form of
     this equation is:

                      N                   M
            y(n) = - SUM c(k+1) y(n-k) + SUM d(k+1) x(n-k)  for 1<=n<=length(x)
                     k=1                 k=0

     where  c = a/a(1) and d = b/a(1).

     If the fourth argument SI is provided, it is taken as the initial
     state of the system and the final state is returned as SF.  The
     state vector is a column vector whose length is equal to the
     length of the longest coefficient vector minus one.  If SI is not
     supplied, the initial state vector is set to all zeros.

     In terms of the z-transform, y is the result of passing the
     discrete- time signal x through a system characterized by the
     following rational system function:

                       M
                      SUM d(k+1) z^(-k)
                      k=0
            H(z) = ----------------------
                         N
                    1 + SUM c(k+1) z(-k)
                        k=1

 - Function File: [H, W] = freqz (B, A, N, "whole")
     Return the complex frequency response H of the rational IIR filter
     whose numerator and denominator coefficients are B and A,
     respectively.  The response is evaluated at N angular frequencies
     between 0 and  2*pi.

     The output value W is a vector of the frequencies.

     If the fourth argument is omitted, the response is evaluated at
     frequencies between 0 and  pi.

     If N is omitted, a value of 512 is assumed.

     If A is omitted, the denominator is assumed to be 1 (this
     corresponds to a simple FIR filter).

     For fastest computation, N should factor into a small number of
     small primes.

 - Function File:  sinc (X)
     Return  sin(pi*x)/(pi*x).

 - Function File: [A, B] = arch_fit (Y, X, P, ITER, GAMMA, A0, B0)
     Fit an ARCH regression model to the time series Y using the
     scoring algorithm in Engle's original ARCH paper.  The model is

          y(t) = b(1) * x(t,1) + ... + b(k) * x(t,k) + e(t),
          h(t) = a(1) + a(2) * e(t-1)^2 + ... + a(p+1) * e(t-p)^2

     in which e(t) is N(0, h(t)), given a time-series vector Y up to
     time t-1 and a matrix of (ordinary) regressors X up to t.  The
     order of the regression of the residual variance is specified by P.

     If invoked as `arch_fit (Y, K, P)' with a positive integer K, fit
     an ARCH(K, P) process, i.e., do the above with the t-th row of X
     given by

          [1, y(t-1), ..., y(t-k)]

     Optionally, one can specify the number of iterations ITER, the
     updating factor GAMMA, and initial values a0 and b0 for the
     scoring algorithm.

 - Function File:  arch_rnd (A, B, T)
     Simulate an ARCH sequence of length T with AR coefficients B and
     CH coefficients A.  I.e., the result y(t) follows the model

          y(t) = b(1) + b(2) * y(t-1) + ... + b(lb) * y(t-lb+1) + e(t),

     where e(t), given Y up to time t-1, is N(0, h(t)), with

          h(t) = a(1) + a(2) * e(t-1)^2 + ... + a(la) * e(t-la+1)^2

 - Function File:  [PVAL, LM] = arch_test (Y, X, P)
     For a linear regression model

          y = x * b + e

     perform a Lagrange Multiplier (LM) test of the null hypothesis of
     no conditional heteroscedascity against the alternative of CH(P).

     I.e., the model is

          y(t) = b(1) * x(t,1) + ... + b(k) * x(t,k) + e(t),

     given Y up to t-1 and X up to t, e(t) is N(0, h(t)) with

          h(t) = v + a(1) * e(t-1)^2 + ... + a(p) * e(t-p)^2,

     and the null is a(1) == ... == a(p) == 0.

     If the second argument is a scalar integer, k, perform the same
     test in a linear autoregression model of order k, i.e., with

          [1, y(t-1), ..., y(t-K)]

     as the t-th row of X.

     Under the null, LM approximately has a chisquare distribution with
     P degrees of freedom and PVAL is the p-value (1 minus the CDF of
     this distribution at LM) of the test.

     If no output argument is given, the p-value is displayed.

 - Function File:  arma_rnd (A, B, V, T, N)
     Return a simulation of the ARMA model

          x(n) = a(1) * x(n-1) + ... + a(k) * x(n-k)
               + e(n) + b(1) * e(n-1) + ... + b(l) * e(n-l)

     in which K is the length of vector A, L is the length of vector B
     and E is gaussian white noise with variance V.  The function
     returns a vector of length T.

     The optional parameter N gives the number of dummy X(I) used for
     initialization, i.e., a sequence of length T+N is generated and
     X(N+1:T+N) is returned.  If N is omitted, N = 100 is used.

 - Function File:  autocor (X, H)
     Return the autocorrelations from lag 0 to H of vector X.  If H is
     omitted, all autocorrelations are computed.  If X is a matrix, the
     autocorrelations of each column are computed.

 - Function File:  autocov (X, H)
     Return the autocovariances from lag 0 to H of vector X.  If H is
     omitted, all autocovariances are computed.  If X is a matrix, the
     autocovariances of each column are computed.

 - Function File:  autoreg_matrix (Y, K)
     Given a time series (vector) Y, return a matrix with ones in the
     first column and the first K lagged values of Y in the other
     columns.  I.e., for T > K, `[1, Y(T-1), ..., Y(T-K)]' is the t-th
     row of the result.  The resulting matrix may be used as a
     regressor matrix in autoregressions.

 - Function File:  bartlett (M)
     Return the filter coefficients of a Bartlett (triangular) window of
     length M.

     For a definition of the Bartlett window, see e.g. A. V. Oppenheim &
     R. W. Schafer, "Discrete-Time Signal Processing".

 - Function File:  blackman (M)
     Return the filter coefficients of a Blackman window of length M.

     For a definition of the  Blackman window, see e.g. A. V. Oppenheim
     & R. W. Schafer, "Discrete-Time Signal Processing".

 - Function File: [D, DD] = diffpara (X, A, B)
     Return the estimator D for the differencing parameter of an
     integrated time series.

     The frequencies from [2*pi*a/t, 2*pi*b/T] are used for the
     estimation.  If B is omitted, the interval [2*pi/T, 2*pi*a/T] is
     used.  If both B and A are omitted then a = 0.5 * sqrt (T) and b =
     1.5 * sqrt (T) is used, where T is the sample size.  If X is a
     matrix, the differencing parameter of each column is estimated.

     The estimators for all frequencies in the intervals described
     above is returned in DD.  The value of D is simply the mean of DD.

     Reference: Brockwell, Peter J. & Davis, Richard A. Time Series:
     Theory and Methods Springer 1987.

 - Function File:  durbinlevinson (C, OLDPHI, OLDV)
     Perform one step of the Durbin-Levinson algorithm.

     The vector C specifies the autocovariances `[gamma_0, ...,
     gamma_t]' from lag 0 to T, OLDPHI specifies the coefficients based
     on C(T-1) and OLDV specifies the corresponding error.

     If OLDPHI and OLDV are omitted, all steps from 1 to T of the
     algorithm are performed.

 - Function File:  fftshift (V)
     Perform a shift of the vector V, for use with the `fft' and `ifft'
     functions, in order the move the frequency 0 to the center of the
     vector or matrix.

     If V is a vector of E elements corresponding to E time samples
     spaced of Dt each, then `fftshift (fft (V))' corresponds to
     frequencies

          f = linspace (-E/(4*Dt), (E/2-1)/(2*Dt), E)

     If V is a matrix, the same holds for rows and columns.

 - Function File:  fractdiff (X, D)
     Compute the fractional differences (1-L)^d x where L denotes the
     lag-operator and d is greater than -1.

 - Function File:  hamming (M)
     Return the filter coefficients of a Hamming window of length M.

     For a definition of the Hamming window, see e.g. A. V. Oppenheim &
     R. W. Schafer, "Discrete-Time Signal Processing".

 - Function File:  hanning (M)
     Return the filter coefficients of a Hanning window of length M.

     For a definition of this window type, see e.g. A. V. Oppenheim &
     R. W. Schafer, "Discrete-Time Signal Processing".

 - Function File:  hurst (X)
     Estimate the Hurst parameter of sample X via the rescaled range
     statistic.  If X is a matrix, the parameter is estimated for every
     single column.

 - Function File:  periodogram (X)
     For a data matrix X from a sample of size N, return the
     periodogram.

 - Function File:  rectangle_lw (N, B)
     Rectangular lag window. Subfunction used for spectral density
     estimation.

 - Function File:  rectangle_sw (N, B)
     Rectangular spectral window.  Subfunction used for spectral density
     estimation.

 - Function File:  sinetone (FREQ, RATE, SEC, AMPL)
     Return a sinetone of frequency FREQ with length of SEC seconds at
     sampling rate RATE and with amplitude AMPL.  The arguments FREQ
     and AMPL may be vectors of common size.

     Defaults are RATE = 8000, SEC = 1 and AMPL = 64.

 - Function File:  sinewave (M, N, D)
     Return an M-element vector with I-th element given by `sin (2 * pi
     * (I+D-1) / N)'.

     The default value for D is 0.

 - Function File:  spectral_adf (C, WIN, B)
     Return the spectral density estimator given a vector of
     autocovariances C, window name WIN, and bandwidth, B.

     The window name, e.g., `"triangle"' or `"rectangle"' is used to
     search for a function called `WIN_sw'.

     If WIN is omitted, the triangle window is used.  If B is omitted,
     `1 / sqrt (length (X))' is used.

 - Function File:  spectral_xdf (X, WIN, B)
     Return the spectral density estimator given a data vector X,
     window name WIN, and bandwidth, B.

     The window name, e.g., `"triangle"' or `"rectangle"' is used to
     search for a function called `WIN_sw'.

     If WIN is omitted, the triangle window is used.  If B is omitted,
     `1 / sqrt (length (X))' is used.

 - Function File:  spencer (X)
     Return Spencer's 15 point moving average of every single column of
     X.

 - Function File: [Y, C] = stft (X, WIN_SIZE, INC, NUM_COEF, W_TYPE)
     Compute the short-term Fourier transform of the vector X with
     NUM_COEF coefficients by applying a window of WIN_SIZE data points
     and an increment of INC points.

     Before computing the Fourier transform, one of the following
     windows is applied:

    hanning
          w_type = 1

    hamming
          w_type = 2

    rectangle
          w_type = 3

     The window names can be passed as strings or by the W_TYPE number.

     If not all arguments are specified, the following defaults are
     used: WIN_SIZE = 80, INC = 24, NUM_COEF = 64, and W_TYPE = 1.

     `Y = stft (X, ...)' returns the absolute values of the Fourier
     coefficients according to the NUM_COEF positive frequencies.

     `[Y, C] = stft (`x', ...)' returns the entire STFT-matrix Y and a
     3-element vector C containing the window size, increment, and
     window type, which is needed by the synthesis function.

 - Function File:  synthesis (Y, C)
     Compute a signal from its short-time Fourier transform Y and a
     3-element vector C specifying window size, increment, and window
     type.

     The values Y and C can be derived by

          [Y, C] = stft (X , ...)

 - Function File:  triangle_lw (N, B)
     Triangular lag window. Subfunction used for spectral density
     estimation.

 - Function File:  triangle_sw (N, B)
     Triangular spectral window.  Subfunction used for spectral density
     estimation.

 - Function File: [A, V] = yulewalker (C)
     Fit an AR (p)-model with Yule-Walker estimates given a vector C of
     autocovariances `[gamma_0, ..., gamma_p]'.

     Returns the AR coefficients, A, and the variance of white noise, V.


Image Processing
****************

   Octave can display images with the X Window System using the
`xloadimage' program.  You do not need to be running X in order to
manipulate images, however, so some of these functions may be useful
even if you are not able to view the results.

   Loading images only works with Octave's image format (a file with a
matrix containing the image data, and a matrix containing the
colormap).  Contributions of robust, well-written functions to read
other image formats are welcome.  If you can provide them, or would like
to improve Octave's image processing capabilities in other ways, please
contact .

 - Function File:  colormap (MAP)
 - Function File:  colormap ("default")
     Set the current colormap.

     `colormap (MAP)' sets the current colormap to MAP.  The color map
     should be an N row by 3 column matrix.  The columns contain red,
     green, and blue intensities respectively.  All entries should be
     between 0 and 1 inclusive.  The new colormap is returned.

     `colormap ("default")' restores the default colormap (a gray scale
     colormap with 64 entries).  The default colormap is returned.

     With no arguments, `colormap' returns the current color map.

 - Function File:  gray (N)
     Return a gray colormap with N entries corresponding to values from
     0 to N-1.  The argument N should be a scalar.  If it is omitted,
     64 is assumed.

 - Function File: [IMG, MAP] = gray2ind ()
     Convert a gray scale intensity image to an Octave indexed image.

 - Function File:  image (X, ZOOM)
 - Function File:  image (X, Y, A, ZOOM)
     Display a matrix as a color image.  The elements of X are indices
     into the current colormap and should have values between 1 and the
     length of the colormap.  If ZOOM is omitted, the image will be
     scaled to fit within 600x350 (to a max of 4).

     It first tries to use `display' from `ImageMagick' then `xv' and
     then `xloadimage'.

     The axis values corresponding to the matrix elements are specified
     in X and Y. At present they are ignored.

 - Function File:  imagesc (A, ZOOM)
 - Function File:  imagesc (X, Y, A, ZOOM)
     Display a scaled version of the matrix A as a color image.  The
     matrix is scaled so that its entries are indices into the current
     colormap.  The scaled matrix is returned.  If ZOOM is omitted, a
     comfortable size is chosen.

     The axis values corresponding to the matrix elements are specified
     in X and Y.  At present they are ignored.

 - Function File:  imshow (X, MAP)
 - Function File:  imshow (X, N)
 - Function File:  imshow (I, N)
 - Function File:  imshow (R, G, B)
     Display images.

     `imshow (X)' displays an indexed image using the current colormap.

     `imshow (X, MAP)' displays an indexed image using the specified
     colormap.

     `imshow (I, N)' displays a gray scale intensity image.

     `imshow (R, G, B)' displays an RGB image.

 - Function File:  ind2gray (X, MAP)
     Convert an Octave indexed image to a gray scale intensity image.
     If MAP is omitted, the current colormap is used to determine the
     intensities.

 - Function File: [R, G, B] = ind2rgb (X, MAP)
     Convert an indexed image to red, green, and blue color components.
     If MAP is omitted, the current colormap is used for the conversion.

 - Function File: [X, MAP] = loadimage (FILE)
     Load an image file and it's associated color map from the specified
     FILE.  The image must be stored in Octave's image format.

 - Function File:  rgb2ntsc (RGB)
     Image format conversion.

 - Function File:  ntsc2rgb (YIQ)
     Image format conversion.

 - Function File:  ocean (N)
     Create color colormap.  The argument N should be a scalar.  If it
     is omitted, 64 is assumed.

 - Function File: [X, MAP] = rgb2ind (R, G, B)
     Convert and RGB image to an Octave indexed image.

 - Function File:  saveimage (FILE, X, FMT, MAP)
     Save the matrix X to FILE in image format FMT.  Valid values for
     FMT are

    `"img"'
          Octave's image format.  The current colormap is also saved in
          the file.

    `"ppm"'
          Portable pixmap format.

    `"ps"'
          PostScript format.  Note that images saved in PostScript
          format can not be read back into Octave with loadimage.

     If the fourth argument is supplied, the specified colormap will
     also be saved along with the image.

     Note: if the colormap contains only two entries and these entries
     are black and white, the bitmap ppm and PostScript formats are
     used.  If the image is a gray scale image (the entries within each
     row of the colormap are equal) the gray scale ppm and PostScript
     image formats are used, otherwise the full color formats are used.

 - Built-in Variable: IMAGEPATH
     A colon separated list of directories in which to search for image
     files.


Audio Processing
****************

   Octave provides a few functions for dealing with audio data.  An
audio `sample' is a single output value from an A/D converter, i.e., a
small integer number (usually 8 or 16 bits), and audio data is just a
series of such samples.  It can be characterized by three parameters:
the sampling rate (measured in samples per second or Hz, e.g. 8000 or
44100), the number of bits per sample (e.g. 8 or 16), and the number of
channels (1 for mono, 2 for stereo, etc.).

   There are many different formats for representing such data.
Currently, only the two most popular, _linear encoding_ and _mu-law
encoding_, are supported by Octave.  There is an excellent FAQ on audio
formats by Guido van Rossum  which can be found at any
FAQ ftp site, in particular in the directory
`/pub/usenet/news.answers/audio-fmts' of the archive site
`rtfm.mit.edu'.

   Octave simply treats audio data as vectors of samples (non-mono data
are not supported yet).  It is assumed that audio files using linear
encoding have one of the extensions `lin' or `raw', and that files
holding data in mu-law encoding end in `au', `mu', or `snd'.

 - Function File:  lin2mu (X)
     If the vector X represents mono audio data in 8- or 16-bit linear
     encoding, `lin2mu (X)' is the corresponding mu-law encoding.

 - Function File:  mu2lin (X, BPS)
     If the vector X represents mono audio data in mu-law encoding,
     `mu2lin' converts it to linear encoding.  The optional argument
     BPS specifies whether the input data uses 8 bit per sample
     (default) or 16 bit.

 - Function File:  loadaudio (NAME, EXT, BPS)
     Loads audio data from the file `NAME.EXT' into the vector X.

     The extension EXT determines how the data in the audio file is
     interpreted;  the extensions `lin' (default) and `raw' correspond
     to linear, the extensions `au', `mu', or `snd' to mu-law encoding.

     The argument BPS can be either 8 (default) or 16, and specifies
     the number of bits per sample used in the audio file.

 - Function File:  saveaudio (NAME, X, EXT, BPS)
     Saves a vector X of audio data to the file `NAME.EXT'.  The
     optional parameters EXT and BPS determine the encoding and the
     number of bits per sample used in the audio file (see
     `loadaudio');  defaults are `lin' and 8, respectively.

   The following functions for audio I/O require special A/D hardware
and operating system support.  It is assumed that audio data in linear
encoding can be played and recorded by reading from and writing to
`/dev/dsp', and that similarly `/dev/audio' is used for mu-law
encoding.  These file names are system-dependent.  Improvements so that
these functions will work without modification on a wide variety of
hardware are welcome.

 - Function File:  playaudio (NAME, EXT)
 - Function File:  playaudio (X)
     Plays the audio file `NAME.EXT' or the audio data stored in the
     vector X.

 - Function File:  record (SEC, SAMPLING_RATE)
     Records SEC seconds of audio input into the vector X.  The default
     value for SAMPLING_RATE is 8000 samples per second, or 8kHz.  The
     program waits until the user types  and then immediately
     starts to record.

 - Function File: setaudio ([W_TYPE [, VALUE]])
     executes the shell command `mixer [W_TYPE [, VALUE]]'


Quaternions
***********

   Quaternions are hypercomplex numbers used to represent spatial
rotations in three dimensions.  This set of routines provides a useful
basis for working with quaternions in Octave.  A tutorial is in the
Octave source, scripts/quaternion/quaternion.ps.

   These functions were written by A. S. Hodel, Associate Professor,
Auburn University.

 - Function File: [A, B, C, D] = quaternion (w)
 - Function File: [VV, THETA] = quaternion (w)
 - Function File: W = quaternion (A, B, C, D)
 - Function File: W = quaternion (VV, THETA)
     Construct or extract a quaternion

          w = a*i + b*j + c*k + d

     from given data.

 - Function File:  qconj (Q)
     Conjugate of a quaternion.

          q = [w, x, y, z] = w*i + x*j + y*k + z
          qconj (q) = -w*i -x*j -y*k + z

 - Function File:  qderiv (omega)
     Derivative of a quaternion.

     Let Q be a quaternion to transform a vector from a fixed frame to
     a rotating frame.  If the rotating frame is rotating about the [x,
     y, z] axes at angular rates [wx, wy, wz], then the derivative of Q
     is given by

          Q' = qderivmat (omega) * Q

     If the passive convention is used (rotate the frame, not the
     vector), then

          Q' = -qderivmat (omega) * Q

 - Function File:  qderivmat (OMEGA)
     Derivative of a quaternion.

     Let Q be a quaternion to transform a vector from a fixed frame to
     a rotating frame.  If the rotating frame is rotating about the [x,
     y, z] axes at angular rates [wx, wy, wz], then the derivative of Q
     is given by

          Q' = qderivmat (omega) * Q

     If the passive convention is used (rotate the frame, not the
     vector), then

          Q' = -qderivmat (omega) * Q.

 - Function File:  qinv (Q)
     Return the inverse of a quaternion.

          q = [w, x, y, z] = w*i + x*j + y*k + z
          qmult (q, qinv (q)) = 1 = [0 0 0 1]

 - Function File:  qmult (A, B)
     Multiply two quaternions.

          [w, x, y, z] = w*i + x*j + y*k + z

     identities:

          i^2 = j^2 = k^2 = -1
          ij = k                 jk = i
          ki = j                 kj = -i
          ji = -k                ik = -j

 - Function File:  qtrans (V, Q)
     Transform the unit quaternion V by the unit quaternion Q.  Returns
     `V = Q*V/Q'.

 - Function File:  qtransv (V, Q)
     Transform the 3-D vector V by the unit quaternion Q.  Return a
     column vector.

          vi = (2*real(q)^2 - 1)*vb + 2*imag(q)*(imag(q)'*vb)
             + 2*real(q)*cross(imag(q),vb)

     Where imag(q) is a column vector of length 3.

 - Function File:  qtransvmat (QIB)
     Construct a 3x3 transformation matrix from quaternion QIB that is
     equivalent to rotation of th radians about axis VV, where `[VV,
     TH] = quaternion (QIB)'.

 - Function File:  qcoordinate_plot (QF, QB, QV)
     Plot in the current figure a set of coordinate axes as viewed from
     the orientation specified by quaternion QV.  Inertial axes are
     also plotted:

    QF
          Quaternion from reference (x,y,z) to inertial.

    QB
          Quaternion from reference to body.

    QV
          Quaternion from reference to view angle.