home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 7.ddi / ROBUST.DI$ / FITD.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  5.4 KB  |  193 lines
  1. function [aw,bw,cw,dw,magfit] = fitd(threl, freq, n, blksz,flag)
  2. % [SS_D,LOGDFIT] = FITD(LOGD,W,N,BLKSZ,FLAG)
  3. % [AD,BD,CD,DD,LOGDFIT] = FITD(LOGD,W,N,BLKSZ,FLAG) produces a stable,
  4. %   minimum-phase state-space realization SS_D of a diagonal transfer
  5. %   function matrix such that the diagonal elements' magnitude bode plots 
  6. %   approximately fit Bode magnitude plot data given in the rows of the
  7. %   matrix LOGD.  Uses the the numerically stable routine "YULEWALK.M".
  8. %
  9. %   Input:     LOGD  matrix whose rows are logs of magnitude bode plots;
  10. %              W frequency vector.
  11. %   Optional:  N vector of desired orders of the approximants (default=0);
  12. %              BLKSZ vector of sizes diagonal blocks in SS_D (default = 1);
  13. %              FLAG (default: display Bode plot; 0: no Bode plot);
  14. %              Bode plot of the fit will be displayed 4 at a time
  15. %              with a "pause" in between every full screen display.
  16. %   Output:    SS_D is a state-space realization of the diagonal transfer
  17. %              function matrix  D(s) =diag(d1(s)I1,...,dn(s)In) where
  18. %              di(s) is the fit to the i-th row of LOGD and I1,...,In are
  19. %              identity matrices of dimensions given in the n-vector BLKSZ.
  20. %              LOGDFIT contains the fit to LOGD, i.e., SS_D's bode plot.
  21.  
  22. % by Weizheng Wang (University of Southern California)
  23. % with enhancements by
  24. % R. Y. Chiang & M. G. Safonov 
  25. % Copyright (c) 1988 by the MathWorks, Inc.
  26. % All Rights Reserved.
  27. % ------------------------------------------------------------------------
  28.  
  29. [mmm,nnn] = size(threl);
  30. if mmm > nnn
  31.    threl = threl';
  32.    [mmm,nnn] = size(threl);
  33. end
  34. thre=exp(threl);  
  35. threl=threl/log(10);     % from log to log10
  36. freql = log10(freq);
  37.  
  38. bodeflag=0;
  39. magfit=[];
  40. if nargout==2 | nargout==5, bodeflag == 1; end
  41.  
  42. nfreq=max(size(freq));
  43. if nnn ~= nfreq,
  44.    error('MAG and W must have the same number of points')
  45. end
  46.  
  47. if nargin<3, n=zeros(1,mmm); end
  48. if max(size(n))==1, n=n*ones(1,mmm); end
  49.  
  50. if nargin <5
  51.    bodeflag=1;% compute bode plot of fit
  52.    flag = 1;  % display bode plot of fit
  53.    if nargin<4,blksz=ones(1,mmm);end
  54. end
  55. if max(size(n))==0, n=zeros(1,mmm); end
  56. if max(size(blksz))==0, blksz=ones(1,mmm); end
  57. if max(size(flag))==0,
  58.    bodeflag=1;% compute bode plot of fit
  59.    flag = 1;  % display bode plot of fit
  60. end
  61.  
  62.  
  63. % if scalar blksz, make it a vector
  64. if max(size(blksz))==1, blksz=blksz*ones(1,mmm);end
  65.  
  66. % Find the best sampling time "ff so that the frequency mapping from "s"
  67. % to "z" is approximately linear between 0 and Pi:
  68. %
  69. mx=max(freql);  mn=min(freql);  
  70. ff=10^(mn+2);
  71. ff=ff/3^(4-mx+mn); % ff is approx. sqrt(max(freq)/min(freq));
  72.  
  73. %
  74. % Bilinear mapping of the frequency axis from "z" to "s":
  75. %
  76. z=[0:126]/126*pi;
  77. j=sqrt(-1);
  78. w=imag(ff*(exp(j*z)-1)./(exp(j*z)+1));
  79. wl(1)=log10(freq(1));
  80. [x,y]=size(freq);
  81. x=max(x,y);
  82. wl(127)=log10(freq(x));
  83. wl(2:126)=log10(w(2:126));
  84. for i=2:126, if wl(i)>wl(127), wl(i)=wl(127); end; end;
  85. for i=2:126, if wl(i)<wl(1),  wl(i)=wl(1); end; end;
  86. z=z/pi;
  87. % 
  88.  
  89.  
  90. % Curve fit the continuous THREL magnitude data, one row at a time:
  91. %
  92. % Now prepare to curve fit THREL magnitude data, one row at a time:
  93.  
  94. loop = mmm;
  95.  
  96. if loop ~= max(size(n)),
  97.    error('Number of rows of MAG must equal the dimension of the vector N')
  98. end
  99. count = 0;
  100.  
  101. for i = 1:loop
  102.   mag =thre(i,:);
  103.   magl = threl(i,:);
  104.   nl = n(1,i);
  105.   count = count + 1;
  106. %      
  107.  
  108.   if nl==0
  109.     ad=1;
  110.     bd = exp(log(10)*polyfit(freql, magl, 0));
  111.   else
  112.   
  113.   % First do a high order least squares polynomial curve fit to the
  114.   % loglog bode magnitude plot data in ith row of THREL:
  115.     y = nl*2;
  116.     if y < 8, y = 8; end
  117.     if y>nfreq/2, y=round(nfreq/2); end
  118.     if y<1, y=1; end; if nfreq>100, y=10; end;
  119.     p = polyfit(freql, magl, y);
  120.  
  121.   % 
  122.   % Now obtain a lower order approximation in z-plane via YULEWALK.M:
  123.   %
  124.     dth = polyval(p, wl);
  125.     dth = exp(log(10)*dth);
  126.     dth(1)=mag(1);
  127.     dth(127)=mag(x);    
  128.     [bd,ad]=yulewalk(nl,z,dth);
  129.     % make sure the fitting is minimum phase
  130.     bd = polystab(bd);
  131.   end
  132.  
  133. % Discrete time state-space (az,bz,cz,dz):
  134. %
  135.   [az,bz,cz,dz]=tf2ss(bd,ad);
  136. %
  137. % Convert the discrete state-space to W-plane:
  138. %
  139.   if nl>0,
  140.     [aw0,bw0,cw0,dw0]=bilin(az,bz,cz,dz,-1,'Tustin',2/ff);
  141.   else
  142.     aw0=[];bw0=[];cw0=[];dw0=dz;
  143.   end
  144.   for j=1:blksz(i);
  145.      if i == 1 & j==1
  146.         aw = aw0; bw = bw0; cw = cw0; dw = dw0;
  147.      else
  148.         [aw,bw,cw,dw] = append(aw,bw,cw,dw,aw0,bw0,cw0,dw0);
  149.      end
  150.   end
  151.  
  152. %
  153. % Compute and display (if flag==1) the frequency response and fit
  154. %
  155.   if bodeflag==1,
  156.      if nl==0,
  157.          aw0=-1;bw0=0; cw0=0;
  158.      end
  159.      [ga_w0,ph0] = bode(aw0,bw0,cw0,dw0,1,freq);
  160.      magfit=[magfit; log(ga_w0')];
  161.      if flag==1,
  162.         if loop > 1
  163.            if count==1, clg; subplot(221); end;
  164.        if count==2, subplot(222); end;
  165.            if count==3, subplot(223); end;
  166.            if count==4, subplot(224); end;
  167.         end;
  168.         loglog(freq,mag,'x',freq,ga_w0);
  169.         title([ 'YULE-WALKER ORDER ' num2str(n(i)) ' FIT']);
  170.         xlabel('R/S (x: data; solid: fit)')
  171.         ylabel(['Mag D' num2str(i)]);
  172.      end
  173.   end
  174.   if count == 4     % reset the counter after four loops
  175.      count = 0;
  176.      pause          % pause the 4-window plot on the screen
  177.      clg
  178.   end
  179. %
  180. end  % of the loop
  181. %
  182. if nargout <3,
  183.     ss_w = mksys(aw,bw,cw,dw);
  184.     aw = ss_w;
  185.     bw=magfit;
  186. end
  187. %
  188. % -------- End of FITD.M % WW/RYC/MGS %
  189.  
  190.  
  191.  
  192.  
  193.