bodeplotgui.m

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

   i=2;
   while (i<=length(term)),
      %If root is at origin, handle it separately
      if (term(i).value==0),
         term(i).type=['Origin' term(i).type];
         %If imaginary part is sufficiently small...
      elseif (abs(imag(term(i).value)/term(i).value)<Acc),
         term(i).type=['Real' term(i).type];  %...Add "Real" to type
         term(i).value=real(term(i).value);   %...And set imaginary part=0
         %If imaginary part is *not* small...
      else
         term(i).type=['Complex' term(i).type];  %...Add "Complex" to type
         term(i+1).multiplicity=0;               %...Remove complex conjugate
         i=i+1;                                  %...And skip conjugate.
      end
      i=i+1;  %Go to next root.
   end

   %Remove all terms with multiplicity 0.
   j=0;
   for i=1:length(term),
      if (term(i).multiplicity~=0)
         j=j+1;
         term(j)=term(i);
         term(j).display=1;
      end
   end

   term=term(1:j);

   DoQuit=0;
   %Check for poles or zeros in right half plane,
   % or on imaginary axis. Poles and zeros at origin are OK.
   if any(real(p)>0),  %Poles in RHP.
      beep;
      waitfor(errordlg('System has poles with positive real part, cannot make plot.'));
      DoQuit=1;
      return;
   end
   if any(real(z)>0),  %Zeros in RHP.
      disp(' ');
      disp(' ');
      disp(' ');
      disp('************Warning******************');
      handles.Sys
      disp('is a nonminimum phase system (zeros in right half plane).');
      disp('The standard rules for phase do not apply.');
      disp(' ');
      disp('Also - The plots produced may be different than the Matlab Bode plot');
      disp('  by a factor of 360 degrees.  So though the graphs don''t look');
      disp('  the same, they are equivalent');
      disp(' ');
      disp('Location(s) of zero(s):');
      disp(z);
      disp('*************************************');
      disp(' ');
      beep;
      waitfor(warndlg('System has zeros with positive real part.  See Command Window for caveats.'));
   end
   %Check for terms near imaginary axis, or multiple poles or zeros at origin.
   for i=2:length(term),
      if (term(i).value~=0),
         if (abs(real(term(i).value)/term(i).value)<Acc),
            disp(' ');
            disp(' ');
            disp(' ');
            disp('************Warning******************');
            handles.Sys
            disp('has a pole or zero near, or on, the imaginary axis.');
            disp('The plots may be inaccurate near that frequency.')
            disp(' ');
            disp('--------');
            disp('Pole(s):');
            disp(p)
            disp('--------');
            disp('Zero(s):');
            disp(z);
            disp('*************************************');
            disp(' ');
            disp(' ');
            disp(' ');
            beep;
            waitfor(warndlg('System has poles or zeros with real part near zero.  See Command Window.'));
         end
      elseif (term(i).multiplicity>1),
         disp(' ');
         disp(' ');
         disp(' ');
         disp('************Warning******************');
         handles.Sys
         disp('has multiple poles or zeros at the origin.');
         disp('Components of the phase plot may appear to disagree.');
         disp('This is because the phase of a complex number is');
         disp('not unique; the phase of -1 could be +180 degrees');
         disp('or -180 or +/-540...  Likewise the phase of 1/s^2');
         disp('could be +180 degrees, or -180 degrees (or +/-540');
         disp('Keep this in mind when looking at the phase plots.');
         disp('*************************************');
         disp(' ');
         beep;
         waitfor(warndlg('System has multiple poles or zeros at origin.  See Command Window.'));
      end
   end

   handles.Terms=term;
   guidata(handles.AsymBodePlot, handles);  %save changes to handles.
end
% ------------------End of function BodePlotTerms --------------------


