home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Languages / RLaB 1.18c / help / filter < prev    next >
Encoding:
Text File  |  1994-11-06  |  2.5 KB  |  111 lines  |  [TEXT/RLAB]
  1.  
  2. Synopsis:    Digital filter implementation
  3.  
  4. Syntax:        filter ( B , A , X )
  5.         filter ( B , A , X , Zi )
  6.  
  7. Description:
  8.  
  9.     Filter is an implementation of the standard difference
  10.     equation:
  11.  
  12.     y[n] = b(1)*x[n] + b(2)*x[n-1] + ... b(nb+1)*x[n-nb]
  13.              - a(2)*y[n-1] - ... a(na+1)*y[n-na]
  14.  
  15.     The filter is implemented using a method described as a
  16.     "Direct Form II Transposed" filter. More for information see
  17.     Chapter 6 of "Discrete-Time Signal Processing" by Oppenheim
  18.     and Schafer.
  19.  
  20.     The inputs to filter are:
  21.  
  22.         B    The numerator coefficients, or zeros of the
  23.             system transfer function. The coefficients are
  24.             specified in a vector like:
  25.  
  26.             [ b(1) , b(2) , ... b(nb) ]
  27.  
  28.         A    The denominator coefficients, or the poles of
  29.             the system transfer function. the coefficients
  30.             are specified in a vector like:
  31.  
  32.             [ a(1) , a(2) , ... a(na) ]
  33.  
  34.         X    A vector of the filter inputs.
  35.  
  36.         Zi    [Optional] The initial delays of the filter.
  37.  
  38.     The filter outputs are in a list with element names:
  39.  
  40.         y    The filter output. y is a vector of the same
  41.             dimension as X.
  42.  
  43.         zf    A vector of the final values of the filter
  44.             delays.
  45.  
  46.     The A(1) coefficient must be non-zero, as the other
  47.     coefficients are divided by A(1).
  48.  
  49.     Below is an implementation of filter() in a r-file - it is
  50.     provided for informational usage only.
  51.  
  52.     #--------------------------------------------------------------
  53.     #  Simplistic version of RLaB's builtin function filter()
  54.     #  Y = filter ( b, a, x )
  55.     #  Y = filter ( b, a, x, zi )
  56.     #
  57.     
  58.     rfilter = function ( b , a , x , zi )
  59.     {
  60.       local ( b , a , x , zi )
  61.  
  62.       ntotal = x.nr * x.nc;
  63.       M = b.nr * b.nc;
  64.       N = a.nr * a.nc;
  65.       NN = max ([ M, N ]);
  66.       y = zeros (x.nr, x.nc); 
  67.     
  68.       # Fix up pole and zero vectors.
  69.       # Make them the same length, this makes
  70.       # filter's job much easier.
  71.     
  72.       if (N < NN) { a[NN] = 0; }
  73.       if (M < NN) { b[NN] = 0; }
  74.       
  75.       # Adjust filter coefficients
  76.       if (a[1] == 0) { error ("rfilter: 1st A term must be non-zero"); }
  77.       a[2:NN] = a[2:NN] ./ a[1];
  78.       b = b ./ a[1];
  79.     
  80.       # Create delay vectors and load inital delays.
  81.       # Add an extra term to vi[] to make filter's 
  82.       # job a little easier. This extra term will
  83.       # always be zero.
  84.     
  85.       v = zeros (NN, 1);
  86.       vi = zeros (NN+1, 1);
  87.     
  88.       if (exist (zi))
  89.       {
  90.         vi[1:NN] = zi;   
  91.       }
  92.     
  93.       #
  94.       # Do the work...
  95.       #
  96.     
  97.       for (n in 1:ntotal)
  98.       {
  99.         v[1] = b[1]*x[n] + vi[2];
  100.         y[n] = v[1];
  101.         for (k in 2:NN)
  102.         {
  103.           v[k] = b[k]*x[n] - a[k]*v[1] + vi[k+1];
  104.           vi[k] = v[k];
  105.         }
  106.       }
  107.     
  108.       return << y = y; zf = v >>;
  109.     };
  110.     #--------------------------------------------------------------
  111.