home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS1.DI$ / FITSYS.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.6 KB  |  89 lines
  1. % function sys = fitsys(resp,ord,wt,code)
  2. %
  3. %   FITSYS fits frequency response data in RESP
  4. %   with a transfer function of order ORD, using a
  5. %   frequency dependent weight in WT (optional).
  6. %
  7. %   To normalize the effect of high freq .vs. low freq
  8. %   on the fit, the program appends an additional
  9. %   frequency dependent weighting to the user supplied
  10. %   weighting WT. The correction is equal to (1/omega.^ord)
  11. %   at each frequency point, omega.
  12. %
  13. %   The fourth argument, CODE, is optional. If set
  14. %   to 0 (default), then the fit is returned unchanged
  15. %   from the INVFREQS function return. If CODE = 1, 
  16. %   as in the mu-synthesis routines, it forces the fit
  17. %   to be stable, minimum phase, simply by reflecting the 
  18. %   poles and zeros if necessary. Typically in this case, 
  19. %   the response RESP comes from the program GENPHASE,
  20. %   and already corresponds to a stable, minimum phase
  21. %   transfer function.
  22. %
  23. %  See also: FITMAG, MAGFIT, GENPHASE, INVFREQS, MUSYNFIT, and MUSYNFLP.
  24.  
  25. %change has been made by Wes Wang at the MathWorks in 5/12/92
  26. %
  27. function sys = fitsys(resp,ord,wt,code)
  28.  if nargin <= 1
  29.    disp(['usage: sys = fitsys(resp,ord,wt,code)']);
  30.    return
  31.  end
  32.  [rtype,rrows,rcols,rnum] = minfo(resp);
  33.  gjw    =  (resp(1:rnum,1)).';        % row vector
  34.  omega  =  (resp(1:rnum,2)).';        % row vector
  35.  if ord == 0
  36.    ordwt = 1;
  37.  else
  38.    ordwt = nd2sys([1],[1 zeros(1,ord)]);
  39.  end
  40.  ordwt_g = vabs(frsp(ordwt,resp));
  41.  wtjw = (ordwt_g(1:rnum,1)).';        % row vector
  42.  
  43.  if nargin == 2
  44.    [b,a] = invfreqs(gjw,omega,ord,ord,wtjw);
  45.    code = 0;
  46.  elseif nargin >= 3
  47.    [wttype,wtrows,wtcols,wtnum] = minfo(wt);
  48.    if wttype == 'vary'
  49.      if wtrows == 1 & wtcols == 1
  50.        wtjw = abs((wt(1:wtnum,1))).' .* wtjw ;
  51.        [b,a] = invfreqs(gjw,omega,ord,ord,wtjw);
  52.      else
  53.        error('weighting function must be 1x1 VARYING matrix')
  54.        return
  55.      end
  56.      if nargin == 3
  57.        code = 0;
  58.      end
  59.    elseif wttype == 'cons'
  60.      [b,a] = invfreqs(gjw,omega,ord,ord,wtjw);
  61. %%%%%%% Change has made here --- Wes Wang %%%%%%%%
  62. % original code
  63. %     code = wt;
  64. %%%%%%% change to be
  65.      code = 0;
  66.    else
  67.      error('third argument shouldn''t be a SYSTEM');
  68.      return
  69.    end
  70.  end
  71.  if code == 1           
  72.    b = refroots(b);
  73.    a = refroots(a);
  74.  elseif code ~= 0
  75.    error('last argument must be 0 or 1')
  76.    return
  77.  end
  78.  if ord == 0
  79.    sys = b/a;
  80.  else
  81.    if max(real(roots(a)))<0
  82.      sys = sysbal(nd2sys(b,a));
  83.    else
  84.      sys = nd2sys(b,a);
  85.    end
  86.  end
  87. %
  88. % Copyright MUSYN INC 1991,  All Rights Reserved
  89.