% --------------------------------------------------------------------
function BodePlotter(handles)

   %Get the constituent terms and the system itself.
   Terms=handles.Terms;
   sys=handles.Sys;

   %Call function to get a system with only included poles and zeros.
   BodePlotSys(handles);
   handles=guidata(handles.AsymBodePlot);  %Reload handles after function call.

   %Get systems with included poles and zeros.
   sysInc=handles.sysInc;

   %Find min and max freq.
   MinF=realmax;
   MaxF=-realmax;
   for i=2:length(Terms),
      if Terms(i).value~=0,
         MinF=min(MinF,abs(Terms(i).value));
         MaxF=max(MaxF,abs(Terms(i).value));
      end
   end

   %If there is exclusively a pole or zero at origin, MinF and MaxF will
   % not have changed.  So set them arbitrarily to be near unity.
   if MaxF==-realmax,
      MinF=0.9;
      MaxF=1.1;
   end

   %MinFreq is a bit more than an order of magnitude below lowest break.
   MinFreq=10^(floor(log10(MinF)-1.01));
   %MaxFreq is a bit more than an order of magnitude above highest break.
   MaxFreq=10^(ceil(log10(MaxF)+1.01));
   %Calculate 500 frequency points for plotting.
   w=logspace(log10(MinFreq),log10(MaxFreq),1000);

   %%%%%%%%%%%%%%%%%%%% Start Magnitude Plot %%%%%%%%%%
   axes(handles.MagPlot);
   cla;

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

   hold on;

   %For each term, plot the magnitude accordingly.
   %The variable mag_a has the combined asymptotic magnitude.
   mag_a=zeros(size(w));
   %The variable peakinds holds the indices at which peaks in underdamped
   %responses occur.
   peakinds=[];
   for i=1:length(Terms),
      if Terms(i).display,
         switch Terms(i).type,
            case 'Constant',
               %A constant term is unchanging from beginning to end.
               f=[MinFreq MaxFreq];
               m=20*log10(abs([Terms(i).value Terms(i).value]));
               semilogx(f,m,...
                  'LineWidth',handles.LnWdth,...
                  'LineStyle',GetLineStyle(i,handles),...
                  'Color',GetLineColor(i,handles));
               mag_a=mag_a+interp1(log(f),m,log(w));  %Build up asymptotic approx.
            case 'RealPole',
               %A real pole has a single break frequency and then
               %decreases at 20 dB per decade (Or more if pole is multiple).
               wo=-Terms(i).value;
               f=[MinFreq wo MaxFreq];
               m=-20*log10([1 1 MaxFreq/wo])*Terms(i).multiplicity;
               semilogx(f,m,...
                  'LineWidth',handles.LnWdth,...
                  'LineStyle',GetLineStyle(i,handles),...
                  'Color',GetLineColor(i,handles));
               mag_a=mag_a+interp1(log(f),m,log(w));  %Build up asymptotic approx.
            case 'RealZero',
               %Similar to real pole, but increases instead of decreasing.
               wo=abs(Terms(i).value);
               f=[MinFreq wo MaxFreq];
               m=20*log10([1 1 MaxFreq/wo])*Terms(i).multiplicity;
               semilogx(f,m,...
                  'LineWidth',handles.LnWdth,...
                  'LineStyle',GetLineStyle(i,handles),...
                  'Color',GetLineColor(i,handles));
               mag_a=mag_a+interp1(log(f),m,log(w));  %Build up asymptotic approx.
            case 'ComplexPole',
               %A complex pole has a single break frequency and then
               %decreases at 40 dB per decade (Or more if pole is multiple).
               %There is also a peak value whose height and location are
               %determined by the natural frequency and damping coefficient.
               %We will plot a circle ('o') at the location of the peak.
               wn=abs(Terms(i).value);
               theta=atan(abs(imag(Terms(i).value)/real(Terms(i).value)));
               zeta=cos(theta);
               if (zeta < (1/sqrt(2))),  %Peaking occurs if zeta<0.707
                  wo=wn*sqrt(1-2*zeta*zeta);
                  peak=2*zeta*sqrt(1-zeta*zeta);
                  f=[MinFreq wo wo wo wn MaxFreq];
                  m=-20*log10([1 1 peak 1 1 (MaxFreq/wn)^2])*Terms(i).multiplicity;
                  semilogx(f,m,...
                     'LineWidth',handles.LnWdth,...
                     'LineStyle',GetLineStyle(i,handles),...
                     'Color',GetLineColor(i,handles));
                  semilogx(wo,-20*log10(peak),'o','Color',GetLineColor(i,handles));
                  % Set up for interpolation (w/o peak)
                  f=[MinFreq wn MaxFreq];
                  m=-20*log10([1  1 (MaxFreq/wn)^2])*Terms(i).multiplicity;
                  mag_a=mag_a+interp1(log(f),m,log(w));  %Build up asymptotic approx.
                  %Find location closest to peak, and adjust its amplitude.
                  index=min(find(w>=wo));
                  mag_a(index)=mag_a(index)-20*log10(peak);
                  peakinds=[peakinds index];    %Save this index.
               else
                  f=[MinFreq wn MaxFreq];
                  m=-20*log10([1 1 (MaxFreq/wn)^2])*Terms(i).multiplicity;
                  semilogx(f,m,...
                     'LineWidth',handles.LnWdth,...
                     'LineStyle',GetLineStyle(i,handles),...
                     'Color',GetLineColor(i,handles));
                  % Set up for interpolation (w/o peak)
                  mag_a=mag_a+interp1(log(f),m,log(w));  %Build up asymptotic approx.
               end
            case 'ComplexZero',
               %Similar to complex pole, but increases instead of decreasing.
               wn=abs(Terms(i).value);
               theta=atan(abs(imag(Terms(i).value)/real(Terms(i).value)));
               zeta=cos(theta);
               if (zeta < (1/sqrt(2))),  %Peaking occurs if zeta<0.707
                  wo=wn*sqrt(1-2*zeta*zeta);
                  peak=2*zeta*sqrt(1-zeta*zeta);
                  f=[MinFreq wo wo wo wn MaxFreq];
                  m=20*log10([1 1 peak 1 1 (MaxFreq/wn)^2])*Terms(i).multiplicity;
                  semilogx(f,m,...
                     'LineWidth',handles.LnWdth,...
                     'LineStyle',GetLineStyle(i,handles),...
                     'Color',GetLineColor(i,handles));
                  semilogx(wo,20*log10(peak),'o','Color',GetLineColor(i,handles));
                  % Set up for interpolation (w/o peak)
                  f=[MinFreq wn MaxFreq];
                  m=20*log10([1  1 (MaxFreq/wn)^2])*Terms(i).multiplicity;
                  mag_a=mag_a+interp1(log(f),m,log(w));  %Build up asymptotic approx.
                  %Find location closest to peak, and adjust its amplitude.
                  index=min(find(w>=wo));
                  mag_a(index)=mag_a(index)+20*log10(peak);
                  peakinds=[peakinds index];    %Save this index.
               else
                  f=[MinFreq wn MaxFreq];
                  m=20*log10([1 1 (MaxFreq/wn)^2])*Terms(i).multiplicity;
                  semilogx(f,m,...
                     'LineWidth',handles.LnWdth,...
                     'LineStyle',GetLineStyle(i,handles),...
                     'Color',GetLineColor(i,handles));
                  % Set up for interpolation (w/o peak)
                  mag_a=mag_a+interp1(log(f),m,log(w));  %Build up asymptotic approx.
               end
            case 'OriginPole',
               %A pole at the origin is a monotonically decreasing straigh line.
               f=[MinFreq MaxFreq];
               m=-20*log10([MinFreq MaxFreq])*Terms(i).multiplicity;
               semilogx(f,m,...
                  'LineWidth',handles.LnWdth,...
                  'LineStyle',GetLineStyle(i,handles),...
                  'Color',GetLineColor(i,handles));
               mag_a=mag_a+interp1(log(f),m,log(w));  %Build up asymptotic approx.
            case 'OriginZero',
               %Similar to pole at origin, but increases instead of decreasing.
               f=[MinFreq MaxFreq];
               m=20*log10([MinFreq MaxFreq])*Terms(i).multiplicity;
               semilogx(f,m,...
                  'LineWidth',handles.LnWdth,...
                  'LineStyle',GetLineStyle(i,handles),...
                  'Color',GetLineColor(i,handles));
               mag_a=mag_a+interp1(log(f),m,log(w));  %Build up asymptotic approx.
         end
      end
   end

   %Set the x-axis limits to the minimum and maximum frequency.
   set(gca,'XLim',[MinFreq MaxFreq]);

   [mg,ph,w]=bode(sysInc,w);   %Calculate the exact bode plot.
   semilogx(w,20*log10(mg(:)),...
      'Color',handles.exactColor,...
      'LineWidth',handles.LnWdth/2);  %Plot it.
   if (get(handles.ShowAsymptoticCheckBox,'Value')~=0),
      semilogx(w,mag_a,...
         'Color',handles.exactColor,...
         'LineStyle',':',...
         'LineWidth',handles.LnWdth);  %Plot asymptotic approx.