scilab/modules/cacsd/macros/bode.sci

   1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   2 // Copyright (C)  1985-2010 - INRIA - Serge Steer
   3 // This file must be used under the terms of the CeCILL.
   4 // This source file is licensed as described in the file COPYING, which
   5 // you should have received as part of this distribution.  The terms
   6 // are also available at
   7 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
   8
   9 function []=bode(varargin)
  10   rhs=size(varargin)
  11   if type(varargin($))==10 then
  12     comments=varargin($),rhs=rhs-1;
  13   else
  14     comments=[];
  15   end
  16   fname="bode";//for error messages
  17   fmax=[]
  18   discr=%f //for shannon limit
  19   if or(typeof(varargin(1))==["state-space" "rational"]) then
  20     //sys,fmin,fmax [,pas] or sys,frq
  21     refdim=1 //for error message
  22     discr=varargin(1).dt<>'c';
  23     if rhs==1 then //sys
  24       [frq,repf]=repfreq(varargin(1),1d-3,1d3)
  25     elseif rhs==2 then //sys,frq
  26       if size(varargin(2),2)<2 then
  27         error(msprintf(_("%s: Wrong size for input argument #%d: A row vector with length>%d expected.\n"),..
  28                        fname,2,1))
  29       end
  30       [frq,repf]=repfreq(varargin(1:rhs))
  31     elseif or(rhs==(3:4)) then //sys,fmin,fmax [,pas]
  32       [frq,repf]=repfreq(varargin(1:rhs))
  33     else
  34       error(msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"),fname,1,5))
  35     end
  36     [phi,d]=phasemag(repf)
  37     if rhs>=3 then fmax=varargin(3),end
  38   elseif  type(varargin(1))==1 then
  39     //frq,db,phi [,comments] or frq, repf [,comments]
  40     refdim=2
  41     select rhs
  42     case 2 then //frq,repf
  43       frq=varargin(1);
  44       if size(frq,2)<2 then
  45         error(msprintf(_("%s: Wrong size for input argument #%d: A row vector with length>%d expected.\n"),..
  46                        fname,1,1))
  47       end
  48       if size(frq,2)<>size(varargin(2),2) then
  49         error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"),..
  50                        fname,1,2))
  51       end
  52       [phi,d]=phasemag(varargin(2))
  53     case 3 then  //frq,db,phi
  54       [frq,d,phi]=varargin(1:rhs)
  55       if size(frq,2)<>size(d,2) then
  56         error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"),..
  57                        fname,1,2))
  58       end
  59       if size(frq,2)<>size(phi,2) then
  60         error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"),..
  61                        fname,1,3))
  62       end
  63     else
  64       error(msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"),fname,2,4))
  65     end
  66   else
  67     error(msprintf(_("%s: Wrong type for input argument #%d: Linear dynamical system or row vector of floats expected.\n"),fname,1))
  68   end;
  69   frq=frq';d=d',phi=phi'
  70   [n,mn]=size(d)
  71
  72   if comments==[] then
  73     hx=0.48
  74   else
  75     if size(comments,"*")<>mn then
  76       error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same number of elements expected.\n"),...
  77                      fname,refdim,rhs+1))
  78     end
  79     hx=0.43
  80   end;
  81
  82   drawlater()
  83   sciCurAxes=gca();
  84   axes=sciCurAxes;
  85   wrect=axes.axes_bounds;
  86
  87
  88   //magnitude
  89   axes.axes_bounds=[wrect(1)+0,wrect(2)+0,wrect(3)*1.0,wrect(4)*hx*0.95]
  90   axes.data_bounds = [min(frq),min(d);max(frq),max(d)];
  91   axes.log_flags = "lnn" ;
  92   axes.grid=color("lightgrey")*ones(1,3);
  93   axes.axes_visible="on";
  94   axes.clip_state = "clipgrf";
  95   if size(d,2)>1&size(frq,2)==1 then
  96     xpolys(frq(:,ones(1,mn)),d,1:mn)
  97   else
  98     xpolys(frq,d,1:mn)
  99   end
 100   //set datatips info
 101   e=gce();
 102
 103   for i=1:size(e.children,"*")
 104     datatipInitStruct(e.children(i),"formatfunction","formatBodeMagTip")
 105   end
 106
 107   if discr&fmax<>[]&max(frq)<fmax then
 108     xpoly(max(frq)*[1;1],axes.y_ticks.locations([1 $]));e=gce();
 109     e.foreground=5;
 110   end
 111   xtitle("",_("Frequency (Hz)"),_("Magnitude (dB)"));
 112
 113   //phase
 114
 115   axes=newaxes();
 116   axes.axes_bounds=[wrect(1)+0,wrect(2)+wrect(4)*hx,wrect(3)*1.0,wrect(4)*hx*0.95];
 117   axes.data_bounds = [min(frq),min(phi);max(frq),max(phi)];
 118   axes.log_flags = "lnn" ;
 119   axes.grid=color("lightgrey")*ones(1,3);
 120   axes.axes_visible="on";
 121   axes.clip_state = "clipgrf";
 122   if size(phi,2)>1&size(frq,2)==1 then
 123     xpolys(frq(:,ones(1,mn)),phi,1:mn)
 124   else
 125     xpolys(frq,phi,1:mn)
 126   end
 127   ephi=gce()
 128   //set datatips info
 129   for i=1:size(ephi.children,"*")
 130     datatipInitStruct(ephi.children(i),"formatfunction","formatBodePhaseTip")
 131   end
 132
 133   if discr&fmax<>[]&max(frq)<fmax then
 134     xpoly(max(frq)*[1;1],axes.y_ticks.locations([1 $]));e=gce();
 135     e.foreground=5;
 136   end
 137   xtitle("",_("Frequency (Hz)"),_("Phase (degree)"));
 138   // create legend
 139   if comments<>[] then
 140     c=captions(ephi.children,comments,"lower_caption")
 141     c.background=get(gcf(),"background")
 142   end
 143   drawnow()
 144   // return to the previous scale
 145   set( "current_axes", sciCurAxes ) ;
 146
 147 endfunction
 148 function str=formatBodeMagTip(curve,pt,index)
 149 //this function is called by the datatips mechanism to format the tip
 150 //string for the magnitude bode curves
 151   str=msprintf("%.4g"+_("Hz")+"\n%.4g"+_("dB"), pt(1),pt(2))
 152 endfunction
 153 function str=formatBodePhaseTip(curve,pt,index)
 154 //this function is called by the datatip mechanism to format the tip
 155 //string for the bode phase curves
 156   str=msprintf("%.4g"+_("Hz")+"\n %.4g"+_("dB"), pt(1),pt(2))
 157 endfunction