home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS1.DI$ / DYPERT.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.9 KB  |  105 lines
  1. %function rat_pert = dypert(pvec,blk,bnds,blks)
  2. %
  3. %  Creates rational, stable perturbation which
  4. %  interpolates the frequency varying perturbation, PVEC,
  5. %  at the peak value of the lower bound in BNDS (frequency
  6. %  where the perturbation is the smallest). PVEC is
  7. %  perturbation vector generated from MU, BLK is the
  8. %  block structure, BNDS is the upper and lower bounds
  9. %  generated from MU, BLKS is the optional numbers of
  10. %  the blocks which are to be included in the rational
  11. %  perturbation output.
  12. %
  13. %   See also: MU, UNWRAPD, UNWRAPP, and SISORAT.
  14.  
  15. function rat_pert = dypert(pertrv,blk,bnds,blks)
  16.   if nargin < 3
  17.     disp('usage: sys = dypert(row_p,blk,bnds)')
  18.     return
  19.   end
  20.   [nblk,dum] = size(blk);
  21.   if nargin == 3
  22.     blks = 1:nblk;
  23.   end
  24.   [btype,brows,bcols,bnum] = minfo(bnds);
  25.   [ptype,prows,pcols,pnum] = minfo(pertrv);
  26.   code = 0;
  27.   if code == 1
  28.     error('different # of frequency points')
  29.     return
  30.   elseif code == 2
  31.     error('bounds and perturbations need to be frequency varying')
  32.     return
  33.   elseif code == 3
  34.     error('frequency points are different')
  35.     return
  36.   end
  37. %
  38.   [peak,peakloc] = max(bnds(1:bnum,2));
  39. %
  40.   rowp = pertrv(peakloc,1:pcols);
  41.   omegap = pertrv(peakloc,pcols+1);
  42.   sys = [];
  43.   loc = 1;
  44.   for i=1:nblk
  45.     if blk(i,2) == 0 | (blk(i,1)*blk(i,2) == 1)
  46.       if any(blks==i)
  47.         gain = rowp(loc);
  48.         block = [];
  49.         delta = sisorat(omegap,gain);
  50.         for j=1:blk(i,1)
  51.           block = daug(block,delta);
  52.         end
  53.         sys = daug(sys,block);
  54.       end
  55.       loc = loc + 1;
  56.     else
  57.       if blk(i,1) == 1
  58.         if any(blks==i)
  59.           block = [];
  60.           for j=1:blk(i,2)
  61.             block = sbs(block,sisorat(omegap,rowp(loc)));
  62.             loc = loc + 1;
  63.           end
  64.           sys = daug(sys,block);
  65.         else
  66.           loc = loc + blk(i,2);
  67.         end
  68.       else
  69.         if blk(i,2) == 1
  70.           if any(blks==i)
  71.             block = [];
  72.             for j=1:blk(i,1)
  73.               block = abv(block,sisorat(omegap,rowp(loc)));
  74.               loc = loc + 1;
  75.             end
  76.             sys = daug(sys,block);
  77.           else
  78.             loc = loc + blk(i,1);
  79.           end
  80.         else
  81.           if any(blks==i)
  82.             tmp = reshape(rowp(loc:loc-1+(blk(i,1)*blk(i,2))),blk(i,2),blk(i,1));
  83.             [u,s,v]=svd(tmp.');
  84.             col = u(:,1)*sqrt(s(1));
  85.             row = (v(:,1))'*sqrt(s(1));
  86.             rtcol = [];
  87.             rtrow = [];
  88.             for j=1:blk(i,2)
  89.               rtrow = sbs(rtrow,sisorat(omegap,row(j)));
  90.             end
  91.             for j=1:blk(i,1)
  92.               rtcol = abv(rtcol,sisorat(omegap,col(j)));
  93.             end
  94.             block = mmult(rtcol,rtrow);
  95.             sys = daug(sys,block);
  96.           end
  97.           loc = loc + blk(i,1)*blk(i,2);
  98.         end
  99.       end
  100.     end
  101.   end
  102.   rat_pert = sys;
  103. %
  104. % Copyright MUSYN INC 1991,  All Rights Reserved
  105.