home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / FITGAIN.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.8 KB  |  104 lines
  1. function [a,b,c,d] = fitgain(mag,w,ord,wt,flag)
  2. % [SS_] = FITGAIN(MAG,W,ORD,WT,FLAG) or
  3. % [A,B,C,D] = FITGAIN(MAG,W,ORD,WT,FLAG) produces a "stable" state-space 
  4. %    realization of the continuous magnitude curve "MAG".
  5. %
  6. %    Inputs: 
  7. %               mag ---- absolute magnitude array
  8. %               w   ---- frequencies at which "mag" are evaluated (rad/sec)
  9. %               ord ---- size of the realization 
  10. %    Optional:  
  11. %               wt  ---- weighting of the curve fit (default = ones(w))
  12. %               flag --- display Bode plot (default), 0: no Bode plot
  13. %
  14. %    Outputs: (a,b,c,d) ---- stable realization of "MAG"
  15. %
  16. %  The algorithm uses the following 3-step procedure:
  17. %                         2
  18. %                   |G(s)|   = G(s) * G(-s)
  19. %    Step 1: take the PSD of G(s), i.e. PSD = |G(s)|^2.
  20. %    Step 2: use the MATLAB function "invfreqs.m" to fit the PSD
  21. %            with a rational transfer function.
  22. %    Step 3: Pull out the stable and minimum phase part of the realization.
  23. %                          (done !!)
  24.  
  25. % R. Y. Chiang & M. G. Safonov 
  26. % Copyright (c) 1988 by the MathWorks, Inc.
  27. % All Rights Reserved.
  28. % ------------------------------------------------------------------------
  29.  
  30. xyz = nargin;
  31. if xyz < 4
  32.    len = length(w);
  33.    wt = ones(1,len);
  34.    flag = 1;
  35. end
  36. %
  37. psd  = mag.*mag;
  38. [num,den] = invfreqs(psd,w,2*ord,2*ord,wt);
  39. k = 1;                           % strip out the leading zeros
  40. numk = num;
  41. while abs(num(1,k)) < 1.e-8
  42.       numk = num(1,k+1:2*ord+1);       
  43.       k = k+1;
  44. end
  45. %
  46. numk = numk/numk(1,1);           % make the polynomial monic
  47. a1 = sqrt(numk(1,1));
  48. [rk,ck] = size(numk);
  49. %
  50. if  ck > 1
  51.    z = roots(numk);
  52.    [rz,cz] = size(z);
  53.    zmin = [];
  54.    for s = 1:rz
  55.       if real(z(s,1)) < 0
  56.          zmin = [zmin;z(s,1)];
  57.       end
  58.    end
  59.    nums = real(poly(zmin));
  60. else
  61.    nums = numk;
  62. end
  63. %
  64. h = 1;                           % strip out the leading zeros
  65. denk = den;
  66. [rkd,ckd] = size(denk);
  67. while abs(den(1,h)) < 1.e-8
  68.       denk = den(1,h+1:2*ord+1);       
  69.       h = h+1;
  70. end
  71. %
  72. b1 = sqrt(denk(1,1));
  73. denk = denk/denk(1,1);           % make the polynomial monic
  74. %
  75. poo = roots(denk);
  76. [rp,cp] = size(poo);
  77. pstable = [];
  78. for q = 1:rp
  79.       if real(poo(q,1)) < 0
  80.          pstable = [pstable;poo(q,1)];
  81.       end
  82. end
  83. %
  84. dens = real(poly(pstable));
  85. if length(nums) > length(dens)
  86.    error('Ill-conditioned result, try more frequency points ..');
  87. end
  88. [a,b,c,d] = tf2ss(a1*nums,b1*dens);
  89. %
  90. if nargout == 1
  91.    ss_ = mksys(a,b,c,d);
  92.    a = ss_;
  93. end
  94. %
  95. if flag == 1
  96.    [magfit,ph] = bode(a,b,c,d,1,w);
  97.    loglog(w,mag,'x',w,magfit);
  98.    title('Equation Error Method')
  99.    xlabel('R/S (x: data; solid: fit)')
  100.    ylabel('Log(Mag)');
  101. end
  102. %
  103. % ------ End of FITGAIN.M % RYC/MGS %
  104.