home *** CD-ROM | disk | FTP | other *** search
/ The Datafile PD-CD 5 / DATAFILE_PDCD5.iso / utilities / r / rlab / CTB / deconv < prev    next >
Encoding:
Text File  |  1995-11-15  |  1.5 KB  |  70 lines
  1. //-------------------------------------------------------------------------------
  2. //
  3. // deconv
  4. //
  5. // Syntax: A=deconv(b,a)
  6. //
  7. // This routine performs "de=convolution." If it is called as A=deconv(b,a)
  8. // then it deconvolves vector a out of vector b such that the following is true:
  9. //
  10. // b=conv(q,a) + r
  11. //
  12. // where q is the result of the deconvolution and b is the remainder.
  13. //
  14. // If b and a are vectors of polynomial coefficients, then using deconv
  15. // is equivalent to polynomial division (q is the quotient and r is the
  16. // remainder).
  17. //
  18. //      The results are returned in a list:
  19. //      A.q = results of a being deconvolved out of b.
  20. //      A.r = remainder
  21. //
  22. // Copyright (C), by Jeffrey B. Layton, 1994
  23. // Version JBL 940918
  24. //-------------------------------------------------------------------------------
  25.  
  26. rfile conv
  27.  
  28. deconv = function(b,a)
  29. {
  30.    local(nargs,mb,nb,na,q,r)
  31.  
  32. // Count number of input arguments
  33.    nargs=0;
  34.    if (exist(b)) {nargs=nargs+1;}
  35.    if (exist(a)) {nargs=nargs+1;}
  36.  
  37.    if (nargs != 2) {
  38.        error("DECONV: Wrong number of input arguments");
  39.    }
  40.  
  41.    mb=b.nr;
  42.    mb=b.nr;
  43.    nb=b.nc;
  44.    nb=a.nc;
  45.    nb=max(mb,nb);
  46.    na=length(a);
  47.  
  48.    if (na > nb) {
  49.        q=0;
  50.        r=b;
  51.        return << q=q; r=r >>
  52.    }
  53.  
  54.    if (a[1] == 0) {
  55.        error("First coefficient of A must be non-zero.")
  56.    }
  57.  
  58. // Deconvolution and polynomial division are the same operations
  59. // as a digital filter's impulse response B(z)/A(z):
  60.    q=filter(b, a, [1,zeros(1,nb-na)]);
  61.  
  62. // if (mb != 1) {
  63. //     q=q(:);
  64. // }
  65.  
  66.    r=b-conv(q,a);
  67.  
  68.    return << q=q; r=r>>
  69. };
  70.