bodeplotgui.m

在Matlab中可以使用的波特图的例子源码

      semilogx(w(peakinds),mag_a(peakinds),'o','Color',handles.exactColor);
   end

   if handles.FirstPlot,
      %If this is the first time, let Matlab do autoscaling, but save
      %the y-axis limits so that they will be unchanged as the plot changes.
      ylims=get(gca,'YLim')/20;
      ylims(1)=min(-20,ceil(ylims(1))*20);
      ylims(2)=max(20,floor(ylims(2))*20);
      handles.MagLims=ylims;
   else
      %If this is not the first time, retrieve the old y-axis limits.
      ylims=handles.MagLims;
   end
   set(gca,'YLim',ylims);
   set(gca,'YTick',ylims(1):20:ylims(2));  %Make ticks every 20 dB.

   ylabel('Magnitude - dB');
   xlabel('Frequency - \omega, rad-sec^{-1}')
   title('Magnitude Plot','color','b','FontWeight','bold');
   if get(handles.GridCheckBox,'Value')==1,
      grid on
   end
   hold off;
   %%%%%%%%%%%%%%%%%%%% End Magnitude Plot %%%%%%%%%%


   %%%%%%%%%%%%%%%%%%%% Start Phase Plot %%%%%%%%%%
   %Much of this section mirrors the previous section and is not commented.
   %One difference is that phase is calculated explicitly, rather than use
   %Matlab's "bode" command.  Since phase is not unique (you can add or
   %subtract multiples of 360 degrees) There were sometimes discrepancies
   %between Matlab's phase calculations and mine
   axes(handles.PhasePlot);
   cla;

   %Plot line at 0 degrees for reference.
   semilogx([MinFreq MaxFreq],[0 0],...
      'Color',handles.zrefColor,...
      'LineWidth',1.5);
   hold on;

   %The variable phs_a has the combined asymptotic phase.
   phs_a=zeros(size(w));
   for i=1:length(Terms),
      if Terms(i).display,
         switch Terms(i).type,
            case 'Constant',
               f=[MinFreq MaxFreq];
               if Terms(i).value>0,
                  p=[0 0];
               else
                  p=[180 180];
               end
               if get(handles.RadianCheckBox,'Value')==1,
                  p=p/180;
               end
               semilogx(f,p,...
                  'LineWidth',handles.LnWdth,...
                  'LineStyle',GetLineStyle(i,handles),...
                  'Color',GetLineColor(i,handles));
            case 'RealPole',
               wo=-Terms(i).value;
               f=[MinFreq wo/10 wo*10 MaxFreq];
               p=[0 0 -90 -90]*Terms(i).multiplicity;
               if get(handles.RadianCheckBox,'Value')==1,
                  p=p/180;
               end
               semilogx(f,p,...
                  'LineWidth',handles.LnWdth,...
                  'LineStyle',GetLineStyle(i,handles),...
                  'Color',GetLineColor(i,handles));
            case 'RealZero',
               if Terms(i).value>0,  %Non-minimum phase
                  wo=Terms(i).value;
                  %Uncomment the next section to get agreement with Matlab plots.
                  %if Terms(1).value>0,  %Choose 0 or 360 to agree with MatLab plots
                  %    p=[0 0 0 0];      % (based on sign of constant term).
                  % else
                  %    p=[360 360 360 360];
                  %end
                  p=[0 0 -90 -90]*Terms(i).multiplicity;
               else
                  %Minimum phase
                  wo=-Terms(i).value;
                  p=[0 0 90 90]*Terms(i).multiplicity;
               end
               f=[MinFreq wo/10 wo*10 MaxFreq];
               if get(handles.RadianCheckBox,'Value')==1,
                  p=p/180;
               end
               semilogx(f,p,...
                  'LineWidth',handles.LnWdth,...
                  'LineStyle',GetLineStyle(i,handles),...
                  'Color',GetLineColor(i,handles));
            case 'ComplexPole',
               wo=abs(Terms(i).value);
               %good for zeta>0.05;
               bf=min(0.5*log10(2/zeta),0.99);  %For small zeta, bf=1
               % bf to 0.99 so f (next line) is four distinct values
               % so that interpolation works properly (it needs distinct
               % values).
               f=[MinFreq wo*bf wo/bf MaxFreq];
               p=[0 0 -180 -180]*Terms(i).multiplicity;
               if get(handles.RadianCheckBox,'Value')==1,
                  p=p/180;
               end
               semilogx(f,p,...
                  'LineWidth',handles.LnWdth,...
                  'LineStyle',GetLineStyle(i,handles),...
                  'Color',GetLineColor(i,handles));
            case 'ComplexZero',
               wo=abs(Terms(i).value);
               bf=min(0.5*log10(2/zeta),0.99);
               f=[MinFreq wo*bf wo/bf MaxFreq];
               if real(Terms(i).value)>0,   %Non-minimum phase
                  p=[0 0 -180 -180]*Terms(i).multiplicity;
               else
                  p=[0 0 180 180]*Terms(i).multiplicity;
               end
               if get(handles.RadianCheckBox,'Value')==1,
                  p=p/180;
               end
               semilogx(f,p,...
                  'LineWidth',handles.LnWdth,...
                  'LineStyle',GetLineStyle(i,handles),...
                  'Color',GetLineColor(i,handles));
            case 'OriginPole',
               f=[MinFreq MaxFreq];
               p=[-90 -90]*Terms(i).multiplicity;
               if get(handles.RadianCheckBox,'Value')==1,
                  p=p/180;
               end
               semilogx(f,p,...
                  'LineWidth',handles.LnWdth,...
                  'LineStyle',GetLineStyle(i,handles),...
                  'Color',GetLineColor(i,handles));

            case 'OriginZero',
               f=[MinFreq MaxFreq];
               p=[90 90]*Terms(i).multiplicity;
               if get(handles.RadianCheckBox,'Value')==1,
                  p=p/180;
               end
               semilogx(f,p,...
                  'LineWidth',handles.LnWdth,...
                  'LineStyle',GetLineStyle(i,handles),...
                  'Color',GetLineColor(i,handles));
         end
         phs_a=phs_a+interp1(log(f),p,log(w));  %Build up asymptotic approx.
      end
   end

   if get(handles.RadianCheckBox,'Value')==1,
      ph=ph/180;
   end

   semilogx(w,ph(:),...
      'Color',handles.exactColor,...
      'LineWidth',handles.LnWdth/2);  %Plot it.
   if (get(handles.ShowAsymptoticCheckBox,'Value')~=0),
      % There can be discrepancies between Matlabs calculation of phase and
      % the quantity "phs_a."  This is because phase is not unique (you can
      % add or subtract multiples of 360 degrees).  The next line shifts the
      % value of phs_a so that phs_a(1)=ph(1).  Ensuring they are the same
      % at the beginning of the plot ensures that they align elsewhere.  Note
      % that if the plots are already aligned (ph(1)=phs_a(1)), that the next
      % line does nothing.
      phs_a=phs_a+(ph(1)-phs_a(1));
      semilogx(w,phs_a,...
         'Color',handles.exactColor,...
         'LineStyle',':',...
         'LineWidth',handles.LnWdth);  %Plot asymptotic approx.
   end

   set(gca,'XLim',[MinFreq MaxFreq]);
   if handles.FirstPlot,
      ylims=get(gca,'YLim')/45;
      ylims(1)=ceil(ylims(1)-1)*45;
      ylims(2)=floor(ylims(2)+1)*45;
      handles.DPhaseLims=ylims;     %Find phase limits in degrees
      handles.RPhaseLims=ylims/180; %Find phase limits in radians/pi
   else
      if get(handles.RadianCheckBox,'Value')==1,
         ylims=handles.RPhaseLims;
      else
         ylims=handles.DPhaseLims;
      end
   end

   if get(handles.RadianCheckBox,'Value')==1,
      set(gca,'YLim',ylims);
      set(gca,'YTick',ylims(1):0.25:ylims(2))
      ylabel('Phase - radians/\pi');
   else
      set(gca,'YLim',ylims);
      set(gca,'YTick',ylims(1):45:ylims(2))
      ylabel('Phase - degrees');
   end

   xlabel('Frequency - \omega, rad-sec^{-1}');
   title('Phase Plot','color','b','FontWeight','bold');
   if get(handles.GridCheckBox,'Value')==1,
      grid on
   end
   hold off;
   %%%%%%%%%%%%%%%%%%%% End Phase Plot %%%%%%%%%%

   BodePlotDispTF(handles);    %Display the transfer function
   BodePlotLegend(handles);    %Display the legend.
   handles=guidata(handles.AsymBodePlot);  %Reload handles after function call.

   %Set first plot to zero (so Matlab won't autoscale on subsequent calls).
   handles.FirstPlot=0;

   guidata(handles.AsymBodePlot, handles);  %save changes to handles.
end
% ----------------- End of function BodePlotter ----------------------


% --------------------------------------------------------------------
function BodePlotSys(handles)
   %This function makes up a transfer function of all the terms that are not in
   %the "Excluded Elements" box in the GUI.

   Terms=handles.Terms;    %Get all terms from original transfer function.
   p=[];                   %Start with no poles, or zeros, and a constant of 1
   z=[];
   k=1;

   for i=1:length(Terms),      %For each term.
      %If the term is not in "Excluded Elements". then we want to display it.
      if Terms(i).display,
         switch Terms(i).type,
            case 'Constant',
               k=Terms(i).value;               %This is the constant.
            case 'RealPole',
               for j=1:Terms(i).multiplicity,
                  p=[p Terms(i).value];       %Add poles.
               end
            case 'RealZero',
               for j=1:Terms(i).multiplicity,
                  z=[z Terms(i).value];       %Add zeros.
               end
            case 'ComplexPole',
               for j=1:Terms(i).multiplicity,
                  p=[p Terms(i).value];       %Add poles.
                  p=[p conj(Terms(i).value)];
               end
            case 'ComplexZero',
               for j=1:Terms(i).multiplicity,
                  z=[z Terms(i).value];       %Add zeros.
                  z=[z conj(Terms(i).value)];
               end
            case 'OriginPole',
               for j=1:Terms(i).multiplicity,
                  p=[p 0];                    %Add poles.
               end
            case 'OriginZero',
               for j=1:Terms(i).multiplicity,
                  z=[z 0];                    %Add zeros.
               end
         end
      end
   end

   %Determine multiplicative constant in standard Bode Plot form.
   for i=1:length(p),
      if p(i)~=0
         k=-k*p(i);
      end
   end
   for i=1:length(z),
      if z(i)~=0
         k=-k/z(i);
      end
   end
   %If poles and/or zeros were complex conjugate pairs, there may be
   %some residual imaginary part due to finite precision.  Remove it.
   k=real(k);

   handles.sysInc=zpk(z,p,k);
   guidata(handles.AsymBodePlot, handles);  %save changes to handles.
end
% ------------------End of function BodePlotSys ----------------------


% --------------------------------------------------------------------
function BodePlotDispTF(handles)
   % This function displays a tranfer function that is a helper function for the
   % BodePlotGui routine.  It takes the transfer function of the numerator and
   % splits it into three lines so that it can be displayed nicely.  For example:
   % "            s + 1"
   % "H(s) = ---------------"
   % "        s^2 + 2 s + 1"
   %
   % The numerator string is in the variable nStr,
   % the second line is in divStr,
   % and the denominator string is in dStr.

   % Get numerator and denominator.
   [n,d]=tfdata(handles.sysInc,'v');