home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / BODEPLOT.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  4.1 KB  |  129 lines
  1. function bodeplot(G,sd,C,mode)
  2. %BODEPLOT  plots the Bode diagram of a transfer function or spectrum
  3. %
  4. %    bodeplot(G)  or  bodeplot(G,SD)
  5. %
  6. %    G is entered as the transfer function(s) or spectrum (spectra)
  7. %    of the standard frequency function format, produced by TH2FF,
  8. %    SPA or ETFE. See HELP FREQFUNC.
  9. %    
  10. %    Several plots in the same diagram are obtained by
  11. %    bodeplot([G1 G2 .. Gn]) 
  12. %    The transfer functions Gk need not be given at the same frequencies.
  13. %    If SD is specified to a number larger than zero, also confidence
  14. %    intervals for the functions (corresponding to SD standard deviations)
  15. %    are shown as dash-dotted  curves.
  16. %
  17. %    The default mode is that amplitude and phase plots are shown simul-
  18. %    taneously, for each input (noise spectrum) present in G. For spectra
  19. %    the phase plots are omitted. ENTER will advance the plot from one
  20. %    input to the next (if any). To show amplitude or phase plots only 
  21. %    use bodeplot(G,SD,'A') and bodeplot(G,SD,'P'), respectively. To ob-
  22. %    tain all plots in the same diagrams use bodeplot(G,SD,C,'same').
  23. %    C is here 'A', 'P' or 'B' for amplitude, phase or both plots.
  24. %    (In these cases SD can be omitted for default)
  25.  
  26. %    L. Ljung 10-1-86, revised 7-7-87,10-2-92
  27. %    Copyright (c) 1986-92 by the MathWorks, Inc.
  28. %    All Rights Reserved.
  29.  
  30. if nargin<4 , mode=[];end
  31. if nargin<3, C=[];end
  32. if nargin<2, sd=[];end
  33. if isempty(mode),mode='sep';end
  34. if isempty(C),C='B';end
  35. if isempty(sd),sd=0;end
  36.  
  37. abp=['a';'A';'b';'B';'p';'P'];
  38. if sum(sd(1)==abp)>0, if length(C)>1,mode=C;end,C=sd;sd=0;end
  39. if length(sd)>1,mode=sd;sd=0;end
  40.  
  41. mode=mode(1:3);
  42. if mode=='SEP', mode='sep';end
  43. if mode=='SAM', mode='sam';end
  44. if mode~='sam', mode='sep';end
  45.  
  46. if C == 'a', C='A'; end
  47. if C == 'p', C='P'; end
  48. if C == 'b', C='B'; end
  49.  
  50. [m,n]=size(G);
  51. inputs1=G(1,:);mo=fix(inputs1/1000);mmo=max(mo);
  52. kko=[]; % The output indices
  53. for k=[0:mmo],if length(find(mo==k))>0,kko=[kko k];end,end
  54. G=G(2:m,:);
  55. for ko=kko+1
  56.  
  57. inputs=inputs1-(ko-1)*1000;
  58. mi=max(inputs(find(inputs<19)));
  59.  
  60. kki=[]; % The input indices
  61. for k=[1:mi 0], if length(find(inputs==k))>0,kki=[kki k];end,end
  62. if sum(kki)==0 & C=='B', C='A';end
  63. if mode=='sep'
  64.   for k=kki
  65.   if k==0 & C=='B',C='A';end
  66.   clf reset,if C=='B',subplot(211),else subplot(111),end
  67.   if C=='A'| C=='B'
  68.  
  69.     loglog(G(:,find(inputs==100+k)),G(:,find(inputs==k)))
  70.     if sd>0,hold on,
  71.     ind=find(k+50==inputs);
  72.     loglog(G(:,ind-2),G(:,ind-1)+sd*G(:,ind),'-.')
  73.     loglog(G(:,ind-2),G(:,ind-1)-sd*G(:,ind),'-.') %Change here 2.04
  74.     hold off,end
  75.     xlabel('frequency (rad/sec)')
  76.  
  77.     if k~=0,title(['AMPLITUDE PLOT, input # ',int2str(k),' output # ',int2str(ko)])
  78.     else title(['SPECTRUM output # ',int2str(ko)]),end
  79.   end
  80.   if ((C=='B')|(C=='P'))&k~=0
  81.         if (C=='B'), subplot(212), end
  82.         if (C=='P'),clf,subplot(111),end
  83.  
  84.     semilogx(G(:,find(inputs==100+k)),G(:,find(inputs==20+k)))
  85.     if sd>0,hold on,
  86.     ind=find(k+70==inputs);
  87.     semilogx(G(:,ind-4),G(:,ind-1)+G(:,ind)*sd,'-.')
  88.     semilogx(G(:,ind-4),G(:,ind-1)-sd*G(:,ind),'-.') %Change here 2.04
  89.     hold off,end
  90.     title(['PHASE PLOT, input # ',int2str(k),' output # ',int2str(ko)])
  91.     xlabel('frequency (rad/sec)')
  92.     ylabel('phase')
  93.   end
  94.  
  95.  
  96.   if k~=kki(length(kki)),pause,end
  97. end
  98. end
  99. if mode=='sam' 
  100.  
  101.     if C=='B',subplot(211),else subplot(111),end
  102.       if C=='B' | C=='A'    
  103.     loglog(G(:,find(inputs>99&inputs<120)),G(:,find(inputs<20&inputs>-1)))
  104.     if sd>0,hold on,
  105.     ind=find(inputs<69 & inputs>49);
  106.     loglog(G(:,ind-2),G(:,ind-1)+G(:,ind)*sd,'-.')
  107.     loglog(G(:,ind-2),G(:,ind-1)-sd*G(:,ind),'-.'),hold off,end
  108.     xlabel('frequency (rad/sec)')
  109.     title(['AMPLITUDE PLOT output #',int2str(ko)])
  110.       end
  111.       if C=='B'| C=='P'
  112.          if C=='B',subplot(212),else subplot(111),end
  113.         semilogx(G(:,find(inputs>100&inputs<120)),G(:,find(inputs<39 & inputs>19)))
  114.     if sd>0,hold on
  115.     ind=find(inputs<89 & inputs>69);
  116.     semilogx(G(:,ind-4),G(:,ind-1)+G(:,ind)*sd,'-.')
  117.     semilogx(G(:,ind-4),G(:,ind-1)-sd*G(:,ind),'-.'),hold off,end
  118.     xlabel('frequency (rad/sec)')
  119.     ylabel('phase')
  120.     title(['PHASE PLOT output #',int2str(ko)])
  121.       end
  122.     
  123. end,if ko~=kko(length(kko))+1,pause,end
  124.  
  125. end
  126. hold off
  127. set(gca,'NextPlot','replace');
  128. set(gcf,'NextPlot','replace');
  129.