# Created by Octave 3.6.1, Thu May 17 20:38:50 2012 UTC <root@brouzouf>
# name: cache
# type: cell
# rows: 3
# columns: 130
# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
ar_psd


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 4320
 Usage:
   [psd,f_out] = ar_psd(a,v,freq,Fs,range,method,plot_type)

  Calculate the power spectrum of the autoregressive model

                         M
  x(n) = sqrt(v).e(n) + SUM a(k).x(n-k)
                        k=1
  where x(n) is the output of the model and e(n) is white noise.
  This function is intended for use with 
    [a,v,k] = arburg(x,poles,criterion)
  which use the Burg (1968) method to calculate a "maximum entropy"
  autoregressive model of "x".  This function runs on octave and matlab.
  
  If the "freq" argument is a vector (of frequencies) the spectrum is
  calculated using the polynomial method and the "method" argument is
  ignored.  For scalar "freq", an integer power of 2, or "method='FFT'",
  causes the spectrum to be calculated by FFT.  Otherwise, the spectrum
  is calculated as a polynomial.  It may be computationally more
  efficient to use the FFT method if length of the model is not much
  smaller than the number of frequency values. The spectrum is scaled so
  that spectral energy (area under spectrum) is the same as the
  time-domain energy (mean square of the signal).

 ARGUMENTS:
     All but the first two arguments are optional and may be empty.

   a      %% [vector] list of M=(order+1) autoregressive model
          %%      coefficients.  The first element of "ar_coeffs" is the
          %%      zero-lag coefficient, which always has a value of 1.

   v      %% [real scalar] square of the moving-average coefficient of
          %%               the AR model.

   freq   %% [real vector] frequencies at which power spectral density
          %%               is calculated
          %% [integer scalar] number of uniformly distributed frequency
          %%          values at which spectral density is calculated.
          %%          [default=256]

   Fs     %% [real scalar] sampling frequency (Hertz) [default=1]

 CONTROL-STRING ARGUMENTS -- each of these arguments is a character string.
   Control-string arguments can be in any order after the other arguments.

   range  %% 'half',  'onesided' : frequency range of the spectrum is
          %%       from zero up to but not including sample_f/2.  Power
          %%       from negative frequencies is added to the positive
          %%       side of the spectrum.
          %% 'whole', 'twosided' : frequency range of the spectrum is
          %%       -sample_f/2 to sample_f/2, with negative frequencies
          %%       stored in "wrap around" order after the positive
          %%       frequencies; e.g. frequencies for a 10-point 'twosided'
          %%       spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1
          %% 'shift', 'centerdc' : same as 'whole' but with the first half
          %%       of the spectrum swapped with second half to put the
          %%       zero-frequency value in the middle. (See "help
          %%       fftshift". If "freq" is vector, 'shift' is ignored.
          %% If model coefficients "ar_coeffs" are real, the default
          %% range is 'half', otherwise default range is 'whole'.

   method %% 'fft':  use FFT to calculate power spectrum.
          %% 'poly': calculate power spectrum as a polynomial of 1/z
          %% N.B. this argument is ignored if the "freq" argument is a
          %%      vector.  The default is 'poly' unless the "freq"
          %%      argument is an integer power of 2.
   
 plot_type%% 'plot', 'semilogx', 'semilogy', 'loglog', 'squared' or 'db':
          %% specifies the type of plot.  The default is 'plot', which
          %% means linear-linear axes. 'squared' is the same as 'plot'.
          %% 'dB' plots "10*log10(psd)".  This argument is ignored and a
          %% spectrum is not plotted if the caller requires a returned
          %% value.

 RETURNED VALUES:
     If returned values are not required by the caller, the spectrum
     is plotted and nothing is returned.
   psd    %% [real vector] estimate of power-spectral density
   f_out  %% [real vector] frequency values 

 N.B. arburg runs in octave and matlab, and does not depend on octave-forge
      or signal-processing-toolbox functions.

 REFERENCE
 [1] Equation 2.28 from Steven M. Kay and Stanley Lawrence Marple Jr.:
   "Spectrum analysis -- a modern perspective",
   Proceedings of the IEEE, Vol 69, pp 1380-1419, Nov., 1981




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 68
 Usage:
   [psd,f_out] = ar_psd(a,v,freq,Fs,range,method,plot_type)




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
arburg


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 4080
 [a,v,k] = arburg(x,poles,criterion)

 Calculate coefficients of an autoregressive (AR) model of complex data
 "x" using the whitening lattice-filter method of Burg (1968).  The inverse
 of the model is a moving-average filter which reduces "x" to white noise.
 The power spectrum of the AR model is an estimate of the maximum
 entropy power spectrum of the data.  The function "ar_psd" calculates the
 power spectrum of the AR model.

 ARGUMENTS:
   x         %% [vector] sampled data

   poles     %% [integer scalar] number of poles in the AR model or
             %%       limit to the number of poles if a
             %%       valid "stop_crit" is provided.

   criterion %% [optional string arg]  model-selection criterion.  Limits
             %%       the number of poles so that spurious poles are not 
             %%       added when the whitened data has no more information
             %%       in it (see Kay & Marple, 1981). Recognised values are
             %%  'AKICc' -- approximate corrected Kullback information
             %%             criterion (recommended),
             %%   'KIC'  -- Kullback information criterion
             %%   'AICc' -- corrected Akaike information criterion
             %%   'AIC'  -- Akaike information criterion
             %%   'FPE'  -- final prediction error" criterion
             %% The default is to NOT use a model-selection criterion

 RETURNED VALUES:
   a         %% [polynomial/vector] list of (P+1) autoregression coeffic-
             %%       ients; for data input x(n) and  white noise e(n),
             %%       the model is
             %%                             P+1
             %%       x(n) = sqrt(v).e(n) + SUM a(k).x(n-k)
             %%                             k=1

   v         %% [real scalar] mean square of residual noise from the
             %%       whitening operation of the Burg lattice filter.

   k         %% [column vector] reflection coefficients defining the
             %%       lattice-filter embodiment of the model

 HINTS:
  (1) arburg does not remove the mean from the data.  You should remove
      the mean from the data if you want a power spectrum.  A non-zero mean
      can produce large errors in a power-spectrum estimate.  See
      "help detrend".
  (2) If you don't know what the value of "poles" should be, choose the
      largest (reasonable) value you could want and use the recommended
      value, criterion='AKICc', so that arburg can find it.
      E.g. arburg(x,64,'AKICc')
      The AKICc has the least bias and best resolution of the available
      model-selection criteria.
  (3) arburg runs in octave and matlab, does not depend on octave forge
      or signal-processing-toolbox functions.
  (4) Autoregressive and moving-average filters are stored as polynomials
      which, in matlab, are row vectors.

 NOTE ON SELECTION CRITERION
   AIC, AICc, KIC and AKICc are based on information theory.  They  attempt
   to balance the complexity (or length) of the model against how well the
   model fits the data.  AIC and KIC are biassed estimates of the asymmetric
   and the symmetric Kullback-Leibler divergence respectively.  AICc and
   AKICc attempt to correct the bias. See reference [4].


 REFERENCES
 [1] John Parker Burg (1968)
   "A new analysis technique for time series data",
   NATO advanced study Institute on Signal Processing with Emphasis on
   Underwater Acoustics, Enschede, Netherlands, Aug. 12-23, 1968.

 [2] Steven M. Kay and Stanley Lawrence Marple Jr.:
   "Spectrum analysis -- a modern perspective",
   Proceedings of the IEEE, Vol 69, pp 1380-1419, Nov., 1981

 [3] William H. Press and Saul A. Teukolsky and William T. Vetterling and
               Brian P. Flannery
   "Numerical recipes in C, The art of scientific computing", 2nd edition,
   Cambridge University Press, 2002 --- Section 13.7.

 [4] Abd-Krim Seghouane and Maiza Bekara
   "A small sample model selection criterion based on Kullback's symmetric
   divergence", IEEE Transactions on Signal Processing,
   Vol. 52(12), pp 3314-3323, Dec. 2004



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 37
 [a,v,k] = arburg(x,poles,criterion)




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
aryule


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 799
 usage:  [a, v, k] = aryule (x, p)
 
 fits an AR (p)-model with Yule-Walker estimates.
 x = data vector to estimate
 a: AR coefficients
 v: variance of white noise
 k: reflection coeffients for use in lattice filter 

 The power spectrum of the resulting filter can be plotted with
 pyulear(x, p), or you can plot it directly with ar_psd(a,v,...).

 See also:
 pyulear, power, freqz, impz -- for observing characteristics of the model
 arburg -- for alternative spectral estimators

 Example: Use example from arburg, but substitute aryule for arburg.

 Note: Orphanidis '85 claims lattice filters are more tolerant of 
 truncation errors, which is why you might want to use them.  However,
 lacking a lattice filter processor, I haven't tested that the lattice
 filter coefficients are reasonable.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 usage:  [a, v, k] = aryule (x, p)
 
 fits an AR (p)-model with Yule-Walker esti



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 11
barthannwin


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 136
 -- Function File: [W] = barthannwin(L)
     Compute the modified Bartlett-Hann window of lenght L.

     See also: rectwin, bartlett





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 54
Compute the modified Bartlett-Hann window of lenght L.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
besselap


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 105
 Return bessel analog filter prototype.

 References: 

 http://en.wikipedia.org/wiki/Bessel_polynomials



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 39
 Return bessel analog filter prototype.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 7
besself


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 804
 Generate a bessel filter.
 Default is a Laplace space (s) filter.
 
 [b,a] = besself(n, Wc)
    low pass filter with cutoff pi*Wc radians

 [b,a] = besself(n, Wc, 'high')
    high pass filter with cutoff pi*Wc radians

 [b,a] = besself(n, [Wl, Wh])
    band pass filter with edges pi*Wl and pi*Wh radians

 [b,a] = besself(n, [Wl, Wh], 'stop')
    band reject filter with edges pi*Wl and pi*Wh radians

 [z,p,g] = besself(...)
    return filter as zero-pole-gain rather than coefficients of the
    numerator and denominator polynomials.
 
 [...] = besself(...,'z')
     return a discrete space (Z) filter, W must be less than 1.
 
 [a,b,c,d] = besself(...)
  return  state-space matrices 

 References: 

 Proakis & Manolakis (1992). Digital Signal Processing. New York:
 Macmillan Publishing Company.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 26
 Generate a bessel filter.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
bilinear


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 2658
 usage: [Zz, Zp, Zg] = bilinear(Sz, Sp, Sg, T)
        [Zb, Za] = bilinear(Sb, Sa, T)

 Transform a s-plane filter specification into a z-plane
 specification. Filters can be specified in either zero-pole-gain or
 transfer function form. The input form does not have to match the
 output form. 1/T is the sampling frequency represented in the z plane.

 Note: this differs from the bilinear function in the signal processing
 toolbox, which uses 1/T rather than T.

 Theory: Given a piecewise flat filter design, you can transform it
 from the s-plane to the z-plane while maintaining the band edges by
 means of the bilinear transform.  This maps the left hand side of the
 s-plane into the interior of the unit circle.  The mapping is highly
 non-linear, so you must design your filter with band edges in the
 s-plane positioned at 2/T tan(w*T/2) so that they will be positioned
 at w after the bilinear transform is complete.

 The following table summarizes the transformation:

 +---------------+-----------------------+----------------------+
 | Transform     | Zero at x             | Pole at x            |
 |    H(S)       |   H(S) = S-x          |    H(S)=1/(S-x)      |
 +---------------+-----------------------+----------------------+
 |       2 z-1   | zero: (2+xT)/(2-xT)   | zero: -1             |
 |  S -> - ---   | pole: -1              | pole: (2+xT)/(2-xT)  |
 |       T z+1   | gain: (2-xT)/T        | gain: (2-xT)/T       |
 +---------------+-----------------------+----------------------+

 With tedious algebra, you can derive the above formulae yourself by
 substituting the transform for S into H(S)=S-x for a zero at x or
 H(S)=1/(S-x) for a pole at x, and converting the result into the
 form:

    H(Z)=g prod(Z-Xi)/prod(Z-Xj)

 Please note that a pole and a zero at the same place exactly cancel.
 This is significant since the bilinear transform creates numerous
 extra poles and zeros, most of which cancel. Those which do not
 cancel have a "fill-in" effect, extending the shorter of the sets to
 have the same number of as the longer of the sets of poles and zeros
 (or at least split the difference in the case of the band pass
 filter). There may be other opportunistic cancellations but I will
 not check for them.

 Also note that any pole on the unit circle or beyond will result in
 an unstable filter.  Because of cancellation, this will only happen
 if the number of poles is smaller than the number of zeros.  The
 analytic design methods all yield more poles than zeros, so this will
 not be a problem.

 References:

 Proakis & Manolakis (1992). Digital Signal Processing. New York:
 Macmillan Publishing Company.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 usage: [Zz, Zp, Zg] = bilinear(Sz, Sp, Sg, T)
        [Zb, Za] = bilinear(Sb, S



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 11
bitrevorder


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 111
 -- Function File: [Y I] = bitrevorder(X)
     Reorder x in the bit reversed order

     See also: fft, ifft





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 36
Reorder x in the bit reversed order




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 14
blackmanharris


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 120
 -- Function File: [W] = blackmanharris(L)
     Compute the Blackman-Harris window.

     See also: rectwin, bartlett





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 35
Compute the Blackman-Harris window.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 15
blackmannuttall


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 123
 -- Function File: [W] = blackmannuttall(L)
     Compute the Blackman-Nuttall window.

     See also: nuttallwin, kaiser





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 36
Compute the Blackman-Nuttall window.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 9
bohmanwin


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 118
 -- Function File: [W] = bohmanwin(L)
     Compute the Bohman window of lenght L.

     See also: rectwin, bartlett





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 38
Compute the Bohman window of lenght L.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
boxcar


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 95
 usage:  w = boxcar (n)

 Returns the filter coefficients of a rectangular window of length n.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 24
 usage:  w = boxcar (n)




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
buffer


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1473
 -- Function File: Y =  buffer (X, N, P, OPT)
 -- Function File: [Y, Z, OPT] =  buffer (...)
     Buffer a signal into a data frame. The arguments to `buffer' are

    X
          The data to be buffered.

    N
          The number of rows in the produced data buffer. This is an
          positive integer value and must be supplied.

    P
          An integer less than N that specifies the under- or overlap
          between column in the data frame. The default value of P is 0.

    OPT
          In the case of an overlap, OPT can be either a vector of
          length P or the string 'nodelay'. If OPT is a vector, then the
          first P entries in Y will be filled with these values. If OPT
          is the string 'nodelay', then the first value of Y
          corresponds to the first value of X.

          In the can of an underlap, OPT must be an integer between 0
          and `-P'. The represents the initial underlap of the first
          column of Y.

          The default value for OPT the vector `zeros (1, P)' in the
          case of an overlap, or 0 otherwise.

     In the case of a single output argument, Y will be padded with
     zeros to fill the missing values in the data frame. With two output
     arguments Z is the remaining data that has not been used in the
     current data frame.

     Likewise, the output OPT is the overlap, or underlap that might be
     used for a future call to `code' to allow continuous buffering.




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 34
Buffer a signal into a data frame.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
butter


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 799
 Generate a butterworth filter.
 Default is a discrete space (Z) filter.
 
 [b,a] = butter(n, Wc)
    low pass filter with cutoff pi*Wc radians

 [b,a] = butter(n, Wc, 'high')
    high pass filter with cutoff pi*Wc radians

 [b,a] = butter(n, [Wl, Wh])
    band pass filter with edges pi*Wl and pi*Wh radians

 [b,a] = butter(n, [Wl, Wh], 'stop')
    band reject filter with edges pi*Wl and pi*Wh radians

 [z,p,g] = butter(...)
    return filter as zero-pole-gain rather than coefficients of the
    numerator and denominator polynomials.
 
 [...] = butter(...,'s')
     return a Laplace space filter, W can be larger than 1.
 
 [a,b,c,d] = butter(...)
  return  state-space matrices 

 References: 

 Proakis & Manolakis (1992). Digital Signal Processing. New York:
 Macmillan Publishing Company.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 31
 Generate a butterworth filter.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 7
buttord


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1107
 Compute butterworth filter order and cutoff for the desired response
 characteristics. Rp is the allowable decibels of ripple in the pass 
 band. Rs is the minimum attenuation in the stop band.

 [n, Wc] = buttord(Wp, Ws, Rp, Rs)
     Low pass (Wp<Ws) or high pass (Wp>Ws) filter design.  Wp is the
     pass band edge and Ws is the stop band edge.  Frequencies are
     normalized to [0,1], corresponding to the range [0,Fs/2].
 
 [n, Wc] = buttord([Wp1, Wp2], [Ws1, Ws2], Rp, Rs)
     Band pass (Ws1<Wp1<Wp2<Ws2) or band reject (Wp1<Ws1<Ws2<Wp2)
     filter design. Wp gives the edges of the pass band, and Ws gives
     the edges of the stop band.

 Theory: |H(W)|^2 = 1/[1+(W/Wc)^(2N)] = 10^(-R/10)
 With some algebra, you can solve simultaneously for Wc and N given
 Ws,Rs and Wp,Rp.  For high pass filters, subtracting the band edges
 from Fs/2, performing the test, and swapping the resulting Wc back
 works beautifully.  For bandpass and bandstop filters this process
 significantly overdesigns.  Artificially dividing N by 2 in this case
 helps a lot, but it still overdesigns.

 See also: butter



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 Compute butterworth filter order and cutoff for the desired response
 character



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 5
cceps


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 195
 usage:  cceps (x [, correct])

 Returns the complex cepstrum of the vector x.
 If the optional argument correct has the value 1, a correction
 method is applied.  The default is not to do this.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 31
 usage:  cceps (x [, correct])




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 4
cheb


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 371
 Usage:  cheb (n, x)

 Returns the value of the nth-order Chebyshev polynomial calculated at
 the point x. The Chebyshev polynomials are defined by the equations:

           / cos(n acos(x),    |x| <= 1
   Tn(x) = |
           \ cosh(n acosh(x),  |x| > 1

 If x is a vector, the output is a vector of the same size, where each
 element is calculated as y(i) = Tn(x(i)).



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 21
 Usage:  cheb (n, x)




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
cheb1ord


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 678
 Compute chebyshev type I filter order and cutoff for the desired response
 characteristics. Rp is the allowable decibels of ripple in the pass 
 band. Rs is the minimum attenuation in the stop band.

 [n, Wc] = cheb1ord(Wp, Ws, Rp, Rs)
     Low pass (Wp<Ws) or high pass (Wp>Ws) filter design.  Wp is the
     pass band edge and Ws is the stop band edge.  Frequencies are
     normalized to [0,1], corresponding to the range [0,Fs/2].
 
 [n, Wc] = cheb1ord([Wp1, Wp2], [Ws1, Ws2], Rp, Rs)
     Band pass (Ws1<Wp1<Wp2<Ws2) or band reject (Wp1<Ws1<Ws2<Wp2)
     filter design. Wp gives the edges of the pass band, and Ws gives
     the edges of the stop band.

 See also: cheby1



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 Compute chebyshev type I filter order and cutoff for the desired response
 char



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
cheb2ord


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 690
 Compute chebyshev type II filter order and cutoff for the desired response
 characteristics. Rp is the allowable decibels of ripple in the pass 
 band. Rs is the minimum attenuation in the stop band.

 [n, Wc] = cheb2ord(Wp, Ws, Rp, Rs)
     Low pass (Wp<Ws) or high pass (Wp>Ws) filter design.  Wp is the
     pass band edge and Ws is the stop band edge.  Frequencies are
     normalized to [0,1], corresponding to the range [0,Fs/2].
 
 [n, Wc] = cheb2ord([Wp1, Wp2], [Ws1, Ws2], Rp, Rs)
     Band pass (Ws1<Wp1<Wp2<Ws2) or band reject (Wp1<Ws1<Ws2<Wp2)
     filter design. Wp gives the edges of the pass band, and Ws gives
     the edges of the stop band.

 Theory: 

 See also: cheby2



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 Compute chebyshev type II filter order and cutoff for the desired response
 cha



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 7
chebwin


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1095
 Usage:  chebwin (L, at)

 Returns the filter coefficients of the L-point Dolph-Chebyshev window
 with at dB of attenuation in the stop-band of the corresponding
 Fourier transform.

 For the definition of the Chebyshev window, see

 * Peter Lynch, "The Dolph-Chebyshev Window: A Simple Optimal Filter",
   Monthly Weather Review, Vol. 125, pp. 655-660, April 1997.
   (http://www.maths.tcd.ie/~plynch/Publications/Dolph.pdf)

 * C. Dolph, "A current distribution for broadside arrays which
   optimizes the relationship between beam width and side-lobe level",
   Proc. IEEE, 34, pp. 335-348.

 The window is described in frequency domain by the expression:

          Cheb(L-1, beta * cos(pi * k/L))
   W(k) = -------------------------------
                 Cheb(L-1, beta)

 with

   beta = cosh(1/(L-1) * acosh(10^(at/20))

 and Cheb(m,x) denoting the m-th order Chebyshev polynomial calculated
 at the point x.

 Note that the denominator in W(k) above is not computed, and after
 the inverse Fourier transform the window is scaled by making its
 maximum value unitary.

 See also: kaiser



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 25
 Usage:  chebwin (L, at)




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
cheby1


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 802
 Generate an Chebyshev type I filter with Rp dB of pass band ripple.
 
 [b, a] = cheby1(n, Rp, Wc)
    low pass filter with cutoff pi*Wc radians

 [b, a] = cheby1(n, Rp, Wc, 'high')
    high pass filter with cutoff pi*Wc radians

 [b, a] = cheby1(n, Rp, [Wl, Wh])
    band pass filter with edges pi*Wl and pi*Wh radians

 [b, a] = cheby1(n, Rp, [Wl, Wh], 'stop')
    band reject filter with edges pi*Wl and pi*Wh radians

 [z, p, g] = cheby1(...)
    return filter as zero-pole-gain rather than coefficients of the
    numerator and denominator polynomials.

 [...] = cheby1(...,'s')
     return a Laplace space filter, W can be larger than 1.
 
 [a,b,c,d] = cheby1(...)
  return  state-space matrices 
 
 References: 

 Parks & Burrus (1987). Digital Filter Design. New York:
 John Wiley & Sons, Inc.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 68
 Generate an Chebyshev type I filter with Rp dB of pass band ripple.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
cheby2


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 808
 Generate an Chebyshev type II filter with Rs dB of stop band attenuation.
 
 [b, a] = cheby2(n, Rs, Wc)
    low pass filter with cutoff pi*Wc radians

 [b, a] = cheby2(n, Rs, Wc, 'high')
    high pass filter with cutoff pi*Wc radians

 [b, a] = cheby2(n, Rs, [Wl, Wh])
    band pass filter with edges pi*Wl and pi*Wh radians

 [b, a] = cheby2(n, Rs, [Wl, Wh], 'stop')
    band reject filter with edges pi*Wl and pi*Wh radians

 [z, p, g] = cheby2(...)
    return filter as zero-pole-gain rather than coefficients of the
    numerator and denominator polynomials.

 [...] = cheby2(...,'s')
     return a Laplace space filter, W can be larger than 1.
 
 [a,b,c,d] = cheby2(...)
  return  state-space matrices 
 
 References: 

 Parks & Burrus (1987). Digital Filter Design. New York:
 John Wiley & Sons, Inc.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 74
 Generate an Chebyshev type II filter with Rs dB of stop band attenuation.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 5
chirp


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 811
 usage: y = chirp(t [, f0 [, t1 [, f1 [, form [, phase]]]]])

 Evaluate a chirp signal at time t.  A chirp signal is a frequency
 swept cosine wave.

 t: vector of times to evaluate the chirp signal
 f0: frequency at time t=0 [ 0 Hz ]
 t1: time t1 [ 1 sec ]
 f1: frequency at time t=t1 [ 100 Hz ]
 form: shape of frequency sweep
    'linear'      f(t) = (f1-f0)*(t/t1) + f0
    'quadratic'   f(t) = (f1-f0)*(t/t1)^2 + f0
    'logarithmic' f(t) = (f1-f0)^(t/t1) + f0
 phase: phase shift at t=0

 Example
    specgram(chirp([0:0.001:5])); # linear, 0-100Hz in 1 sec
    specgram(chirp([-2:0.001:15], 400, 10, 100, 'quadratic'));
    soundsc(chirp([0:1/8000:5], 200, 2, 500, "logarithmic"),8000);

 If you want a different sweep shape f(t), use the following:
    y = cos(2*pi*integral(f(t)) + 2*pi*f0*t + phase);



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 61
 usage: y = chirp(t [, f0 [, t1 [, f1 [, form [, phase]]]]])




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
cmorwavf


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 96
 -- Function File: [PSI,X] = cmorwavf (LB,UB,N,FB,FC)
     Compute the Complex Morlet wavelet.




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 35
Compute the Complex Morlet wavelet.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
cohere


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 377
 Usage:
   [Pxx,freq] = cohere(x,y,Nfft,Fs,window,overlap,range,plot_type,detrend)

     Estimate (mean square) coherence of signals "x" and "y".
     Use the Welch (1967) periodogram/FFT method.
     Compatible with Matlab R11 cohere and earlier.
     See "help pwelch" for description of arguments, hints and references
     --- especially hint (7) for Matlab R11 defaults.




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 Usage:
   [Pxx,freq] = cohere(x,y,Nfft,Fs,window,overlap,range,plot_type,detren



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 7
convmtx


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 498
 -- Function File:  convmtx (A, N)
     If A is a column vector and X is a column vector of length N, then

     `convmtx(A, N) * X'

     gives the convolution of of A and X and is the same as `conv(A,
     X)'. The difference is if many vectors are to be convolved with
     the same vector, then this technique is possibly faster.

     Similarly, if A is a row vector and X is a row vector of length N,
     then

     `X * convmtx(A, N)'

     is the same as `conv(X, A)'.

   See also: conv




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 67
If A is a column vector and X is a column vector of length N, then




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
cplxreal


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 710
 -- Function File: [ZC, ZR] = cplxreal (Z, THRESH)
     Split the vector z into its complex (ZC) and real (ZR) elements,
     eliminating one of each complex-conjugate pair.

     INPUTS:
        *   Z      = row- or column-vector of complex numbers
        *   THRESH = tolerance threshold for numerical comparisons
          (default = 100*eps)

     RETURNED:
        * ZC = elements of Z having positive imaginary parts
        * ZR = elements of Z having zero imaginary part

     Each complex element of Z is assumed to have a complex-conjugate
     counterpart elsewhere in Z as well.  Elements are declared real if
     their imaginary parts have magnitude less than THRESH.

     See also: cplxpair





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
Split the vector z into its complex (ZC) and real (ZR) elements,
eliminating one



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 4
cpsd


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 260
 Usage:
   [Pxx,freq] = cpsd(x,y,Nfft,Fs,window,overlap,range,plot_type,detrend)

     Estimate cross power spectrum of data "x" and "y" by the Welch (1967)
     periodogram/FFT method.
     See "help pwelch" for description of arguments, hints and references



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 Usage:
   [Pxx,freq] = cpsd(x,y,Nfft,Fs,window,overlap,range,plot_type,detrend)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 3
csd


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 358
 Usage:
   [Pxx,freq] = csd(x,y,Nfft,Fs,window,overlap,range,plot_type,detrend)

     Estimate cross power spectrum of data "x" and "y" by the Welch (1967)
     periodogram/FFT method.  Compatible with Matlab R11 csd and earlier.
     See "help pwelch" for description of arguments, hints and references
     --- especially hint (7) for Matlab R11 defaults.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 Usage:
   [Pxx,freq] = csd(x,y,Nfft,Fs,window,overlap,range,plot_type,detrend)




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 3
czt


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 828
 usage y=czt(x, m, w, a)

 Chirp z-transform.  Compute the frequency response starting at a and
 stepping by w for m steps.  a is a point in the complex plane, and
 w is the ratio between points in each step (i.e., radius increases
 exponentially, and angle increases linearly).

 To evaluate the frequency response for the range f1 to f2 in a signal
 with sampling frequency Fs, use the following:
     m = 32;                               ## number of points desired
     w = exp(-j*2*pi*(f2-f1)/((m-1)*Fs));  ## freq. step of f2-f1/m
     a = exp(j*2*pi*f1/Fs);                ## starting at frequency f1
     y = czt(x, m, w, a);

 If you don't specify them, then the parameters default to a fourier 
 transform:
     m=length(x), w=exp(-j*2*pi/m), a=1

 If x is a matrix, the transform will be performed column-by-column.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 25
 usage y=czt(x, m, w, a)




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
data2fun


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1111
 -- Function File: [FHANDLE, FULLNAME] =  data2fun (TI, YI)
 -- Function File: [ ... ] =  data2fun (TI, YI,PROPERTY,VALUE)
     Creates a vectorized function based on data samples using
     interpolation.

     The values given in YI (N-by-k matrix) correspond to evaluations
     of the function y(t) at the points TI (N-by-1 matrix).  The data
     is interpolated and the function handle to the generated
     interpolant is returned.

     The function accepts property-value pairs described below.

    `file'
          Code is generated and .m file is created. The VALUE contains
          the name of the function. The returned function handle is a
          handle to that file. If VALUE is empty, then a name is
          automatically generated using `tmpnam' and the file is
          created in the current directory. VALUE must not have an
          extension, since .m will be appended.  Numerical value used
          in the function are stored in a .mat file with the same name
          as the function.

    `interp'
          Type of interpolation. See `interp1'.


     See also: interp1





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 72
Creates a vectorized function based on data samples using interpolation.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 3
dct


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 687
 y = dct (x, n)
    Computes the discrete cosine transform of x.  If n is given, then
    x is padded or trimmed to length n before computing the transform.
    If x is a matrix, compute the transform along the columns of the
    the matrix. The transform is faster if x is real-valued and even
    length.

 The discrete cosine transform X of x can be defined as follows:

               N-1
   X[k] = w(k) sum x[n] cos (pi (2n+1) k / 2N ),  k = 0, ..., N-1
               n=0

 with w(0) = sqrt(1/N) and w(k) = sqrt(2/N), k = 1, ..., N-1.  There
 are other definitions with different scaling of X[k], but this form
 is common in image processing.

 See also: idct, dct2, idct2, dctmtx



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 64
 y = dct (x, n)
    Computes the discrete cosine transform of x.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 4
dct2


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 202
 y = dct2 (x)
   Computes the 2-D discrete cosine transform of matrix x

 y = dct2 (x, m, n) or y = dct2 (x, [m n])
   Computes the 2-D DCT of x after padding or trimming rows to m and
   columns to n.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 72
 y = dct2 (x)
   Computes the 2-D discrete cosine transform of matrix x




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
dctmtx


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 599
 T = dctmtx (n)
 Return the DCT transformation matrix of size n x n.

 If A is an n x n matrix, then the following are true:
     T*A    == dct(A),  T'*A   == idct(A)
     T*A*T' == dct2(A), T'*A*T == idct2(A)

 A dct transformation matrix is useful for doing things like jpeg
 image compression, in which an 8x8 dct matrix is applied to
 non-overlapping blocks throughout an image and only a subblock on the
 top left of each block is kept.  During restoration, the remainder of
 the block is filled with zeros and the inverse transform is applied
 to the block.

 See also: dct, idct, dct2, idct2



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 68
 T = dctmtx (n)
 Return the DCT transformation matrix of size n x n.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
decimate


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1021
 usage: y = decimate(x, q [, n] [, ftype])

 Downsample the signal x by a factor of q, using an order n filter
 of ftype 'fir' or 'iir'.  By default, an order 8 Chebyshev type I
 filter is used or a 30 point FIR filter if ftype is 'fir'.  Note
 that q must be an integer for this rate change method.

 Example
    ## Generate a signal that starts away from zero, is slowly varying
    ## at the start and quickly varying at the end, decimate and plot.
    ## Since it starts away from zero, you will see the boundary
    ## effects of the antialiasing filter clearly.  Next you will see
    ## how it follows the curve nicely in the slowly varying early
    ## part of the signal, but averages the curve in the quickly
    ## varying late part of the signal.
    t=0:0.01:2; x=chirp(t,2,.5,10,'quadratic')+sin(2*pi*t*0.4); 
    y = decimate(x,4);   # factor of 4 decimation
    stem(t(1:121)*1000,x(1:121),"-g;Original;"); hold on; # plot original
    stem(t(1:4:121)*1000,y(1:31),"-r;Decimated;"); hold off; # decimated



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 43
 usage: y = decimate(x, q [, n] [, ftype])




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
dftmtx


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 343
 -- Function File: D =  dftmtx (N)
     If N is a scalar, produces a N-by-N matrix D such that the Fourier
     transform of a column vector of length N is given by `dftmtx(N) *
     x' and the inverse Fourier transform is given by `inv(dftmtx(N)) *
     x'. In general this is less efficient than calling the "fft" and
     "ifft" directly.




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
If N is a scalar, produces a N-by-N matrix D such that the Fourier
transform of 



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 5
diric


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 116
 -- Function File: [Y] = diric(X,N)
     Compute the dirichlet function.

     See also: sinc, gauspuls, sawtooth





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 31
Compute the dirichlet function.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 10
downsample


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 511
 -- Function File: Y = downsample (X, N)
 -- Function File: Y = downsample (X, N, OFFSET)
     Downsample the signal, selecting every nth element.  If X is a
     matrix, downsample every column.

     For most signals you will want to use `decimate' instead since it
     prefilters the high frequency components of the signal and avoids
     aliasing effects.

     If OFFSET is defined, select every nth element starting at sample
     OFFSET.

     See also: decimate, interp, resample, upfirdn, upsample





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 51
Downsample the signal, selecting every nth element.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 3
dst


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 482
 -- Function File: Y = dst (X)
 -- Function File: Y = dst (X, N)
     Computes the type I discrete sine transform of X.  If N is given,
     then X is padded or trimmed to length N before computing the
     transform.  If X is a matrix, compute the transform along the
     columns of the the matrix.

     The discrete sine transform X of x can be defined as follows:

            N
     X[k] = sum x[n] sin (pi n k / (N+1) ),  k = 1, ..., N
            n=1

     See also: idst





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 49
Computes the type I discrete sine transform of X.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 3
dwt


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 111
 -- Function File: [CA CD] = dwt(X,LO_D,HI_D)
     Comupte de discrete wavelet transform of x with one level.




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 58
Comupte de discrete wavelet transform of x with one level.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 5
ellip


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1050
 N-ellip 0.2.1
usage: [Zz, Zp, Zg] = ellip(n, Rp, Rs, Wp, stype,'s')

 Generate an Elliptic or Cauer filter (discrete and contnuious).
 
 [b,a] = ellip(n, Rp, Rs, Wp)
  low pass filter with order n, cutoff pi*Wp radians, Rp decibels 
  of ripple in the passband and a stopband Rs decibels down.

 [b,a] = ellip(n, Rp, Rs, Wp, 'high')
  high pass filter with cutoff pi*Wp...

 [b,a] = ellip(n, Rp, Rs, [Wl, Wh])
  band pass filter with band pass edges pi*Wl and pi*Wh ...

 [b,a] = ellip(n, Rp, Rs, [Wl, Wh], 'stop')
  band reject filter with edges pi*Wl and pi*Wh, ...

 [z,p,g] = ellip(...)
  return filter as zero-pole-gain.

 [...] = ellip(...,'s')
     return a Laplace space filter, W can be larger than 1.
 
 [a,b,c,d] = ellip(...)
  return  state-space matrices 

 References: 

 - Oppenheim, Alan V., Discrete Time Signal Processing, Hardcover, 1999.
 - Parente Ribeiro, E., Notas de aula da disciplina TE498 -  Processamento 
   Digital de Sinais, UFPR, 2001/2002.
 - Kienzle, Paul, functions from Octave-Forge, 1999 (http://octave.sf.net).



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 11
 N-ellip 0.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
ellipord


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 346
 usage: [n,wp] = ellipord(wp,ws, rp,rs)

 Calculate the order for the elliptic filter (discrete)
 wp: Cutoff frequency
 ws: Stopband edge
 rp: decibels of ripple in the passband.
 rs: decibels of ripple in the stopband.

 References: 

 - Lamar, Marcus Vinicius, Notas de aula da disciplina TE 456 - Circuitos 
   Analogicos II, UFPR, 2001/2002.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 40
 usage: [n,wp] = ellipord(wp,ws, rp,rs)




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 3
fht


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 788
 -- Function File: m =  fht ( d, n, dim )
     The function fht calculates  Fast Hartley Transform  where D is
     the real input vector (matrix), and M is the real-transform
     vector. For matrices the hartley transform is calculated along the
     columns by default. The options N,and DIM are similar to the
     options of FFT function.

     The forward and inverse hartley transforms are the same (except
     for a scale factor of 1/N for the inverse hartley transform), but
     implemented using different functions .

     The definition of the forward hartley transform for vector d,

     m[K] = \sum_i=0^N-1 d[i]*(cos[K*2*pi*i/N] + sin[K*2*pi*i/N]), for
     0 <= K < N.  m[K] = \sum_i=0^N-1 d[i]*CAS[K*i], for  0 <= K < N.

          fht(1:4)

     See also: ifht, fft





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
The function fht calculates  Fast Hartley Transform  where D is the
real input v



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
filtfilt


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 673
 usage: y = filtfilt(b, a, x)

 Forward and reverse filter the signal. This corrects for phase
 distortion introduced by a one-pass filter, though it does square the
 magnitude response in the process. That's the theory at least.  In
 practice the phase correction is not perfect, and magnitude response
 is distorted, particularly in the stop band.

 Example
    [b, a]=butter(3, 0.1);                   % 10 Hz low-pass filter
    t = 0:0.01:1.0;                         % 1 second sample
    x=sin(2*pi*t*2.3)+0.25*randn(size(t));  % 2.3 Hz sinusoid+noise
    y = filtfilt(b,a,x); z = filter(b,a,x); % apply filter
    plot(t,x,';data;',t,y,';filtfilt;',t,z,';filter;')



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 30
 usage: y = filtfilt(b, a, x)




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
filtic


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 718
 Set initial condition vector for filter function
 The vector zf has the same values that would be obtained 
 from function filter given past inputs x and outputs y

 The vectors x and y contain the most recent inputs and outputs
 respectively, with the newest values first:

 x = [x(-1) x(-2) ... x(-nb)], nb = length(b)-1
 y = [y(-1) y(-2) ... y(-na)], na = length(a)-a

 If length(x)<nb then it is zero padded
 If length(y)<na then it is zero padded

 zf = filtic(b, a, y)
    Initial conditions for filter with coefficients a and b
    and output vector y, assuming input vector x is zero

 zf = filtic(b, a, y, x)
    Initial conditions for filter with coefficients a and b
    input vector x and output vector y



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 Set initial condition vector for filter function
 The vector zf has the same va



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 4
fir1


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1180
 usage: b = fir1(n, w [, type] [, window] [, noscale])

 Produce an order n FIR filter with the given frequency cutoff,
 returning the n+1 filter coefficients in b.  

 n: order of the filter (1 less than the length of the filter)
 w: band edges
    strictly increasing vector in range [0, 1]
    singleton for highpass or lowpass, vector pair for bandpass or
    bandstop, or vector for alternating pass/stop filter.
 type: choose between pass and stop bands
    'high' for highpass filter, cutoff at w
    'stop' for bandstop filter, edges at w = [lo, hi]
    'DC-0' for bandstop as first band of multiband filter
    'DC-1' for bandpass as first band of multiband filter
 window: smoothing window
    defaults to hamming(n+1) row vector
    returned filter is the same shape as the smoothing window
 noscale: choose whether to normalize or not
    'scale': set the magnitude of the center of the first passband to 1
    'noscale': don't normalize

 To apply the filter, use the return vector b:
       y=filter(b,1,x);

 Examples:
   freqz(fir1(40,0.3));
   freqz(fir1(15,[0.2, 0.5], 'stop'));  # note the zero-crossing at 0.1
   freqz(fir1(15,[0.2, 0.5], 'stop', 'noscale'));



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 55
 usage: b = fir1(n, w [, type] [, window] [, noscale])




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 4
fir2


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1197
 usage: b = fir2(n, f, m [, grid_n [, ramp_n]] [, window])

 Produce an FIR filter of order n with arbitrary frequency response,
 returning the n+1 filter coefficients in b.

 n: order of the filter (1 less than the length of the filter)
 f: frequency at band edges
    f is a vector of nondecreasing elements in [0,1]
    the first element must be 0 and the last element must be 1
    if elements are identical, it indicates a jump in freq. response
 m: magnitude at band edges
    m is a vector of length(f)
 grid_n: length of ideal frequency response function
    defaults to 512, should be a power of 2 bigger than n
 ramp_n: transition width for jumps in filter response
    defaults to grid_n/20; a wider ramp gives wider transitions
    but has better stopband characteristics.
 window: smoothing window
    defaults to hamming(n+1) row vector
    returned filter is the same shape as the smoothing window

 To apply the filter, use the return vector b:
       y=filter(b,1,x);
 Note that plot(f,m) shows target response.

 Example:
   f=[0, 0.3, 0.3, 0.6, 0.6, 1]; m=[0, 0, 1, 1/2, 0, 0];
   [h, w] = freqz(fir2(100,f,m));
   plot(f,m,';target response;',w/pi,abs(h),';filter response;');



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 59
 usage: b = fir2(n, f, m [, grid_n [, ramp_n]] [, window])




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 5
firls


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 908
 b = firls(N, F, A);
 b = firls(N, F, A, W);

  FIR filter design using least squares method. Returns a length N+1
  linear phase filter such that the integral of the weighted mean
  squared error in the specified bands is minimized.

  F specifies the frequencies of the band edges, normalized so that
  half the sample frequency is equal to 1.  Each band is specified by
  two frequencies, to the vector must have an even length.

  A specifies the amplitude of the desired response at each band edge.

  W is an optional weighting function that contains one value for each
  band that weights the mean squared error in that band. A must be the
  same length as F, and W must be half the length of F.

 The least squares optimization algorithm for computing FIR filter
 coefficients is derived in detail in:

 I. Selesnick, "Linear-Phase FIR Filter Design by Least Squares,"
 http://cnx.org/content/m10577



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 45
 b = firls(N, F, A);
 b = firls(N, F, A, W);




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 10
flattopwin


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 679
 flattopwin(L, [periodic|symmetric])

 Return the window f(w):

   f(w) = 1 - 1.93 cos(2 pi w) + 1.29 cos(4 pi w)
            - 0.388 cos(6 pi w) + 0.0322cos(8 pi w)

 where w = i/(L-1) for i=0:L-1 for a symmetric window, or
 w = i/L for i=0:L-1 for a periodic window.  The default
 is symmetric.  The returned window is normalized to a peak
 of 1 at w = 0.5.

 This window has low pass-band ripple, but high bandwidth.

 According to [1]:

    The main use for the Flat Top window is for calibration, due
    to its negligible amplitude errors.

 [1] Gade, S; Herlufsen, H; (1987) "Use of weighting functions in DFT/FFT
 analysis (Part I)", Bruel & Kjaer Technical Review No.3.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 37
 flattopwin(L, [periodic|symmetric])




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 9
fracshift


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 279
 -- Function File: [Y H]= fracshift(X,D)
 -- Function File: Y = fracshift(X,D,H)
     Shift the series X by a (possibly fractional) number of samples D.
     The interpolator H is either specified or either designed with a
     Kaiser-windowed sinecard.

   See also: circshift




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 66
Shift the series X by a (possibly fractional) number of samples D.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 5
freqs


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 300
 Usage: H = freqs(B,A,W);

 Compute the s-plane frequency response of the IIR filter B(s)/A(s) as 
 H = polyval(B,j*W)./polyval(A,j*W).  If called with no output
 argument, a plot of magnitude and phase are displayed.

 Example:
    B = [1 2]; A = [1 1];
    w = linspace(0,4,128);
    freqs(B,A,w);



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 26
 Usage: H = freqs(B,A,W);




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 10
freqs_plot


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 90
 -- Function File: freqs_plot (W, H)
     Plot the amplitude and phase of the vector H.





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 45
Plot the amplitude and phase of the vector H.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 4
fwhm


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1318
 Compute peak full-width at half maximum (FWHM) or at another level of peak
 maximum for vector or matrix data y, optionally sampled as y(x). If y is
 a matrix, return FWHM for each column as a row vector.
   Syntax:
      f = fwhm({x, } y {, 'zero'|'min' {, 'rlevel', rlevel}})
      f = fwhm({x, } y {, 'alevel', alevel})
   Examples:
      f = fwhm(y)
      f = fwhm(x, y)
      f = fwhm(x, y, 'zero')
      f = fwhm(x, y, 'min')
      f = fwhm(x, y, 'alevel', 15.3)
      f = fwhm(x, y, 'zero', 'rlevel', 0.5)
      f = fwhm(x, y, 'min',  'rlevel', 0.1)

 The default option 'zero' computes fwhm at half maximum, i.e. 0.5*max(y).
 The option 'min' computes fwhm at the middle curve, i.e. 0.5*(min(y)+max(y)).

 The option 'rlevel' computes full-width at the given relative level of peak
 profile, i.e. at rlevel*max(y) or rlevel*(min(y)+max(y)), respectively.
 For example, fwhm(..., 'rlevel', 0.1) computes full width at 10 % of peak
 maximum with respect to zero or minimum; FWHM is equivalent to
 fwhm(..., 'rlevel', 0.5).

 The option 'alevel' computes full-width at the given absolute level of y.

 Return 0 if FWHM does not exist (e.g. monotonous function or the function
 does not cut horizontal line at rlevel*max(y) or rlevel*(max(y)+min(y)) or
 alevel, respectively).

 Compatibility: Octave 3.x, Matlab



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 Compute peak full-width at half maximum (FWHM) or at another level of peak
 max



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
gauspuls


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 97
 -- Function File: [Y] = gauspuls(T,FC,BW)
     Return the Gaussian modulated sinusoidal pulse.




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 47
Return the Gaussian modulated sinusoidal pulse.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
gaussian


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 442
 usage: w = gaussian(n, a)

 Generate an n-point gaussian convolution window of the given
 width.  Use larger a for a narrower window.  Use larger n for
 longer tails.

     w = exp ( -(a*x)^2/2 )

 for x = linspace ( -(n-1)/2, (n-1)/2, n ).

 Width a is measured in frequency units (sample rate/num samples). 
 It should be f when multiplying in the time domain, but 1/f when 
 multiplying in the frequency domain (for use in convolutions).



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 27
 usage: w = gaussian(n, a)




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
gausswin


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 227
 usage: w = gausswin(L, a)

 Generate an L-point gaussian window of the given width. Use larger a
 for a narrow window.  Use larger L for a smoother curve. 

     w = exp ( -(a*x)^2/2 )

 for x = linspace(-(L-1)/L, (L-1)/L, L)



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 27
 usage: w = gausswin(L, a)




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 9
gmonopuls


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 78
 -- Function File: [Y] = gmonopuls(T,FC)
     Return the gaussian monopulse.




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 30
Return the gaussian monopulse.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
grpdelay


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 2463
 Compute the group delay of a filter.

 [g, w] = grpdelay(b)
   returns the group delay g of the FIR filter with coefficients b.
   The response is evaluated at 512 angular frequencies between 0 and
   pi. w is a vector containing the 512 frequencies.
   The group delay is in units of samples.  It can be converted
   to seconds by multiplying by the sampling period (or dividing by
   the sampling rate fs).

 [g, w] = grpdelay(b,a)
   returns the group delay of the rational IIR filter whose numerator
   has coefficients b and denominator coefficients a.

 [g, w] = grpdelay(b,a,n)
   returns the group delay evaluated at n angular frequencies.  For fastest
   computation n should factor into a small number of small primes.

 [g, w] = grpdelay(b,a,n,'whole')
   evaluates the group delay at n frequencies between 0 and 2*pi.

 [g, f] = grpdelay(b,a,n,Fs)
   evaluates the group delay at n frequencies between 0 and Fs/2.

 [g, f] = grpdelay(b,a,n,'whole',Fs)
   evaluates the group delay at n frequencies between 0 and Fs.

 [g, w] = grpdelay(b,a,w)
   evaluates the group delay at frequencies w (radians per sample).

 [g, f] = grpdelay(b,a,f,Fs)
   evaluates the group delay at frequencies f (in Hz).

 grpdelay(...)
   plots the group delay vs. frequency.

 If the denominator of the computation becomes too small, the group delay
 is set to zero.  (The group delay approaches infinity when
 there are poles or zeros very close to the unit circle in the z plane.)

 Theory: group delay, g(w) = -d/dw [arg{H(e^jw)}],  is the rate of change of
 phase with respect to frequency.  It can be computed as:

               d/dw H(e^-jw)
        g(w) = -------------
                 H(e^-jw)

 where
         H(z) = B(z)/A(z) = sum(b_k z^k)/sum(a_k z^k).

 By the quotient rule,
                    A(z) d/dw B(z) - B(z) d/dw A(z)
        d/dw H(z) = -------------------------------
                               A(z) A(z)
 Substituting into the expression above yields:
                A dB - B dA
        g(w) =  ----------- = dB/B - dA/A
                    A B

 Note that,
        d/dw B(e^-jw) = sum(k b_k e^-jwk)
        d/dw A(e^-jw) = sum(k a_k e^-jwk)
 which is just the FFT of the coefficients multiplied by a ramp.

 As a further optimization when nfft>>length(a), the IIR filter (b,a)
 is converted to the FIR filter conv(b,fliplr(conj(a))).
 For further details, see
 http://ccrma.stanford.edu/~jos/filters/Numerical_Computation_Group_Delay.html



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 37
 Compute the group delay of a filter.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 4
hann


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 28
 w = hann(n)
   see hanning



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 28
 w = hann(n)
   see hanning




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 7
hilbert


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 647
 -- Function File: H = hilbert (F,N,DIM)
     Analytic extension of real valued signal

     `H=hilbert(F)' computes the extension of the real valued signal F
     to an analytic signal. If F is a matrix, the transformation is
     applied to each column. For N-D arrays, the transformation is
     applied to the first non-singleton dimension.

     `real(H)' contains the original signal F.  `imag(H)' contains the
     Hilbert transform of F.

     `hilbert(F,N)' does the same using a length N Hilbert transform.
     The result will also have length N.

     `hilbert(F,[],DIM)' or `hilbert(F,N,DIM)' does the same along
     dimension dim.




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 41
Analytic extension of real valued signal




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 4
idct


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 584
 y = dct (x, n)
    Computes the inverse discrete cosine transform of x.  If n is
    given, then x is padded or trimmed to length n before computing
    the transform. If x is a matrix, compute the transform along the
    columns of the the matrix. The transform is faster if x is
    real-valued and even length.

 The inverse discrete cosine transform x of X can be defined as follows:

          N-1
   x[n] = sum w(k) X[k] cos (pi (2n+1) k / 2N ),  n = 0, ..., N-1
          k=0

 with w(0) = sqrt(1/N) and w(k) = sqrt(2/N), k = 1, ..., N-1

 See also: idct, dct2, idct2, dctmtx



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 72
 y = dct (x, n)
    Computes the inverse discrete cosine transform of x.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 5
idct2


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 221
 y = idct2 (x)
   Computes the inverse 2-D discrete cosine transform of matrix x

 y = idct2 (x, m, n) or y = idct2 (x, [m n])
   Computes the 2-D inverse DCT of x after padding or trimming rows to m and
   columns to n.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 y = idct2 (x)
   Computes the inverse 2-D discrete cosine transform of matrix x



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 4
idst


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 333
 -- Function File: Y = idst (X)
 -- Function File: Y = idst (X, N)
     Computes the inverse type I discrete sine transform of Y.  If N is
     given, then Y is padded or trimmed to length N before computing
     the transform.  If Y is a matrix, compute the transform along the
     columns of the the matrix.

     See also: dst





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 57
Computes the inverse type I discrete sine transform of Y.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 4
ifht


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 805
 -- Function File: m =  ifht ( d, n, dim )
     The function ifht calculates  Fast Hartley Transform  where D is
     the real input vector (matrix), and M is the real-transform
     vector. For matrices the hartley transform is calculated along the
     columns by default. The options N, and DIM are similar to the
     options of FFT function.

     The forward and inverse hartley transforms are the same (except
     for a scale factor of 1/N for the inverse hartley transform), but
     implemented using different functions .

     The definition of the forward hartley transform for vector d,

     m[K] = 1/N \sum_i=0^N-1 d[i]*(cos[K*2*pi*i/N] + sin[K*2*pi*i/N]),
     for  0 <= K < N.  m[K] = 1/N \sum_i=0^N-1 d[i]*CAS[K*i], for  0 <=
     K < N.

          ifht(1:4)

     See also: fht, fft





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
The function ifht calculates  Fast Hartley Transform  where D is the
real input 



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
iirlp2mb


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1305
 IIR Low Pass Filter to Multiband Filter Transformation

 [Num,Den,AllpassNum,AllpassDen] = iirlp2mb(B,A,Wo,Wt)
 [Num,Den,AllpassNum,AllpassDen] = iirlp2mb(B,A,Wo,Wt,Pass)

 Num,Den:               numerator,denominator of the transformed filter
 AllpassNum,AllpassDen: numerator,denominator of allpass transform,
 B,A:                   numerator,denominator of prototype low pass filter
 Wo:                    normalized_angular_frequency/pi to be transformed
 Wt:                    [phi=normalized_angular_frequencies]/pi target vector
 Pass:                  This parameter may have values 'pass' or 'stop'.  If
                        not given, it defaults to the value of 'pass'.

 With normalized ang. freq. targets 0 < phi(1) <  ... < phi(n) < pi radians

 for Pass == 'pass', the target multiband magnitude will be:
       --------       ----------        -----------...
      /        \     /          \      /            .
 0   phi(1) phi(2)  phi(3)   phi(4)   phi(5)   (phi(6))    pi

 for Pass == 'stop', the target multiband magnitude will be:
 -------      ---------        ----------...
        \    /         \      /           .
 0   phi(1) phi(2)  phi(3)   phi(4)  (phi(5))              pi

 Example of use:
 [B, A] = butter(6, 0.5);
 [Num, Den] = iirlp2mb(B, A, 0.5, [.2 .4 .6 .8]);



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 56
 IIR Low Pass Filter to Multiband Filter Transformation




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
impinvar


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 877
 -- Function File: [B_OUT, A_OUT] = impinvar (B, A, FS, TOL)
 -- Function File: [B_OUT, A_OUT] = impinvar (B, A, FS)
 -- Function File: [B_OUT, A_OUT] = impinvar (B, A)
     Converts analog filter with coefficients B and A to digital,
     conserving impulse response.

     If FS is not specificied, or is an empty vector, it defaults to
     1Hz.

     If TOL is not specified, it defaults to 0.0001 (0.1%) This
     function does the inverse of impinvar so that the following
     example should restore the original values of A and B.

     `invimpinvar' implements the reverse of this function.
          [b, a] = impinvar (b, a);
          [b, a] = invimpinvar (b, a);

     Reference: Thomas J. Cavicchi (1996) "Impulse invariance and
     multiple-order poles". IEEE transactions on signal processing, Vol
     40 (9): 2344-2347

     See also: bilinear, invimpinvar





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
Converts analog filter with coefficients B and A to digital, conserving
impulse 



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 4
impz


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 545
 usage: [x, t] = impz(b [, a, n, fs])

 Generate impulse-response characteristics of the filter. The filter
 coefficients correspond to the the z-plane rational function with
 numerator b and denominator a.  If a is not specified, it defaults to
 1. If n is not specified, or specified as [], it will be chosen such
 that the signal has a chance to die down to -120dB, or to not explode
 beyond 120dB, or to show five periods if there is no significant
 damping. If no return arguments are requested, plot the results.

 See also: freqz, zplane



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 38
 usage: [x, t] = impz(b [, a, n, fs])




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
interp


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 631
 usage: y = interp(x, q [, n [, Wc]])

 Upsample the signal x by a factor of q, using an order 2*q*n+1 FIR
 filter. Note that q must be an integer for this rate change method.
 n defaults to 4 and Wc defaults to 0.5.

 Example
                                          # Generate a signal.
    t=0:0.01:2; x=chirp(t,2,.5,10,'quadratic')+sin(2*pi*t*0.4); 
    y = interp(x(1:4:length(x)),4,4,1);   # interpolate a sub-sample
    stem(t(1:121)*1000,x(1:121),"-g;Original;"); hold on;
    stem(t(1:121)*1000,y(1:121),"-r;Interpolated;");
    stem(t(1:4:121)*1000,x(1:4:121),"-b;Subsampled;"); hold off;

 See also: decimate, resample



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 38
 usage: y = interp(x, q [, n [, Wc]])




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 7
invfreq


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1687
 usage: [B,A] = invfreq(H,F,nB,nA)
        [B,A] = invfreq(H,F,nB,nA,W)
        [B,A] = invfreq(H,F,nB,nA,W,[],[],plane)
        [B,A] = invfreq(H,F,nB,nA,W,iter,tol,plane)

 Fit filter B(z)/A(z) or B(s)/A(s) to complex frequency response at 
 frequency points F. A and B are real polynomial coefficients of order 
 nA and nB respectively.  Optionally, the fit-errors can be weighted vs 
 frequency according to the weights W. Also, the transform plane can be
 specified as either 's' for continuous time or 'z' for discrete time. 'z'
 is chosen by default.  Eventually, Steiglitz-McBride iterations will be
 specified by iter and tol.

 H: desired complex frequency response
     It is assumed that A and B are real polynomials, hence H is one-sided.
 F: vector of frequency samples in radians
 nA: order of denominator polynomial A
 nB: order of numerator polynomial B
 plane='z': F on unit circle (discrete-time spectra, z-plane design)
 plane='s': F on jw axis     (continuous-time spectra, s-plane design)
 H(k) = spectral samples of filter frequency response at points zk,
  where zk=exp(sqrt(-1)*F(k)) when plane='z' (F(k) in [0,.5])
     and zk=(sqrt(-1)*F(k)) when plane='s' (F(k) nonnegative)
 Example:
     [B,A] = butter(12,1/4);
     [H,w] = freqz(B,A,128);
     [Bh,Ah] = invfreq(H,F,4,4);
     Hh = freqz(Bh,Ah);
     disp(sprintf('||frequency response error|| = %f',norm(H-Hh)));

 References: J. O. Smith, "Techniques for Digital Filter Design and System 
  	Identification with Application to the Violin, Ph.D. Dissertation, 
 	Elec. Eng. Dept., Stanford University, June 1983, page 50; or,

 http://ccrma.stanford.edu/~jos/filters/FFT_Based_Equation_Error_Method.html



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 usage: [B,A] = invfreq(H,F,nB,nA)
        [B,A] = invfreq(H,F,nB,nA,W)
        



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
invfreqs


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 943
 Usage: [B,A] = invfreqs(H,F,nB,nA)
        [B,A] = invfreqs(H,F,nB,nA,W)
        [B,A] = invfreqs(H,F,nB,nA,W,iter,tol,'trace')

 Fit filter B(s)/A(s)to the complex frequency response H at frequency
 points F.  A and B are real polynomial coefficients of order nA and nB.
 Optionally, the fit-errors can be weighted vs frequency according to
 the weights W.
 Note: all the guts are in invfreq.m

 H: desired complex frequency response
 F: frequency (must be same length as H)
 nA: order of the denominator polynomial A
 nB: order of the numerator polynomial B
 W: vector of weights (must be same length as F)

 Example:
       B = [1/2 1];
       A = [1 1];
       w = linspace(0,4,128);
       H = freqs(B,A,w);
       [Bh,Ah] = invfreqs(H,w,1,1);
       Hh = freqs(Bh,Ah,w);
       plot(w,[abs(H);abs(Hh)])
       legend('Original','Measured');
       err = norm(H-Hh);
       disp(sprintf('L2 norm of frequency response error = %f',err));



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 Usage: [B,A] = invfreqs(H,F,nB,nA)
        [B,A] = invfreqs(H,F,nB,nA,W)
      



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
invfreqz


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 819
 usage: [B,A] = invfreqz(H,F,nB,nA)
        [B,A] = invfreqz(H,F,nB,nA,W)
        [B,A] = invfreqz(H,F,nB,nA,W,iter,tol,'trace')

 Fit filter B(z)/A(z)to the complex frequency response H at frequency
 points F.  A and B are real polynomial coefficients of order nA and nB.
 Optionally, the fit-errors can be weighted vs frequency according to
 the weights W.
 Note: all the guts are in invfreq.m

 H: desired complex frequency response
 F: normalized frequncy (0 to pi) (must be same length as H)
 nA: order of the denominator polynomial A
 nB: order of the numerator polynomial B
 W: vector of weights (must be same length as F)

 Example:
     [B,A] = butter(4,1/4);
     [H,F] = freqz(B,A);
     [Bh,Ah] = invfreq(H,F,4,4);
     Hh = freqz(Bh,Ah);
     disp(sprintf('||frequency response error|| = %f',norm(H-Hh)));



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 usage: [B,A] = invfreqz(H,F,nB,nA)
        [B,A] = invfreqz(H,F,nB,nA,W)
      



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 11
invimpinvar


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 823
 -- Function File: [B_OUT, A_OUT] = invimpinvar (B, A, FS, TOL)
 -- Function File: [B_OUT, A_OUT] = invimpinvar (B, A, FS)
 -- Function File: [B_OUT, A_OUT] = invimpinvar (B, A)
     Converts digital filter with coefficients B and A to analog,
     conserving impulse response.

     This function does the inverse of impinvar so that the following
     example should restore the original values of A and B.
          [b, a] = impinvar (b, a);
          [b, a] = invimpinvar (b, a);

     If FS is not specificied, or is an empty vector, it defaults to
     1Hz.

     If TOL is not specified, it defaults to 0.0001 (0.1%)

     Reference: Thomas J. Cavicchi (1996) "Impulse invariance and
     multiple-order poles". IEEE transactions on signal processing, Vol
     40 (9): 2344-2347

     See also: bilinear, impinvar





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
Converts digital filter with coefficients B and A to analog, conserving
impulse 



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
kaiser


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 454
 usage:  kaiser (L, beta)

 Returns the filter coefficients of the L-point Kaiser window with
 parameter beta.

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

 The continuous version of width L centered about x=0 is:

         besseli(0, beta * sqrt(1-(2*x/L).^2))
 k(x) =  -------------------------------------,  L/2 <= x <= L/2
                besseli(0, beta)

 See also: kaiserord



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 26
 usage:  kaiser (L, beta)




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 9
kaiserord


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1809
 usage: [n, Wn, beta, ftype] = kaiserord(f, m, dev [, fs])

 Returns the parameters needed for fir1 to produce a filter of the
 desired specification from a kaiser window:
       n: order of the filter (length of filter minus 1)
       Wn: band edges for use in fir1
       beta: parameter for kaiser window of length n+1
       ftype: choose between pass and stop bands
       b = fir1(n,Wn,kaiser(n+1,beta),ftype,'noscale');

 f: frequency bands, given as pairs, with the first half of the
    first pair assumed to start at 0 and the last half of the last
    pair assumed to end at 1.  It is important to separate the
    band edges, since narrow transition regions require large order
    filters.
 m: magnitude within each band.  Should be non-zero for pass band
    and zero for stop band.  All passbands must have the same
    magnitude, or you will get the error that pass and stop bands
    must be strictly alternating.
 dev: deviation within each band.  Since all bands in the resulting
    filter have the same deviation, only the minimum deviation is
    used.  In this version, a single scalar will work just as well.
 fs: sampling rate.  Used to convert the frequency specification into
    the [0, 1], where 1 corresponds to the Nyquist frequency, fs/2.

 The Kaiser window parameters n and beta are computed from the
 relation between ripple (A=-20*log10(dev)) and transition width
 (dw in radians) discovered empirically by Kaiser:

           / 0.1102(A-8.7)                        A > 50
    beta = | 0.5842(A-21)^0.4 + 0.07886(A-21)     21 <= A <= 50
           \ 0.0                                  A < 21

    n = (A-8)/(2.285 dw)

 Example
    [n, w, beta, ftype] = kaiserord([1000,1200], [1,0], [0.05,0.05], 11025);
    freqz(fir1(n,w,kaiser(n+1,beta),ftype,'noscale'),1,[],11025);



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 59
 usage: [n, Wn, beta, ftype] = kaiserord(f, m, dev [, fs])




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
levinson


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 818
 usage:  [a, v, ref] = levinson (acf [, p])

 Use the Durbin-Levinson algorithm to solve:
    toeplitz(acf(1:p)) * x = -acf(2:p+1).
 The solution [1, x'] is the denominator of an all pole filter
 approximation to the signal x which generated the autocorrelation
 function acf.  

 acf is the autocorrelation function for lags 0 to p.
 p defaults to length(acf)-1.
 Returns 
   a=[1, x'] the denominator filter coefficients. 
   v= variance of the white noise = square of the numerator constant
   ref = reflection coefficients = coefficients of the lattice
         implementation of the filter
 Use freqz(sqrt(v),a) to plot the power spectrum.

 REFERENCE
 [1] Steven M. Kay and Stanley Lawrence Marple Jr.:
   "Spectrum analysis -- a modern perspective",
   Proceedings of the IEEE, Vol 69, pp 1380-1419, Nov., 1981



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 44
 usage:  [a, v, ref] = levinson (acf [, p])




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 7
marcumq


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 963
 -- Function File: Q =  marcumq (A, B)
 -- Function File: Q =  marcumq (A, B, M)
 -- Function File: Q =  marcumq (A, B, M, TOL)
     Compute the generalized Marcum Q function of order M with
     noncentrality parameter A and argument B.  If the order M is
     omitted it defaults to 1.  An optional relative tolerance TOL may
     be included, the default is `eps'.

     If the input arguments are commensurate vectors, this function
     will produce a table of values.

     This function computes Marcum's Q function using the infinite
     Bessel series, truncated when the relative error is less than the
     specified tolerance.  The accuracy is limited by that of the
     Bessel functions, so reducing the tolerance is probably not useful.

     Reference: Marcum, "Tables of Q Functions", Rand Corporation.

     Reference: R.T. Short, "Computation of Noncentral Chi-squared and
     Rice Random Variables", www.phaselockedsystems.com/publications





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
Compute the generalized Marcum Q function of order M with noncentrality
paramete



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 7
mexihat


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 85
 -- Function File: [PSI,X] = mexihat(LB,UB,N)
     Compute the Mexican hat wavelet.




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 32
Compute the Mexican hat wavelet.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
meyeraux


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 89
 -- Function File: [Y] = meyeraux(X)
     Compute the Meyer wavelet auxiliary function.




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 45
Compute the Meyer wavelet auxiliary function.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
morlet


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 79
 -- Function File: [PSI,X] = morlet(LB,UB,N)
     Compute the Morlet wavelet.




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 27
Compute the Morlet wavelet.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
mscohere


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 270
 Usage:
   [Pxx,freq]=mscohere(x,y,Nfft,Fs,window,overlap,range,plot_type,detrend)

     Estimate (mean square) coherence of signals "x" and "y".
     Use the Welch (1967) periodogram/FFT method.
     See "help pwelch" for description of arguments, hints and references



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 Usage:
   [Pxx,freq]=mscohere(x,y,Nfft,Fs,window,overlap,range,plot_type,detren



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
ncauer


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 383
 usage: [Zz, Zp, Zg] = ncauer(Rp, Rs, n)

 Analog prototype for Cauer filter.
 [z, p, g]=ncauer(Rp, Rs, ws)
 Rp = Passband ripple
 Rs = Stopband ripple
 Ws = Desired order

 References: 

 - Serra, Celso Penteado, Teoria e Projeto de Filtros, Campinas: CARTGRAF, 
   1983.
 - Lamar, Marcus Vinicius, Notas de aula da disciplina TE 456 - Circuitos 
   Analogicos II, UFPR, 2001/2002.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 41
 usage: [Zz, Zp, Zg] = ncauer(Rp, Rs, n)




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 10
nuttallwin


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 154
 -- Function File: [W] = nuttallwin(L)
     Compute the Blackman-Harris window defined by Nuttall of length L.

     See also: blackman, blackmanharris





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 66
Compute the Blackman-Harris window defined by Nuttall of length L.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 9
parzenwin


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 118
 -- Function File: [W] = parzenwin(L)
     Compute the Parzen window of lenght L.

     See also: rectwin, bartlett





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 38
Compute the Parzen window of lenght L.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 5
pburg


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 3783
 usage:
    [psd,f_out] = pburg(x,poles,freq,Fs,range,method,plot_type,criterion)

 Calculate Burg maximum-entropy power spectral density.
 The functions "arburg" and "ar_psd" do all the work.
 See "help arburg" and "help ar_psd" for further details.

 ARGUMENTS:
     All but the first two arguments are optional and may be empty.
   x       %% [vector] sampled data

   poles   %% [integer scalar] required number of poles of the AR model

   freq    %% [real vector] frequencies at which power spectral density
           %%               is calculated
           %% [integer scalar] number of uniformly distributed frequency
           %%          values at which spectral density is calculated.
           %%          [default=256]

   Fs      %% [real scalar] sampling frequency (Hertz) [default=1]


 CONTROL-STRING ARGUMENTS -- each of these arguments is a character string.
   Control-string arguments can be in any order after the other arguments.


   range   %% 'half',  'onesided' : frequency range of the spectrum is
           %%       from zero up to but not including sample_f/2.  Power
           %%       from negative frequencies is added to the positive
           %%       side of the spectrum.
           %% 'whole', 'twosided' : frequency range of the spectrum is
           %%       -sample_f/2 to sample_f/2, with negative frequencies
           %%       stored in "wrap around" order after the positive
           %%       frequencies; e.g. frequencies for a 10-point 'twosided'
           %%       spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1
           %% 'shift', 'centerdc' : same as 'whole' but with the first half
           %%       of the spectrum swapped with second half to put the
           %%       zero-frequency value in the middle. (See "help
           %%       fftshift". If "freq" is vector, 'shift' is ignored.
           %% If model coefficients "ar_coeffs" are real, the default
           %% range is 'half', otherwise default range is 'whole'.

   method  %% 'fft':  use FFT to calculate power spectral density.
           %% 'poly': calculate spectral density as a polynomial of 1/z
           %% N.B. this argument is ignored if the "freq" argument is a
           %%      vector.  The default is 'poly' unless the "freq"
           %%      argument is an integer power of 2.
   
 plot_type %% 'plot', 'semilogx', 'semilogy', 'loglog', 'squared' or 'db':
           %% specifies the type of plot.  The default is 'plot', which
           %% means linear-linear axes. 'squared' is the same as 'plot'.
           %% 'dB' plots "10*log10(psd)".  This argument is ignored and a
           %% spectrum is not plotted if the caller requires a returned
           %% value.

 criterion %% [optional string arg]  model-selection criterion.  Limits
           %%       the number of poles so that spurious poles are not 
           %%       added when the whitened data has no more information
           %%       in it (see Kay & Marple, 1981). Recognised values are
           %%  'AKICc' -- approximate corrected Kullback information
           %%             criterion (recommended),
           %%   'KIC'  -- Kullback information criterion
           %%   'AICc' -- corrected Akaike information criterion
           %%   'AIC'  -- Akaike information criterion
           %%   'FPE'  -- final prediction error" criterion
           %% The default is to NOT use a model-selection criterion

 RETURNED VALUES:
     If return values are not required by the caller, the spectrum
     is plotted and nothing is returned.
   psd       %% [real vector] power-spectral density estimate 
   f_out     %% [real vector] frequency values 

 HINTS
   This function is a wrapper for arburg and ar_psd.
   See "help arburg", "help ar_psd".



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 usage:
    [psd,f_out] = pburg(x,poles,freq,Fs,range,method,plot_type,criterion



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 15
pei_tseng_notch


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 707
 -- Function File:  [ B, A ]  = pei_tseng_notch ( FREQUENCIES,
          BANDWIDTHS
     Return coefficients for an IIR notch-filter with one or more
     filter frequencies and according (very narrow) bandwidths to be
     used with `filter' or `filtfilt'.  The filter construction is
     based on an allpass which performs a reversal of phase at the
     filter frequencies.  Thus, the mean of the phase-distorted and the
     original signal has the respective frequencies removed.  See the
     demo for an illustration.

     Original source: Pei, Soo-Chang, and Chien-Cheng Tseng "IIR
     Multiple Notch Filter Design Based on Allpass Filter" 1996 IEEE
     Tencon doi: 10.1109/TENCON.1996.608814)




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
Return coefficients for an IIR notch-filter with one or more filter
frequencies 



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
polystab


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 156
 b = polystab(a)

 Stabalize the polynomial transfer function by replacing all roots
 outside the unit circle with their reflection inside the unit circle.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 17
 b = polystab(a)




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
pulstran


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1187
 usage: y=pulstran(t,d,'func',...)
        y=pulstran(t,d,p,Fs,'interp')

 Generate the signal y=sum(func(t+d,...)) for each d.  If d is a
 matrix of two columns, the first column is the delay d and the second
 column is the amplitude a, and y=sum(a*func(t+d)) for each d,a.
 Clearly, func must be a function which accepts a vector of times.
 Any extra arguments needed for the function must be tagged on the end.

 Example
   fs = 11025;  # arbitrary sample rate
   f0 = 100;    # pulse train sample rate
   w = 0.001;   # pulse width of 1 millisecond
   auplot(pulstran(0:1/fs:0.1, 0:1/f0:0.1, 'rectpuls', w), fs);

 If instead of a function name you supply a pulse shape sampled at
 frequency Fs (default 1 Hz),  an interpolated version of the pulse
 is added at each delay d.  The interpolation stays within the the
 time range of the delayed pulse.  The interpolation method defaults
 to linear, but it can be any interpolation method accepted by the
 function interp1.

 Example
   fs = 11025;  # arbitrary sample rate
   f0 = 100;    # pulse train sample rate
   w = boxcar(10);  # pulse width of 1 millisecond at 10 kHz
   auplot(pulstran(0:1/fs:0.1, 0:1/f0:0.1, w, 10000), fs);



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 31
 usage: y=pulstran(t,d,'func',.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
pwelch


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 7140
 USAGE:
   [spectra,freq] = pwelch(x,window,overlap,Nfft,Fs,
                           range,plot_type,detrend,sloppy)
     Estimate power spectral density of data "x" by the Welch (1967)
     periodogram/FFT method.  All arguments except "x" are optional.
         The data is divided into segments.  If "window" is a vector, each
     segment has the same length as "window" and is multiplied by "window"
     before (optional) zero-padding and calculation of its periodogram. If
     "window" is a scalar, each segment has a length of "window" and a
     Hamming window is used.
         The spectral density is the mean of the periodograms, scaled so that
     area under the spectrum is the same as the mean square of the
     data.  This equivalence is supposed to be exact, but in practice there
     is a mismatch of up to 0.5% when comparing area under a periodogram
     with the mean square of the data.

  [spectra,freq] = pwelch(x,y,window,overlap,Nfft,Fs,
                          range,plot_type,detrend,sloppy,results)
     Two-channel spectrum analyser.  Estimate power spectral density, cross-
     spectral density, transfer function and/or coherence functions of time-
     series input data "x" and output data "y" by the Welch (1967)
     periodogram/FFT method.
       pwelch treats the second argument as "y" if there is a control-string
     argument "cross", "trans", "coher" or "ypower"; "power" does not force
     the 2nd argument to be treated as "y".  All other arguments are
     optional.  All spectra are returned in matrix "spectra". 

  [spectra,Pxx_ci,freq] = pwelch(x,window,overlap,Nfft,Fs,conf,
                                 range,plot_type,detrend,sloppy)
  [spectra,Pxx_ci,freq] = pwelch(x,y,window,overlap,Nfft,Fs,conf,
                                 range,plot_type,detrend,sloppy,results)
     Estimates confidence intervals for the spectral density.
     See Hint (7) below for compatibility options.  Confidence level "conf"
     is the 6th or 7th numeric argument.  If "results" control-string 
     arguments are used, one of them must be "power" when the "conf"
     argument is present; pwelch can estimate confidence intervals only for
     the power spectrum of the "x" data.  It does not know how to estimate
     confidence intervals of the cross-power spectrum, transfer function or
     coherence; if you can suggest a good method, please send a bug report.

 ARGUMENTS
 All but the first argument are optional and may be empty, except that
 the "results" argument may require the second argument to be "y".

 x           %% [non-empty vector] system-input time-series data
 y           %% [non-empty vector] system-output time-series data

 window      %% [real vector] of window-function values between 0 and 1; the
             %%       data segment has the same length as the window.
             %%       Default window shape is Hamming.
             %% [integer scalar] length of each data segment.  The default
             %%       value is window=sqrt(length(x)) rounded up to the
             %%       nearest integer power of 2; see 'sloppy' argument.

 overlap     %% [real scalar] segment overlap expressed as a multiple of
             %%       window or segment length.   0 <= overlap < 1,
             %%       The default is overlap=0.5 .

 Nfft        %% [integer scalar] Length of FFT.  The default is the length
             %%       of the "window" vector or has the same value as the
             %%       scalar "window" argument.  If Nfft is larger than the
             %%       segment length, "seg_len", the data segment is padded
             %%       with "Nfft-seg_len" zeros.  The default is no padding.
             %%       Nfft values smaller than the length of the data
             %%       segment (or window) are ignored silently.

 Fs          %% [real scalar] sampling frequency (Hertz); default=1.0

 conf        %% [real scalar] confidence level between 0 and 1.  Confidence
             %%       intervals of the spectral density are estimated from
             %%       scatter in the periodograms and are returned as Pxx_ci.
             %%       Pxx_ci(:,1) is the lower bound of the confidence
             %%       interval and Pxx_ci(:,2) is the upper bound.  If there
             %%       are three return values, or conf is an empty matrix,
             %%       confidence intervals are calculated for conf=0.95 .
             %%       If conf is zero or is not given, confidence intervals
             %%       are not calculated. Confidence intervals can be 
             %%       obtained only for the power spectral density of x;
             %%       nothing else.

 CONTROL-STRING ARGUMENTS -- each of these arguments is a character string.
   Control-string arguments must be after the other arguments but can be in
   any order.
  
 range     %% 'half',  'onesided' : frequency range of the spectrum is
           %%       zero up to but not including Fs/2.  Power from
           %%       negative frequencies is added to the positive side of
           %%       the spectrum, but not at zero or Nyquist (Fs/2)
           %%       frequencies.  This keeps power equal in time and
           %%       spectral domains.  See reference [2].
           %% 'whole', 'twosided' : frequency range of the spectrum is
           %%       -Fs/2 to Fs/2, with negative frequencies
           %%       stored in "wrap around" order after the positive
           %%       frequencies; e.g. frequencies for a 10-point 'twosided'
           %%       spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1
           %% 'shift', 'centerdc' : same as 'whole' but with the first half
           %%       of the spectrum swapped with second half to put the
           %%       zero-frequency value in the middle. (See "help
           %%       fftshift".
           %% If data (x and y) are real, the default range is 'half',
           %% otherwise default range is 'whole'.

 plot_type %% 'plot', 'semilogx', 'semilogy', 'loglog', 'squared' or 'db':
           %% specifies the type of plot.  The default is 'plot', which
           %% means linear-linear axes. 'squared' is the same as 'plot'.
           %% 'dB' plots "10*log10(psd)".  This argument is ignored and a
           %% spectrum is not plotted if the caller requires a returned
           %% value.

 detrend   %% 'no-strip', 'none' -- do NOT remove mean value from the data
           %% 'short', 'mean' -- remove the mean value of each segment from
           %%                    each segment of the data.
           %% 'linear',       -- remove linear trend from each segment of
           %%                    the data.
           %% 'long-mean'     -- remove the mean value from the data before
           %%              splitting it into segments.  This is the default.

   sloppy  %% 'sloppy': FFT length is rounded up to the nearest integer
           %%       power of 2 by zero padding.  FFT length is adjusted
           %%       after addition of padding by explicit Nfft argument.
           %%       The default is to use exactly the FFT and window/



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 USAGE:
   [spectra,freq] = pwelch(x,window,overlap,Nfft,Fs,
                   



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 7
pyulear


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 3140
 usage:
    [psd,f_out] = pyulear(x,poles,freq,Fs,range,method,plot_type)

 Calculates a Yule-Walker autoregressive (all-pole) model of the data "x"
 and computes the power spectrum of the model.  This is a wrapper for
 functions "aryule" and "ar_psd" which perform the argument checking.
 See "help aryule" and "help ar_psd" for further details.

 ARGUMENTS:
     All but the first two arguments are optional and may be empty.
   x       %% [vector] sampled data

   poles   %% [integer scalar] required number of poles of the AR model

   freq    %% [real vector] frequencies at which power spectral density
           %%               is calculated
           %% [integer scalar] number of uniformly distributed frequency
           %%          values at which spectral density is calculated.
           %%          [default=256]

   Fs      %% [real scalar] sampling frequency (Hertz) [default=1]


 CONTROL-STRING ARGUMENTS -- each of these arguments is a character string.
   Control-string arguments can be in any order after the other arguments.


   range   %% 'half',  'onesided' : frequency range of the spectrum is
           %%       from zero up to but not including sample_f/2.  Power
           %%       from negative frequencies is added to the positive
           %%       side of the spectrum.
           %% 'whole', 'twosided' : frequency range of the spectrum is
           %%       -sample_f/2 to sample_f/2, with negative frequencies
           %%       stored in "wrap around" order after the positive
           %%       frequencies; e.g. frequencies for a 10-point 'twosided'
           %%       spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1
           %% 'shift', 'centerdc' : same as 'whole' but with the first half
           %%       of the spectrum swapped with second half to put the
           %%       zero-frequency value in the middle. (See "help
           %%       fftshift". If "freq" is vector, 'shift' is ignored.
           %% If model coefficients "ar_coeffs" are real, the default
           %% range is 'half', otherwise default range is 'whole'.

   method  %% 'fft':  use FFT to calculate power spectrum.
           %% 'poly': calculate power spectrum as a polynomial of 1/z
           %% N.B. this argument is ignored if the "freq" argument is a
           %%      vector.  The default is 'poly' unless the "freq"
           %%      argument is an integer power of 2.
   
 plot_type %% 'plot', 'semilogx', 'semilogy', 'loglog', 'squared' or 'db':
           %% specifies the type of plot.  The default is 'plot', which
           %% means linear-linear axes. 'squared' is the same as 'plot'.
           %% 'dB' plots "10*log10(psd)".  This argument is ignored and a
           %% spectrum is not plotted if the caller requires a returned
           %% value.

 RETURNED VALUES:
     If return values are not required by the caller, the spectrum
     is plotted and nothing is returned.
   psd     %% [real vector] power-spectrum estimate 
   f_out   %% [real vector] frequency values 

 HINTS
   This function is a wrapper for aryule and ar_psd.
   See "help aryule", "help ar_psd".



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 74
 usage:
    [psd,f_out] = pyulear(x,poles,freq,Fs,range,method,plot_type)




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 9
qp_kaiser


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 602
 Usage:  qp_kaiser (nb, at, linear)

 Computes a finite impulse response (FIR) filter for use with a
 quasi-perfect reconstruction polyphase-network filter bank. This
 version utilizes a Kaiser window to shape the frequency response of
 the designed filter. Tha number nb of bands and the desired
 attenuation at in the stop-band are given as parameters.

 The Kaiser window is multiplied by the ideal impulse response
 h(n)=a.sinc(a.n) and converted to its minimum-phase version by means
 of a Hilbert transform.

 By using a third non-null argument, the minimum-phase calculation is
 ommited at all.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 36
 Usage:  qp_kaiser (nb, at, linear)




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 5
rceps


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 693
 usage: [y, xm] = rceps(x)
   Produce the cepstrum of the signal x, and if desired, the minimum
   phase reconstruction of the signal x.  If x is a matrix, do so
   for each column of the matrix.

 Example
   f0=70; Fs=10000;           # 100 Hz fundamental, 10kHz sampling rate
   a=poly(0.985*exp(1i*pi*[0.1, -0.1, 0.3, -0.3])); # two formants
   s=0.005*randn(1024,1);      # Noise excitation signal
   s(1:Fs/f0:length(s)) = 1;   # Impulse glottal wave
   x=filter(1,a,s);            # Speech signal in x
   [y, xm] = rceps(x.*hanning(1024)); # cepstrum and min phase reconstruction

 Reference
    Programs for digital signal processing. IEEE Press.
    New York: John Wiley & Sons. 1979.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 usage: [y, xm] = rceps(x)
   Produce the cepstrum of the signal x, and if desir



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
rectpuls


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 429
 usage: y = rectpuls(t, w)

 Generate a rectangular pulse over the interval [-w/2,w/2), sampled at
 times t.  This is useful with the function pulstran for generating a
 series pulses.

 Example
   fs = 11025;  # arbitrary sample rate
   f0 = 100;    # pulse train sample rate
   w = 0.3/f0;  # pulse width 3/10th the distance between pulses
   auplot(pulstran(0:1/fs:4/f0, 0:1/f0:4/f0, 'rectpuls', w), fs);

 See also: pulstran



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 27
 usage: y = rectpuls(t, w)




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 7
rectwin


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 142
 -- Function File: [W] = rectwin(L)
     Return the filter coefficients of a rectangle window of length L.

     See also: hamming, hanning





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 65
Return the filter coefficients of a rectangle window of length L.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
resample


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 637
 -- Function File: [Y H]= resample(X,P,Q)
 -- Function File: Y = resample(X,P,Q,H)
     Change the sample rate of X by a factor of P/Q. This is performed
     using a polyphase algorithm. The impulse response H of the
     antialiasing filter is either specified or either designed with a
     Kaiser-windowed sinecard.

     Ref [1] J. G. Proakis and D. G. Manolakis, Digital Signal
     Processing: Principles, Algorithms, and Applications, 4th ed.,
     Prentice Hall, 2007. Chap. 6

     Ref [2] A. V. Oppenheim, R. W. Schafer and J. R. Buck,
     Discrete-time signal processing, Signal processing series,
     Prentice-Hall, 1999




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 47
Change the sample rate of X by a factor of P/Q.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
residued


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1061
 -- Function File: [R, P, F, M] = residued (B, A)
     Compute the partial fraction expansion (PFE) of filter H(z) =
     B(z)/A(z).  In the usual PFE function `residuez', the IIR part
     (poles P and residues R) is driven _in parallel_ with the FIR part
     (F).  In this variant (`residued') the IIR part is driven by the
     _output_ of the FIR part.  This structure can be more accurate in
     signal modeling applications.

     INPUTS: B and A are vectors specifying the digital filter H(z) =
     B(z)/A(z).  Say `help filter' for documentation of the B and A
     filter coefficients.

     RETURNED:
        * R = column vector containing the filter-pole residues
        * P = column vector containing the filter poles
        * F = row vector containing the FIR part, if any
        * M = column vector of pole multiplicities

     EXAMPLES:
            Say `test residued verbose' to see a number of examples.

     For the theory of operation, see
     <http://ccrma.stanford.edu/~jos/filters/residued.html>

     See also: residue residued





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 72
Compute the partial fraction expansion (PFE) of filter H(z) = B(z)/A(z).



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
residuez


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 750
 -- Function File: [R, P, F, M] = residuez (B, A)
     Compute the partial fraction expansion of filter H(z) = B(z)/A(z).

     INPUTS: B and A are vectors specifying the digital filter H(z) =
     B(z)/A(z).  Say `help filter' for documentation of the B and A
     filter coefficients.

     RETURNED:
        * R = column vector containing the filter-pole residues
        * P = column vector containing the filter poles
        * F = row vector containing the FIR part, if any
        * M = column vector of pole multiplicities

     EXAMPLES:
            Say `test residuez verbose' to see a number of examples.

     For the theory of operation, see
     <http://ccrma.stanford.edu/~jos/filters/residuez.html>

     See also: residue residued





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 66
Compute the partial fraction expansion of filter H(z) = B(z)/A(z).



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 18
sampled2continuous


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 402
 Usage:
 
 xt = sampled2continuous( xn , T, t )
 
 Calculate the x(t) reconstructed
 from samples x[n] sampled at a rate 1/T samples
 per unit time.
 
 t is all the instants of time when you need x(t) 
 from x[n]; this time is relative to x[0] and not
 an absolute time.
 
 This function can be used to calculate sampling rate
 effects on aliasing, actual signal reconstruction
 from discrete samples.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 Usage:
 
 xt = sampled2continuous( xn , T, t )
 
 Calculate the x(t) reconstruc



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
sawtooth


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 664
 -- Function File: [Y] = sawtooth(T)
 -- Function File: [Y] = sawtooth(T,WIDTH)
     Generates a sawtooth wave of period `2 * pi' with limits `+1/-1'
     for the elements of T.

     WIDTH is a real number between `0' and `1' which specifies the
     point between `0' and `2 * pi' where the maximum is. The function
     increases linearly from `-1' to `1' in  `[0, 2 * pi * WIDTH]'
     interval, and decreases linearly from `1' to `-1' in the interval
     `[2 * pi * WIDTH, 2 * pi]'.

     If WIDTH is 0.5, the function generates a standard triangular wave.

     If WIDTH is not specified, it takes a value of 1, which is a
     standard sawtooth function.




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
Generates a sawtooth wave of period `2 * pi' with limits `+1/-1'  for
the elemen



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 7
sftrans


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 3640
 usage: [Sz, Sp, Sg] = sftrans(Sz, Sp, Sg, W, stop)

 Transform band edges of a generic lowpass filter (cutoff at W=1)
 represented in splane zero-pole-gain form.  W is the edge of the
 target filter (or edges if band pass or band stop). Stop is true for
 high pass and band stop filters or false for low pass and band pass
 filters. Filter edges are specified in radians, from 0 to pi (the
 nyquist frequency).

 Theory: Given a low pass filter represented by poles and zeros in the
 splane, you can convert it to a low pass, high pass, band pass or 
 band stop by transforming each of the poles and zeros individually.
 The following table summarizes the transformation:

 Transform         Zero at x                  Pole at x
 ----------------  -------------------------  ------------------------
 Low Pass          zero: Fc x/C               pole: Fc x/C
 S -> C S/Fc       gain: C/Fc                 gain: Fc/C 
 ----------------  -------------------------  ------------------------
 High Pass         zero: Fc C/x               pole: Fc C/x
 S -> C Fc/S       pole: 0                    zero: 0
                   gain: -x                   gain: -1/x
 ----------------  -------------------------  ------------------------
 Band Pass         zero: b ħ sqrt(b^2-FhFl)   pole: b ħ sqrt(b^2-FhFl)
        S^2+FhFl   pole: 0                    zero: 0
 S -> C --------   gain: C/(Fh-Fl)            gain: (Fh-Fl)/C
        S(Fh-Fl)   b=x/C (Fh-Fl)/2            b=x/C (Fh-Fl)/2
 ----------------  -------------------------  ------------------------
 Band Stop         zero: b ħ sqrt(b^2-FhFl)   pole: b ħ sqrt(b^2-FhFl)
        S(Fh-Fl)   pole: ħsqrt(-FhFl)         zero: ħsqrt(-FhFl)
 S -> C --------   gain: -x                   gain: -1/x
        S^2+FhFl   b=C/x (Fh-Fl)/2            b=C/x (Fh-Fl)/2
 ----------------  -------------------------  ------------------------
 Bilinear          zero: (2+xT)/(2-xT)        pole: (2+xT)/(2-xT)
      2 z-1        pole: -1                   zero: -1
 S -> - ---        gain: (2-xT)/T             gain: (2-xT)/T
      T z+1
 ----------------  -------------------------  ------------------------

 where C is the cutoff frequency of the initial lowpass filter, Fc is
 the edge of the target low/high pass filter and [Fl,Fh] are the edges
 of the target band pass/stop filter.  With abundant tedious algebra,
 you can derive the above formulae yourself by substituting the
 transform for S into H(S)=S-x for a zero at x or H(S)=1/(S-x) for a
 pole at x, and converting the result into the form:

    H(S)=g prod(S-Xi)/prod(S-Xj)

 The transforms are from the references.  The actual pole-zero-gain
 changes I derived myself.

 Please note that a pole and a zero at the same place exactly cancel.
 This is significant for High Pass, Band Pass and Band Stop filters
 which create numerous extra poles and zeros, most of which cancel.
 Those which do not cancel have a "fill-in" effect, extending the 
 shorter of the sets to have the same number of as the longer of the
 sets of poles and zeros (or at least split the difference in the case
 of the band pass filter).  There may be other opportunistic
 cancellations but I will not check for them.

 Also note that any pole on the unit circle or beyond will result in
 an unstable filter.  Because of cancellation, this will only happen
 if the number of poles is smaller than the number of zeros and the
 filter is high pass or band pass.  The analytic design methods all
 yield more poles than zeros, so this will not be a problem.
 
 References: 

 Proakis & Manolakis (1992). Digital Signal Processing. New York:
 Macmillan Publishing Company.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 52
 usage: [Sz, Sp, Sg] = sftrans(Sz, Sp, Sg, W, stop)




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
sgolay


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1079
 F = sgolay (p, n [, m [, ts]])
   Computes the filter coefficients for all Savitzsky-Golay smoothing
   filters of order p for length n (odd). m can be used in order to
   get directly the mth derivative. In this case, ts is a scaling factor. 

 The early rows of F smooth based on future values and later rows
 smooth based on past values, with the middle row using half future
 and half past.  In particular, you can use row i to estimate x(k)
 based on the i-1 preceding values and the n-i following values of x
 values as y(k) = F(i,:) * x(k-i+1:k+n-i).

 Normally, you would apply the first (n-1)/2 rows to the first k
 points of the vector, the last k rows to the last k points of the
 vector and middle row to the remainder, but for example if you were
 running on a realtime system where you wanted to smooth based on the
 all the data collected up to the current time, with a lag of five
 samples, you could apply just the filter on row n-5 to your window
 of length n each time you added a new sample.

 Reference: Numerical recipes in C. p 650

 See also: sgolayfilt



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 F = sgolay (p, n [, m [, ts]])
   Computes the filter coefficients for all Savi



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 10
sgolayfilt


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 857
 y = sgolayfilt (x, p, n [, m [, ts]])
    Smooth the data in x with a Savitsky-Golay smoothing filter of
    polynomial order p and length n, n odd, n > p.  By default, p=3
    and n=p+2 or n=p+3 if p is even.

 y = sgolayfilt (x, F)
    Smooth the data in x with smoothing filter F computed by sgolay.

 These filters are particularly good at preserving lineshape while
 removing high frequency squiggles. Particularly, compare a 5 sample
 averager, an order 5 butterworth lowpass filter (cutoff 1/3) and
 sgolayfilt(x, 3, 5), the best cubic estimated from 5 points:

    [b, a] = butter(5,1/3);
    x=[zeros(1,15), 10*ones(1,10), zeros(1,15)];
    plot(sgolayfilt(x),"r;sgolayfilt;",...
         filtfilt(ones(1,5)/5,1,x),"g;5 sample average;",...
         filtfilt(b,a,x),"c;order 5 butterworth;",...
         x,"+b;original data;");

 See also: sgolay



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 y = sgolayfilt (x, p, n [, m [, ts]])
    Smooth the data in x with a Savitsky-



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
shanwavf


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 97
 -- Function File: [PSI,X] = shanwavf (LB,UB,N,FB,FC)
     Compute the Complex Shannon wavelet.




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 36
Compute the Complex Shannon wavelet.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 13
sigmoid_train


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 555
 -- Function File: Y = sigmoid_train(T, RANGES, RC)
     Evaluates a train of sigmoid functions at T.

     The number and duration of each sigmoid is determined from RANGES.
     Each row of RANGES represents a real interval, e.g. if sigmod `i'
     starts at `t=0.1' and ends at `t=0.5', then `RANGES(i,:) = [0.1
     0.5]'.  The input RC is a array that defines the rising and
     falling time constants of each sigmoids. Its size must equal the
     size of RANGES.

     Run `demo sigmoid_train' to some examples of the use of this
     function.





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 44
Evaluates a train of sigmoid functions at T.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
sos2tf


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 808
 -- Function File: [B, A] = sos2tf (SOS, BSCALE)
     Convert series second-order sections to direct form H(z) =
     B(z)/A(z).

     INPUTS:
        * SOS = matrix of series second-order sections, one per row:
          SOS = [B1.' A1.'; ...; BN.' AN.'], where
          `B1.'==[b0 b1 b2] and A1.'==[1 a1 a2]' for section 1, etc.
          b0 must be nonzero for each section.
          See `filter()' for documentation of the second-order
          direct-form filter coefficients Bi and Ai.

        * BSCALE is an overall gain factor that effectively scales the
          output B vector (or any one of the input Bi vectors).

     RETURNED: B and A are vectors specifying the digital filter H(z) =
     B(z)/A(z).  See `filter()' for further details.

     See also: tf2sos zp2sos sos2pz zp2tf tf2zp





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 69
Convert series second-order sections to direct form H(z) = B(z)/A(z).



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
sos2zp


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1009
 -- Function File: [Z, P, G] = sos2zp (SOS, BSCALE)
     Convert series second-order sections to zeros, poles, and gains
     (pole residues).

     INPUTS:
        * SOS = matrix of series second-order sections, one per row:
          SOS = [B1.' A1.'; ...; BN.' AN.'], where
          `B1.'==[b0 b1 b2] and A1.'==[1 a1 a2]' for section 1, etc.
          b0 must be nonzero for each section.  See `filter()' for
          documentation of the second-order direct-form filter
          coefficients Bi and Ai.

        * BSCALE is an overall gain factor that effectively scales any
          one of the input Bi vectors.

     RETURNED:
        *   Z = column-vector containing all zeros (roots of B(z))
        *   P = column-vector containing all poles (roots of A(z))
        *   G = overall gain = B(Inf)

     EXAMPLE:
            [z,p,g] = sos2zp([1 0 1, 1 0 -0.81; 1 0 0, 1 0 0.49])
            => z = [i; -i; 0; 0], p = [0.9, -0.9, 0.7i, -0.7i], g=1

     See also: zp2sos sos2tf tf2sos zp2tf tf2zp





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
Convert series second-order sections to zeros, poles, and gains (pole
residues).



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
specgram


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 4570
 usage: [S [, f [, t]]] = specgram(x [, n [, Fs [, window [, overlap]]]])

 Generate a spectrogram for the signal. This chops the signal into
 overlapping slices, windows each slice and applies a Fourier
 transform to determine the frequency components at that slice.

 x: vector of samples
 n: size of fourier transform window, or [] for default=256
 Fs: sample rate, or [] for default=2 Hz
 window: shape of the fourier transform window, or [] for default=hanning(n)
    Note: window length can be specified instead, in which case
    window=hanning(length)
 overlap: overlap with previous window, or [] for default=length(window)/2

 Return values
    S is complex output of the FFT, one row per slice
    f is the frequency indices corresponding to the rows of S.
    t is the time indices corresponding to the columns of S.
    If no return value is requested, the spectrogram is displayed instead.

 Example
    x = chirp([0:0.001:2],0,2,500);  # freq. sweep from 0-500 over 2 sec.
    Fs=1000;                  # sampled every 0.001 sec so rate is 1 kHz
    step=ceil(20*Fs/1000);    # one spectral slice every 20 ms
    window=ceil(100*Fs/1000); # 100 ms data window
    specgram(x, 2^nextpow2(window), Fs, window, window-step);

    ## Speech spectrogram
    [x, Fs] = auload(file_in_loadpath("sample.wav")); # audio file
    step = fix(5*Fs/1000);     # one spectral slice every 5 ms
    window = fix(40*Fs/1000);  # 40 ms data window
    fftn = 2^nextpow2(window); # next highest power of 2
    [S, f, t] = specgram(x, fftn, Fs, window, window-step);
    S = abs(S(2:fftn*4000/Fs,:)); # magnitude in range 0<f<=4000 Hz.
    S = S/max(S(:));           # normalize magnitude so that max is 0 dB.
    S = max(S, 10^(-40/10));   # clip below -40 dB.
    S = min(S, 10^(-3/10));    # clip above -3 dB.
    imagesc(t, f, flipud(log(S)));   # display in log scale

 The choice of window defines the time-frequency resolution.  In
 speech for example, a wide window shows more harmonic detail while a
 narrow window averages over the harmonic detail and shows more
 formant structure. The shape of the window is not so critical so long
 as it goes gradually to zero on the ends.

 Step size (which is window length minus overlap) controls the
 horizontal scale of the spectrogram. Decrease it to stretch, or
 increase it to compress. Increasing step size will reduce time
 resolution, but decreasing it will not improve it much beyond the
 limits imposed by the window size (you do gain a little bit,
 depending on the shape of your window, as the peak of the window
 slides over peaks in the signal energy).  The range 1-5 msec is good
 for speech.

 FFT length controls the vertical scale.  Selecting an FFT length
 greater than the window length does not add any information to the
 spectrum, but it is a good way to interpolate between frequency
 points which can make for prettier spectrograms.

 After you have generated the spectral slices, there are a number of
 decisions for displaying them.  First the phase information is
 discarded and the energy normalized:

     S = abs(S); S = S/max(S(:));

 Then the dynamic range of the signal is chosen.  Since information in
 speech is well above the noise floor, it makes sense to eliminate any
 dynamic range at the bottom end.  This is done by taking the max of
 the magnitude and some minimum energy such as minE=-40dB. Similarly,
 there is not much information in the very top of the range, so
 clipping to a maximum energy such as maxE=-3dB makes sense:

     S = max(S, 10^(minE/10)); S = min(S, 10^(maxE/10));

 The frequency range of the FFT is from 0 to the Nyquist frequency of
 one half the sampling rate.  If the signal of interest is band
 limited, you do not need to display the entire frequency range. In
 speech for example, most of the signal is below 4 kHz, so there is no
 reason to display up to the Nyquist frequency of 10 kHz for a 20 kHz
 sampling rate.  In this case you will want to keep only the first 40%
 of the rows of the returned S and f.  More generally, to display the
 frequency range [minF, maxF], you could use the following row index:

     idx = (f >= minF & f <= maxF);

 Then there is the choice of colormap.  A brightness varying colormap
 such as copper or bone gives good shape to the ridges and valleys. A
 hue varying colormap such as jet or hsv gives an indication of the
 steepness of the slopes.  The final spectrogram is displayed in log
 energy scale and by convention has low frequencies on the bottom of
 the image:

     imagesc(t, f, flipud(log(S(idx,:))));



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 74
 usage: [S [, f [, t]]] = specgram(x [, n [, Fs [, window [, overlap]]]])




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
square


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 379
 -- Function File: S = square(T, DUTY)
 -- Function File: S = square(T)
     Generate a square wave of period 2 pi with limits +1/-1.

     If DUTY is specified, the square wave is +1 for that portion of
     the time.

                         on time
        duty cycle = ------------------
                     on time + off time

     See also: cos, sawtooth, sin, tripuls





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 56
Generate a square wave of period 2 pi with limits +1/-1.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 5
ss2tf


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 355
 -- Function File: [NUM, DEN] = ss2tf (A, B, C, D)
     Conversion from transfer 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)





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 49
Conversion from transfer function to state-space.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 5
ss2zp


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 172
 -- Function File: [POL, ZER, K] = ss2zp (A, B, C, D)
     Converts a state space representation to a set of poles and zeros;
     K is a gain associated with the zeros.





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
Converts a state space representation to a set of poles and zeros; K is
a gain a



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
tf2sos


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1051
 -- Function File: [SOS, G] = tf2sos (B, A)
     Convert direct-form filter coefficients to series second-order
     sections.

     INPUTS: B and A are vectors specifying the digital filter H(z) =
     B(z)/A(z).  See `filter()' for documentation of the B and A filter
     coefficients.

     RETURNED: SOS = matrix of series second-order sections, one per
     row:
     SOS = [B1.' A1.'; ...; BN.' AN.'], where
     `B1.'==[b0 b1 b2] and A1.'==[1 a1 a2]' for section 1, etc.
     b0 must be nonzero for each section (zeros at infinity not
     supported).  BSCALE is an overall gain factor that effectively
     scales any one of the Bi vectors.

     EXAMPLE:
          B=[1 0 0 0 0 1];
          A=[1 0 0 0 0 .9];
          [sos,g] = tf2sos(B,A)

          sos =

             1.00000   0.61803   1.00000   1.00000   0.60515   0.95873
             1.00000  -1.61803   1.00000   1.00000  -1.58430   0.95873
             1.00000   1.00000  -0.00000   1.00000   0.97915  -0.00000

          g = 1

     See also: sos2tf zp2sos sos2pz zp2tf tf2zp





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 72
Convert direct-form filter coefficients to series second-order sections.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 5
tf2ss


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 518
 -- Function File: [A, B, C, D] = tf2ss (NUM, DEN)
     Conversion from transfer 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)

     The state space system matrices obtained from this function will
     be in observable companion form as Wolovich's Observable Structure
     Theorem is used.




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 49
Conversion from transfer function to state-space.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 5
tf2zp


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 241
 -- Function File: [ZER, POL, K] = tf2zp (NUM, DEN)
     Converts transfer functions to poles-and-zero representations.

     Returns the zeros and poles of the system defined by NUM/DEN.  K
     is a gain associated with the system zeros.




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 62
Converts transfer functions to poles-and-zero representations.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 3
tfe


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 381
 Usage:
   [Pxx,freq] = tfe(x,y,Nfft,Fs,window,overlap,range,plot_type,detrend)

     Estimate transfer function of system with input "x" and output "y".
     Use the Welch (1967) periodogram/FFT method.
     Compatible with Matlab R11 tfe and earlier.
     See "help pwelch" for description of arguments, hints and references
     --- especially hint (7) for Matlab R11 defaults.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 Usage:
   [Pxx,freq] = tfe(x,y,Nfft,Fs,window,overlap,range,plot_type,detrend)




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 10
tfestimate


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 284
 Usage:
   [Pxx,freq]=tfestimate(x,y,Nfft,Fs,window,overlap,range,plot_type,detrend)

     Estimate transfer function of system with input "x" and output "y".
     Use the Welch (1967) periodogram/FFT method.
     See "help pwelch" for description of arguments, hints and references.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
 Usage:
   [Pxx,freq]=tfestimate(x,y,Nfft,Fs,window,overlap,range,plot_type,detr



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
triang


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 277
 usage:  w = triang (L)

 Returns the filter coefficients of a triangular window of length L.
 Unlike the bartlett window, triang does not go to zero at the edges
 of the window.  For odd L, triang(L) is equal to bartlett(L+2) except
 for the zeros at the edges of the window.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 24
 usage:  w = triang (L)




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 7
tripuls


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 629
 usage: y = tripuls(t, w, skew)

 Generate a triangular pulse over the interval [-w/2,w/2), sampled at
 times t.  This is useful with the function pulstran for generating a
 series pulses.

 skew is a value between -1 and 1, indicating the relative placement
 of the peak within the width.  -1 indicates that the peak should be
 at -w/2, and 1 indicates that the peak should be at w/2.

 Example
   fs = 11025;  # arbitrary sample rate
   f0 = 100;    # pulse train sample rate
   w = 0.3/f0;  # pulse width 3/10th the distance between pulses
   auplot(pulstran(0:1/fs:4/f0, 0:1/f0:4/f0, 'tripuls', w), fs);

 See also: pulstran



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 32
 usage: y = tripuls(t, w, skew)




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
tukeywin


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 633
 -- Function File: W = tukeywin (L, R)
     Return the filter coefficients of a Tukey window (also known as the
     cosine-tapered window) of length L. R defines the ratio between
     the constant section and and the cosine section. It has to be
     between 0 and 1. The function returns a Hanning window for R egals
     0 and a full box for R egals 1. By default R is set to 1/2.

     For a definition of the Tukey window, see e.g. Fredric J. Harris,
     "On the Use of Windows for Harmonic Analysis with the Discrete
     Fourier Transform, Proceedings of the IEEE", Vol. 66, No. 1,
     January 1978, Page 67, Equation 38.




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
Return the filter coefficients of a Tukey window (also known as the
cosine-taper



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
upsample


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 372
 -- Function File: Y = upsample (X, N)
 -- Function File: Y = upsample (X, N, OFFSET)
     Upsample the signal, inserting n-1 zeros between every element.

     If X is a matrix, upsample every column.

     If OFFSET is specified, control the position of the inserted
     sample in the block of n zeros.

     See also: decimate, downsample, interp, resample, upfirdn





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 63
Upsample the signal, inserting n-1 zeros between every element.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 8
welchwin


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 696
 -- Function File: [W] = welchwin(L,C)
     Returns a row vector containing a Welch window, given by
     W(n)=1-(n/N-1)^2,   n=[0,1, ... L-1].  Argument L is the length of
     the window.  Optional argument C specifies a "symmetric" window
     (the default), or a "periodic" window.

     A symmetric window has zero at each end and maximum in the middle;
     L must be an integer larger than 2.  `if c=="symmetric", N=(L-1)/2'

     A periodic window wraps around the cyclic interval [0,1, ... L-1],
     and is intended for use with the DFT  (functions fft(),
     periodogram() etc).  L must be an integer larger than 1.  `if
     c=="periodic", N=L/2'.

     See also: blackman, kaiser





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
Returns a row vector containing a Welch window, given by
W(n)=1-(n/N-1)^2,   n=[



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
window


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 221
 -- Function File: W = window (F, N, OPTS)
     Create a N-point windowing from the function F. The function F can
     be for example `@blackman'. Any additional arguments OPT are
     passed to the windowing function.




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 47
Create a N-point windowing from the function F.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 5
wkeep


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 127
 -- Function File: [Y] = wkeep(X,L,OPT)
     Extract the elements of x of size l from the center, the right or
     the left.




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 75
Extract the elements of x of size l from the center, the right or the
left.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 4
wrev


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 121
 -- Function File: [Y] = wrev(X)
     Reverse the order of the element of the vector x.

     See also: flipud, fliplr





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 49
Reverse the order of the element of the vector x.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 5
xcorr


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 2493
 usage: [R, lag] = xcorr (X [, Y] [, maxlag] [, scale])

 Estimate the cross correlation R_xy(k) of vector arguments X and Y
 or, if Y is omitted, estimate autocorrelation R_xx(k) of vector X,
 for a range of lags k specified by  argument "maxlag".  If X is a
 matrix, each column of X is correlated with itself and every other
 column.

 The cross-correlation estimate between vectors "x" and "y" (of
 length N) for lag "k" is given by 
    R_xy(k) = sum_{i=1}^{N}{x_{i+k} conj(y_i),
 where data not provided (for example x(-1), y(N+1)) is zero.

 ARGUMENTS
  X       [non-empty; real or complex; vector or matrix] data

  Y       [real or complex vector] data
          If X is a matrix (not a vector), Y must be omitted.
          Y may be omitted if X is a vector; in this case xcorr
          estimates the autocorrelation of X.

  maxlag  [integer scalar] maximum correlation lag
          If omitted, the default value is N-1, where N is the
          greater of the lengths of X and Y or, if X is a matrix,
          the number of rows in X.

  scale   [character string] specifies the type of scaling applied
          to the correlation vector (or matrix). is one of:
    'none'      return the unscaled correlation, R,
    'biased'    return the biased average, R/N, 
    'unbiased'  return the unbiassed average, R(k)/(N-|k|), 
    'coeff'     return the correlation coefficient, R/(rms(x).rms(y)),
          where "k" is the lag, and "N" is the length of X.
          If omitted, the default value is "none".
          If Y is supplied but does not have the ame length as X,
          scale must be "none".
          
 RETURNED VARIABLES
  R       array of correlation estimates
  lag     row vector of correlation lags [-maxlag:maxlag]

  The array of correlation estimates has one of the following forms.
    (1) Cross-correlation estimate if X and Y are vectors.
    (2) Autocorrelation estimate if is a vector and Y is omitted,
    (3) If X is a matrix, R is an matrix containing the cross-
        correlation estimate of each column with every other column.
        Lag varies with the first index so that R has 2*maxlag+1
        rows and P^2 columns where P is the number of columns in X.
        If Rij(k) is the correlation between columns i and j of X
             R(k+maxlag+1,P*(i-1)+j) == Rij(k)
        for lag k in [-maxlag:maxlag], or
             R(:,P*(i-1)+j) == xcorr(X(:,i),X(:,j)).
        "reshape(R(k,:),P,P)" is the cross-correlation matrix for X(k,:).




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 56
 usage: [R, lag] = xcorr (X [, Y] [, maxlag] [, scale])




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
xcorr2


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1046
 -- Function File: C = xcorr2 (A)
 -- Function File: C = xcorr2 (A, B)
 -- Function File: C = xcorr2 (..., SCALE)
     Compute the 2D cross-correlation of matrices A and B.

     If B is not specified, it defaults to the same matrix as A, i.e.,
     it's the same as `xcorr(A, A)'.

     The optional argument SCALE, defines the type of scaling applied
     to the cross-correlation matrix (defaults to "none"). Possible
     values are:
        * "biased"

          Scales the raw cross-correlation by the maximum number of
          elements of A and B involved in the generation of any element
          of C.

        * "unbiased"

          Scales the raw correlation by dividing each element in the
          cross-correlation matrix by the number of products A and B
          used to generate that element

        * "coeff"

          Normalizes the sequence so that the largest cross-correlation
          element is identically 1.0.

        * "none"

          No scaling (this is the default).

     See also: conv2, corr2, xcorr





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 53
Compute the 2D cross-correlation of matrices A and B.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 4
xcov


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 630
 usage: [c, lag] = xcov (X [, Y] [, maxlag] [, scale])

 Compute covariance at various lags [=correlation(x-mean(x),y-mean(y))].

 X: input vector
 Y: if specified, compute cross-covariance between X and Y,
 otherwise compute autocovariance of X.
 maxlag: is specified, use lag range [-maxlag:maxlag], 
 otherwise use range [-n+1:n-1].
 Scale:
    'biased'   for covariance=raw/N, 
    'unbiased' for covariance=raw/(N-|lag|), 
    'coeff'    for covariance=raw/(covariance at lag 0),
    'none'     for covariance=raw
 'none' is the default.

 Returns the covariance for each lag in the range, plus an 
 optional vector of lags.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 55
 usage: [c, lag] = xcov (X [, Y] [, maxlag] [, scale])




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 12
zerocrossing


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 186
 -- Function File: X0 = zerocrossing (X, Y)
     Estimates the points at which a given waveform y=y(x) crosses the
     x-axis using linear interpolation.

     See also: fzero, roots





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 80
Estimates the points at which a given waveform y=y(x) crosses the
x-axis using l



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
zp2sos


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1179
 -- Function File: [SOS, G] = zp2sos (Z, P)
     Convert filter poles and zeros to second-order sections.

     INPUTS:
        *   Z = column-vector containing the filter zeros
        *   P = column-vector containing the filter poles
        *   G = overall filter gain factor

     RETURNED:
        * SOS = matrix of series second-order sections, one per row:
          SOS = [B1.' A1.'; ...; BN.' AN.'], where
          `B1.'==[b0 b1 b2] and A1.'==[1 a1 a2]' for section 1, etc.
          b0 must be nonzero for each section.
          See `filter()' for documentation of the second-order
          direct-form filter coefficients Bi and %Ai, i=1:N.

        * BSCALE is an overall gain factor that effectively scales any
          one of the Bi vectors.

     EXAMPLE:
            [z,p,g] = tf2zp([1 0 0 0 0 1],[1 0 0 0 0 .9]);
            [sos,g] = zp2sos(z,p,g)

          sos =
             1.0000    0.6180    1.0000    1.0000    0.6051    0.9587
             1.0000   -1.6180    1.0000    1.0000   -1.5843    0.9587
             1.0000    1.0000         0    1.0000    0.9791         0

          g =
              1

     See also: sos2pz sos2tf tf2sos zp2tf tf2zp





# name: <cell-element>
# type: sq_string
# elements: 1
# length: 56
Convert filter poles and zeros to second-order sections.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 5
zp2ss


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 556
 -- 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, in the form:
                    .
                    x = Ax + Bu
                    y = Cx + Du




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 43
Conversion from zero / pole to state space.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 5
zp2tf


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 326
 -- 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 must appear in conjugate pairs.

    K
          Real scalar (leading coefficient).




# name: <cell-element>
# type: sq_string
# elements: 1
# length: 46
Converts zeros / poles to a transfer function.



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 6
zplane


# name: <cell-element>
# type: sq_string
# elements: 1
# length: 1223
 usage: zplane(b [, a]) or zplane(z [, p])

 Plot the poles and zeros.  If the arguments are row vectors then they
 represent filter coefficients (numerator polynomial b and denominator
 polynomial a), but if they are column vectors or matrices then they
 represent poles and zeros.

 This is a horrid interface, but I didn't choose it; better would be
 to accept b,a or z,p,g like other functions.  The saving grace is
 that poly(x) always returns a row vector and roots(x) always returns
 a column vector, so it is usually right.  You must only be careful
 when you are creating filters by hand.

 Note that due to the nature of the roots() function, poles and zeros
 may be displayed as occurring around a circle rather than at a single
 point.

 The transfer function is
                         
        B(z)   b0 + b1 z^(-1) + b2 z^(-2) + ... + bM z^(-M)
 H(z) = ---- = --------------------------------------------
        A(z)   a0 + a1 z^(-1) + a2 z^(-2) + ... + aN z^(-N)

               b0          (z - z1) (z - z2) ... (z - zM)
             = -- z^(-M+N) ------------------------------
               a0          (z - p1) (z - p2) ... (z - pN)

 The denominator a defaults to 1, and the poles p defaults to [].



# name: <cell-element>
# type: sq_string
# elements: 1
# length: 43
 usage: zplane(b [, a]) or zplane(z [, p])