home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS1.DI$ / MAGFIT.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  8.6 KB  |  227 lines
  1. % function [sys,fit,extracons] = magfit(vmagdata,dim,wt)
  2. %
  3. %  MAGFIT fits a stable, minimum phase transfer function
  4. %  to magnitude data, VMAGDATA, with a supplied frequency
  5. %  domain weighting function, WT. VMAGDATA is a VARYING 
  6. %  matrix giving the magnitude and frequency points to be 
  7. %  matched, 
  8. %
  9. %     [mag,rowpoint,omega,err] = vunpck(vmagdata).
  10. %
  11. %  DIM is a 1 by 4 vector with dim= [hmax,htol,nmin,nmax].
  12. %  The output system, SYS, is a stable minimum phase SISO
  13. %  system that satisfies:
  14. %
  15. %             1/(1+ h/W) < g^2/mag^2 < 1+ h/W
  16. %
  17. %  where g is the magnitude of the frequency response of SYS 
  18. %  for frequencies omega, [W,rowpoint,omega,err] = vunpck(WT).
  19. %  The degree of SYS is N, where 0<=NMIN<=N<=NMAX is the minimum 
  20. %  degree to meet the specificaton with H<=HMAX. If no such n 
  21. %  exists then N=NMAX. The value of H is approximately minimized 
  22. %  between upper and lower bounds terminating when 
  23. %  hupper-hlower<=HTOL. FIT is the achieved value of H. 
  24. %
  25. % See Also: FITMAG, FITMAGLP, MUSYNFIT, and MUSYNFLP.
  26.  
  27.  
  28. % The method used is note that g^2 is the ratio of
  29. % polynomials in omega^2 ,  g^2=p/q, with p and q having no 
  30. % positive real roots.  A scaled version of the  following linear
  31. % program is then solved
  32. %  max epp, s.t.
  33. % epp/(W*(1+h/W)^2)+q/(qbar*(1+h/W))<= p/mag^2*qbar
  34. %      <= (1+h/W)*q/qbar-epp/W
  35. % (the coefs. of  p and q are named a and b in the code),
  36. % and W is initially ones if not specified and qbar is an estimate of q.
  37. % If epp>0 and p and q have no real positive roots then a solution 
  38. % exists.  If epp>0 but p or q have positive real roots then  
  39. % extra cutting plane constraints are incorporated, with additional
  40. % frequency points having to be within a factor, sqrt(exfact),
  41. % (exfact=16 in the code at present) of the adjacent given magnitude values.
  42. % The dual LP is solved via modifications to LP and LINP called
  43. % MFLINP and MFLP.
  44.  
  45. function [sys,fit,extracons] = magfit(vmagdata,dim,Wt)
  46.  
  47. if ((nargin>3)|(nargin<2)),
  48.    disp('usage: [sys,fit,extracons] = magfit(vmagdata,dim)')
  49.    return
  50.  elseif (nargin>=2 & size(dim)~=[1 4]),
  51.    disp('usage [sys,fit,extracons] = magfit(vmagdata,dim)')
  52.    return
  53.  end
  54.     
  55.  
  56. hmax=dim(1);htol=dim(2);nmin=dim(3);nmax=dim(4);
  57. [mattype,rows,cols,pts] = minfo(vmagdata);
  58. if mattype~='vary', error('data must be of type varying'),end
  59. [magdata,rowpoint,omega,err]=vunpck(vmagdata);
  60. [sz_mg1,sz_mg2]=size(magdata);
  61. if any(sign(magdata)-ones(sz_mg1, sz_mg2))|any(sign(omega)-ones(sz_mg1, sz_mg2)),
  62.        error('magnitude and frequency data must be positive'), end,
  63. om=omega.^2;mag=magdata.^2;
  64. if nargin==2, Wdata=ones(pts,1);Winv=Wdata; 
  65.    else Wt=vabs(Wt);
  66.        [Wdata,rowpoint,indv,err] = vunpck(Wt); Winv=Wdata.^(-1);
  67.   end % if nargin==2
  68. [om,index]=sort(om); mag=mag(index);
  69. %define average frequency value.
  70. n=nmin;
  71. if om(1)==0, omav=sqrt(om(2)*om(pts));
  72.               else omav=sqrt(om(1)*om(pts));end,
  73. % rescale frequency list for numerical reasons.
  74.  om=om/omav;bold=ones(n+1,1);
  75. maxmag=max(mag);  minmag=min(mag);
  76. %initialize h-iteration variables.
  77. h=hmax;hupper=hmax;hlower=0;
  78. if htol>=0.49*hmax, htol=0.49*hmax;end,
  79. %initialize matrices for the LP with n=nmin.
  80. foundn=0;extracons=0;
  81. M=(om*ones(1,n+1)).^(ones(pts,1)*[n:-1:0]);
  82. N=(mag*ones(1,n+1)).\M;
  83. extracons=0;exfact=16; %this factor gives the bounds on extra points >1.
  84. inc=floor((pts-1)/(2*n+1));
  85. startbasic=[2:2*inc:2*n*inc+2,pts+inc+2:2*inc:pts+2+(2*n+1)*inc];
  86. while((hupper-hlower)>htol*max(hupper/hmax,1)),
  87.   qbar=polyval(bold,om);
  88.   ophw=1+h*Winv;ophwinv=ophw.^(-1);
  89.   P=(ophw*ones(1,n)).*M(:,2:n+1);
  90.   R=(ophwinv*ones(1,n)).*M(:,2:n+1); 
  91.   Winv2=Winv.*(ophwinv.^2);    
  92.   if n==0; extraA=[];lngthexa=0; 
  93.      else extraA=[zeros(1,n) 1 zeros(1,n-1) -maxmag*100 0;
  94.                   zeros(1,n) -1 zeros(1,n-1) minmag/100 0];
  95.           lngthexa=2;
  96.      end; %if n==0
  97.   A=[ zeros(1,2*n+1) 1; -1 zeros(1,2*n+1); zeros(1,n) -1 zeros(1,n+1);
  98.       zeros(1,2*n) -1 0; extraA;
  99. [qbar*ones(1,2*n+1);qbar*ones(1,2*n+1)].\[ N -P ; -N  R],[Winv;Winv2]];
  100.   B=[zeros(1,2*n+1) 1 ];
  101.   C=[h; -minmag/100; zeros(2+lngthexa,1); 
  102.     [qbar;qbar].\[M(:,1).*ophw;-M(:,1)./ophw]];
  103. % scale A, B and C.
  104.  % D1=max(abs([A C]'));
  105. %  A=(D1'*ones(1,2*n+2)).\A;
  106. %  C=D1'.\C;
  107.   D2=max(abs([A;B])); 
  108.   A=A./(ones(2*(pts+extracons)+4+lngthexa,1)*D2);
  109.   B=B./D2;
  110.   [basic,sol,epp,lambda,tnpiv,flopcount] = mflinp(A',B',C',startbasic);
  111.   x=D2'.\(A(basic,:)\C(basic));
  112.   oldh=h;
  113.   if epp<=0,  
  114.       if foundn, hlower=h;h=(0.5*hlower+0.5*hupper);
  115.          else
  116.            if n==nmax, 
  117.              hlower=h; h=2*h;hupper=h;
  118.            else n=n+1;
  119.              M=[om.*M(:,1),M];
  120.              N=[mag.\M(:,1),N];
  121.              bold=ones(n+1,1); 
  122.            end, %if n==nmax
  123.           inc=floor((pts-1)/(2*n+1));
  124.           startbasic=[5:2*inc:2*n*inc+5,...
  125.               extracons+pts+inc+5:2*inc:extracons+pts+5+(2*n+1)*inc];
  126.         end,%if foundn
  127.      end, %if epp<=0,
  128.  
  129.  
  130.   if epp>0,   %must check that a and b do not have any positive real roots.
  131.     a=x(1:n+1);
  132.     b=[1; x(n+2:2*n+1)];
  133.     rootsb=roots(b);
  134.     rootsa=roots(a);
  135.     realposinxb=[];bnotpos=0;
  136.     for i=1:length(rootsb), if ((imag(rootsb(i))==0)&(real(rootsb(i))>0)),
  137.       bnotpos=1;  realposinxb=[realposinxb, i];
  138.       end,end, %for i=1:length(rootsb)
  139.     realposinxa=[];anotpos=0;
  140.     for i=1:length(rootsa), if ((imag(rootsa(i))==0)&(real(rootsa(i))>0)),
  141.       anotpos=1; realposinxa=[realposinxa, i];
  142.      end,end,%for i=1:length(rootsa)
  143.  
  144.     if (bnotpos),
  145.       pass=0;
  146.       realposb=sort(rootsb(realposinxb));
  147.       if (max(size(realposb))==1), newomb=realposb(1)*[.9 ;1;1.1];
  148.          else  rb1=realposb(1);rb2=realposb(2);
  149.             newomb=[(rb1^0.25)*(rb2^0.75);sqrt(rb1*rb2);(rb2^0.25)*(rb1^0.75)];
  150.         end, %if (max(size(realposb))==1)
  151.       placeb=(pts-sum(sign(om(1:pts)-newomb(2)*ones(pts,1))))/2;
  152.       if placeb==0, newmagb=mag(1);newW=(exfact-1)/h;
  153.          elseif placeb==pts, newmagb=mag(pts);newW=(exfact-1)/h;
  154.            else newmagb=sqrt(mag(placeb)*mag(placeb+1));
  155.                 newW= (exfact*max(mag(placeb),mag(placeb+1))/newmagb-1)/h;
  156.         end % if placeb==0,
  157.       lob=length(newomb);
  158.       tmp1=(newomb*ones(1,n+1)).^(ones(lob,1)*(n:-1:0));
  159.       M=[M;tmp1];
  160.       N=[N;newmagb\tmp1];
  161.       Winv=[Winv;newW*ones(lob,1)];
  162.       om=[om;newomb];
  163.       mag=[mag;newmagb*ones(lob,1)];
  164.       extracons=extracons+lob;
  165.     end, %if (bnotpos)
  166.  
  167.     if (anotpos),
  168.       pass=0;
  169.           realposa=sort(rootsa(realposinxa));
  170.           if (max(size(realposa))==1), newoma=realposa(1)*[.9 ;1;1.1];                      else   ra1=realposa(1);ra2=realposa(2);
  171.              newoma=[(ra1^0.25)*(ra2^0.75);sqrt(ra1*ra2);(ra2^0.25)*(ra1^0.75)];
  172.             end; %if (max(size(realposa))==1)
  173.           placea=(pts-sum(sign(om(1:pts)-newoma(2)*ones(pts,1))))/2;
  174.           if placea==0, newmaga=mag(1);newW=(exfact-1)/h;
  175.               elseif placea==pts, newmaga=mag(pts);newW=(exfact-1)/h;
  176.                else newmaga=sqrt(mag(placea)*mag(placea+1));
  177.                 newW= (exfact*max(mag(placea),mag(placea+1))/newmaga-1)/h;
  178.           end % placea==0
  179.       loa=length(newoma);
  180.       tmp1=(newoma*ones(1,n+1)).^(ones(loa,1)*(n:-1:0));
  181.       M=[M;tmp1];
  182.       N=[N;newmaga\tmp1];
  183.       Winv=[Winv;newW*ones(loa,1)];
  184.       om=[om;newoma];
  185.       mag=[mag;newmaga*ones(loa,1)];
  186.       extracons=extracons+loa; 
  187.     end, %if (anotpos)
  188.  
  189.  
  190.    if (~anotpos&~bnotpos),
  191.     foundn=1;
  192.     g=abs(polyval(a,om(1:pts))./polyval(b,om(1:pts)));
  193.     fit=max([((g(1:pts)./mag(1:pts)-1)./Winv(1:pts));...
  194.              ((mag(1:pts)./g(1:pts)-1)./Winv(1:pts))]);
  195.     hupper=fit; h=(0.5*hupper+0.5*hlower);
  196.     agood=a; bgood=b; startbasic=basic;bold=b;
  197.     end, %if (~anotpos&~bnotpos)
  198.    end, %if epp>0
  199.   if (h~=oldh & extracons~=0),
  200.     Winv(pts+1:pts+extracons)=(oldh/h)*Winv(pts+1:pts+extracons);
  201.     end, %if (h~=oldh &
  202.  
  203.  end,% while((hupper-hlower)
  204.  
  205.            
  206. sqrootsa=-sqrt(-roots(agood))*sqrt(omav);
  207. [absrootsa,inda]=sort(abs(sqrootsa));
  208. sqrootsa=sqrootsa(inda);
  209. sqrootsb=-sqrt(-roots(bgood))*sqrt(omav);
  210. [absrootsb,indb]=sort(abs(sqrootsb));
  211. sqrootsb=sqrootsb(indb);
  212.  
  213. if agood(1)==0, agood(1)=1e-12; end  %something has gone numerically wrong.
  214. dega=n;for i=1:n, if absrootsa(i)==0, dega=dega-1;end;end
  215. degb=n;for i=1:dega, if absrootsb(i)==0 , degb=degb-1;end;end
  216. num=real(sqrt(abs(agood(1)))*poly(sqrootsa(1:dega)));
  217. den=real(poly(sqrootsb(1:degb)));
  218. g=abs(polyval(num,sqrt(-1)*omega)./polyval(den,sqrt(-1)*omega));
  219. fit=max([(((g.^2)./(magdata.^2)-1).*Wdata);...
  220.              (((magdata.^2)./(g.^2)-1).*Wdata)]);
  221.  
  222. if n>0, sys=nd2sys(num,den);
  223.    else sys=  num/den;
  224.  end
  225. %
  226. % Copyright MUSYN INC 1991,  All Rights Reserved
  227.