home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 8.ddi / SIGNAL.DI$ / INVFREQZ.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  4.6 KB  |  139 lines
  1. function [b,a]=invfreqz(g,w,nb,na,wf,maxiter,tol,pf)
  2. %INVFREQZ  Discrete filter least sqaures fit to frequency response data.
  3. %
  4. %    [B,A] = INVFREQZ(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]=INVFREQZ(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] = INVFREQZ(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]=INVFREQZ(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]=INVFREQZ(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 INVFREQS.
  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-19-92 Trace mode made optional TPK
  29. %    (c) Copyright 1986-92 by The MathWorks, Inc.
  30.  
  31. nk=0;T=1; % The code is prepared for arbitrary sampling interval T and for
  32.           % constraining the numerator to 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. if any(w>pi),disp('WARNING: W contains values over the Nyquist frequency!'),end
  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,nb+nk-1);
  54. OM=exp(-i*[0:nm]'*w*T);end
  55. %
  56. % Estimation in the least squares case:
  57. %
  58.     Dva=(OM(2:na+1,:).').*(g*ones(1,na));
  59.     Dvb=-(OM(nk+1:nk+nb,:).');
  60.     D=[Dva Dvb].*(wf*ones(1,na+nb));
  61.     R=real(D'*D);
  62.     Vd=real(D'*(-g.*wf));
  63.     th=R\Vd;
  64.     a=[1 th(1:na)'];b=[zeros(1,nk) th(na+1:na+nb)'];
  65.  
  66. if nargin<6,return,end
  67.  
  68. % Now for the iterative minimization
  69.  
  70. if nargin<7,tol=0.01;end
  71. indb=1:length(b);indg=1:length(a);
  72. a=polystab(a); % Stabilizing the denominator
  73.  
  74. % The initial estimate:
  75.  
  76.     GC=((b*OM(indb,:))./(a*OM(indg,:))).';
  77.     e=(GC-g).*wf;
  78.     Vcap=e'*e; t=[a(2:na+1) b(nk+1:nk+nb)]';
  79.     if (verb),
  80.         clc, disp(['  INITIAL ESTIMATE'])    
  81.         disp(['Current fit: ' num2str(Vcap)])
  82.         disp(['par-vector'])
  83.         disp(t)
  84.     end;
  85. %
  86. % ** the minimization loop **
  87. %
  88.     gndir=2*tol+1;l=0;st=0;
  89.     while [norm(gndir)>tol l<maxiter st~=1]
  90.         l=l+1;
  91.  
  92. %     * compute gradient *
  93.  
  94.         D31=(OM(2:na+1,:).').*(-GC./((a*OM(1:na+1,:)).')*ones(1,na));
  95.         D32=(OM(nk+1:nk+nb,:).')./((a*OM(1:na+1,:)).'*ones(1,nb));
  96.         D3=[D31 D32].*(wf*ones(1,na+nb));
  97. %     * compute Gauss-Newton search direction
  98.  
  99.         e=(GC-g).*wf;
  100.         R=real(D3'*D3);
  101.         Vd=real(D3'*e);
  102.         gndir=R\Vd;
  103.  
  104. %     * search along the gndir-direction *
  105.  
  106.         ll=0;,k=1;V1=Vcap+1;
  107.         while [V1 > Vcap ll<20],
  108.  
  109.             t1=t-k*gndir; if ll==19,t1=t;end
  110.                a=polystab([1 t1(1:na)']);
  111.             t1(1:na)=a(2:na+1)';   %Stabilizing denominator
  112.                b=[zeros(1,nk) t1(na+1:na+nb)'];
  113.                GC=((b*OM(indb,:))./(a*OM(indg,:))).';
  114.                V1=((GC-g).*wf)'*((GC-g).*wf); t1=[a(2:na+1) b(nk+1:nk+nb)]';
  115.             if (verb),
  116.                 home, disp(int2str(ll))
  117.             end;
  118.             k=k/2;
  119.             ll=ll+1; if ll==20, st=1;end
  120.             if ll==10,gndir=Vd/norm(R)*length(R);k=1;end
  121.         end
  122.  
  123.         if (verb),
  124.             home
  125.             disp(['      ITERATION # ' int2str(l)])
  126.             disp(['Current fit:  ' num2str(V1) '  Previous fit:  ',...
  127.                   num2str(Vcap)])
  128.             disp(['Current par prev.par   GN-dir'])
  129.             disp([t1 t gndir])
  130.             disp(['Norm of GN-vector: ' num2str(norm(gndir))])
  131.             if st==1, 
  132.                 disp(['No improvement of the criterion possible along ',...
  133.                       'the search direction']), 
  134.                 disp('Iterations therefore terminated'),
  135.             end
  136.         end
  137.         t=t1; Vcap=V1;
  138.     end
  139.