home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 8.ddi / SIGNAL.DI$ / INVFREQS.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  4.4 KB  |  145 lines
  1. function [b,a]=invfreqs(g,w,nb,na,wf,maxiter,tol,pf)
  2. %INVFREQS Analog filter least squares fit to frequency response data.
  3. %
  4. %    [B,A] = INVFREQS(H,W,nb,na) gives real numerator and denominator 
  5. %    coefficients B and A of orders nb and na respectively, where
  6. %    H is the desired complex frequency response of the system at frequency
  7. %    points W, and W contains the frequency values in radians/s.
  8. %
  9. %    [B,A]=INVFREQS(H,W,nb,na,Wt) allows the fit-errors to the weighted
  10. %    versus frequency.  LENGTH(Wt)=LENGTH(W)=LENGTH(H).
  11. %    Determined by minimization of sum |B-H*A|^2*Wt over the freqs in W.
  12. %
  13. %    [B,A] = INVFREQS(H,W,nb,na,Wt,ITER) does another type of fit:
  14. %    Sum |B/A-H|^2*Wt is minimized with respect to the coefficients in B and
  15. %    A by numerical search in at most ITER iterations.  The A-polynomial is 
  16. %    then constrained to be stable.  [B,A]=INVFREQS(H,W,nb,na,Wt,ITER,TOL)
  17. %    stops the iterations when the norm of the gradient is less than TOL.
  18. %    The default value of TOL is 0.01.  The default value of Wt is all ones.
  19. %    This default value is also obtained by Wt=[].
  20. %
  21. %     [B,A]=INVFREQS(H,W,nb,na,Wt,ITER,TOL,'trace') will provide a textual
  22. %     progress report of the iteration.
  23. %
  24. %    See also FREQZ, FREQS, and INVFREQZ.
  25.  
  26. %    J.O. Smith and J.N. Little 4-23-86
  27. %    Revised 4-27-88 JNL, Rewritten 9-21-92 Lennart Ljung
  28. %     10-22-92 Trace mode made optional TPK
  29. %    (c) Copyright 1986-92 by The MathWorks, Inc.
  30.  
  31. nk=0;     % The code is prepared for constraining the numerator to 
  32.       % begin with nk zeros.
  33.  
  34. nb=nb+nk+1;
  35. if nargin<5,wf=[];end
  36. if nargin<8,
  37.     verb=0;
  38. elseif (strcmp(pf,'trace')),
  39.     verb=1;
  40. else
  41.     error(['Trace flag ''' pf ''' not recognizeable']);
  42. end
  43. if isempty(wf),wf=ones(length(w),1);end
  44. wf=sqrt(wf);
  45.  
  46. if length(g)~=length(w),error('The lengths of H and W must coincide'),end
  47. if length(wf)~=length(w),error('The lengths of Wt and W must coincide'),end
  48.  
  49. [rw,cw]=size(w);if rw>cw w=w';end
  50. [rg,cg]=size(g);if cg>rg g=g.';end
  51. [rwf,cwf]=size(wf);if cwf>rwf wf=wf';end
  52.  
  53. i=sqrt(-1);nm=max(na+1,nb+nk);
  54. indb=nb:-1:1; indg=na+1:-1:1; inda=na:-1:1;
  55.  
  56.    OM=ones(1,length(w));
  57.    for kom=1:nm-1
  58.        OM=[OM;(i*w).^kom];
  59.    end
  60.  
  61. %
  62. % Estimation in the least squares case:
  63. %
  64.     Dva=(OM(inda,:).').*(g*ones(1,na));
  65.     Dvb=-(OM(indb,:).');
  66.     D=[Dva Dvb].*(wf*ones(1,na+nb));
  67.     R=real(D'*D);
  68.     Vd=real(D'*((-g.*OM(na+1,:).').*wf));
  69.     th=R\Vd;
  70.     a=[1 th(1:na)'];b=[zeros(1,nk) th(na+1:na+nb)'];
  71.  
  72. if nargin<6,return,end
  73.  
  74. % Now for the iterative minimization
  75.  
  76. if nargin<7,tol=0.01;end
  77. % Stabilizing the denominator:
  78.  
  79. if na>1,v=roots(a);vind=(real(v)>0);v(vind)=-v(vind);a=real(poly(v));end
  80.  
  81. % The initial estimate:
  82.  
  83.     GC=((b*OM(indb,:))./(a*OM(indg,:))).';
  84.     e=(GC-g).*wf;
  85.     Vcap=e'*e; t=[a(2:na+1) b(nk+1:nk+nb)]';
  86.     if (verb),
  87.         clc, disp(['  INITIAL ESTIMATE'])    
  88.         disp(['Current fit: ' num2str(Vcap)])
  89.         disp(['par-vector'])
  90.         disp(t)
  91.     end
  92. %
  93. % ** the minimization loop **
  94. %
  95. gndir=2*tol+1;l=0;st=0;
  96. while [norm(gndir)>tol l<maxiter st~=1]
  97.     l=l+1;
  98.  
  99. %     * compute gradient *
  100.  
  101.     D31=(OM(inda,:).').*(-GC./((a*OM(indg,:)).')*ones(1,na));
  102.     D32=(OM(indb,:).')./((a*OM(indg,:)).'*ones(1,nb));
  103.     D3=[D31 D32].*(wf*ones(1,na+nb));
  104. %     * compute Gauss-Newton search direction
  105.  
  106.     e=(GC-g).*wf;
  107.     R=real(D3'*D3);
  108.     Vd=real(D3'*e);
  109.     gndir=R\Vd;
  110.  
  111. %     * search along the gndir-direction *
  112.  
  113.     ll=0;,k=1;V1=Vcap+1;
  114.     while [V1 > Vcap ll<20],
  115.  
  116.         t1=t-k*gndir; if ll==19,t1=t;end
  117.      a=[1 t1(1:na)'];
  118.     b=[zeros(1,nk) t1(na+1:na+nb)'];
  119.         if na>1,v=roots(a);vind=(real(v)>0);v(vind)=-v(vind);a=real(poly(v));
  120.     end,t1(1:na)=a(2:na+1)'; % Stabilizing the denominator
  121.     GC=((b*OM(indb,:))./(a*OM(indg,:))).';
  122.     V1=((GC-g).*wf)'*((GC-g).*wf);
  123.         if (verb),
  124.         home, disp(int2str(ll))
  125.         end;
  126.     k=k/2;
  127.     ll=ll+1; if ll==10, gndir=Vd/norm(R)*length(R);k=1;end
  128.               if ll==20,st=1;end
  129.     end
  130.     if (verb),
  131.         home
  132.         disp(['      ITERATION # ' int2str(l)])
  133.         disp(['Current fit:  ' num2str(V1) '  Previous fit:  ' num2str(Vcap)])
  134.         disp(['Current par prev.par   GN-dir'])
  135.         disp([t1 t gndir])
  136.         disp(['Norm of GN-vector: ' num2str(norm(gndir))])
  137.         if st==1, 
  138.             disp(['No improvement of the criterion possible along the '...
  139.                    'search ''direction''']), 
  140.             disp('Iterations therefore terminated'),
  141.         end
  142.     end
  143.     t=t1; Vcap=V1;
  144. end
  145.