dfilter.m

离散控制系统设计的MATLAB 代码

% Interactive plot of frequency response for a linear fixed
% discrete-time system
%
% To use program, enter dfilter in MATLAB command window and 
% response to prompts.

%%%%%%%%%%%%%%%%%%% dfilter.m %%%%%%%%%%%%%%%%%%%
%   Discrete-Time Control Problems using        %
%       MATLAB and the Control System Toolbox   %
%   by J.H. Chow, D.K. Frederick, & N.W. Chbat  %
%         Brooks/Cole Publishing Company        %
%                September 2002                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc
clear
%
% some definitions
%
  stpaus = ['                         ' ...
           '  Hit any key to continue'];
  sttitl = 'Enter title > ';
% prepare unit circle
%
  c = [0:pi/180:2*pi];
  c = c';
  cxy = [cos(c) sin(c)];
  cxy2 = cxy(1:181,:);
%
% create graphic window
%
  fign = figure;
%  set(fign,'Position',[46 440 560 420])
  set(fign,'Position',[46 100 500 420])
  axis('equal')
  axis([-1.2 1.2 -1.2 1.2])
  axis(axis)
  grid
  hold;
%
% unit circle
%
  plot(cxy(:,1),cxy(:,2),'m')
%
% input poles and zeroes
%
%  nfrq2

clc
%
% data statements
%
  streal = '  Real component of ';
  stimag = '  Imaginary component of ';
  stzero = '#';
  stpole = '#';
  stdef0 = ' (default = 0)>  ';
  stdef1 = ' (default = 1)>  ';
  stpaus = '  Hit any key to continue';
  sttyin = [ ...
             '      1) prompt' 
             '      2) mouse '
             '      3) file  '];
%
  zn = 1;
  pn = 1;
%
% get type of input for discrete frequency response
%
  disp(' ')
  disp('  Enter the number corresponding to the type of input:')
  disp(' ')
  disp(sttyin)
  disp(' ')
  disp(' ')
%echo on
%
% Enter the number corresponding to the type of input
% 
%	1) prompt
%	2) mouse
%	3) file
%
%echo off
%
  type1 = input('  Type of input (default = 1)> ');
%
  ansz = size(type1);
  if (ansz(1,1) == 0)
    type1 = 1;
  end
  type1 = fix(type1);
  if ((type1 < 1) | (type1 > 3))
    error('Choice out of bounds')
  end
  type = type1;
%
  clc
%
% grab input
%
% number of zeroes
%
  zn1 = input(['  Number of zeroes' stdef1]);
  if (zn1 > 0.5)
    zn = fix(zn1);
  end
%
%   set all initial zeroes to the origin
%
  zarr = zeros(zn,2);
%
% if file input, open file
%
  if (type == 3)
    fid = fopen('zero.dat','r');
  end
%
% initialize flag
%
  parflg = 0;
%
% get the zeroes
%
  for n = 1:zn
%
%   if conjugate pair
%
    if (parflg == 1)
      plot(zr1,-zi1,'go')
      zarr(n,1) = zr1;
      zarr(n,2) = -zi1;
      parflg = 0;
    else
%
%     keyboard input
%
      if (type ==1)
%
        stnumb = int2str(n);
        zr1 = input([streal stzero stnumb stdef0]);
        ansz = size(zr1);
        if (ansz(1,1) == 0)
          zr1 = 0;
        end
        zi1 = input([stimag stpole stnumb stdef0]);
        ansz = size(zi1);
        if (ansz(1,1) == 0)
          zi1 = 0;
        end
%
      end
%
%     mouse input
%
      if (type == 2)
%
        z11 = ginput(1);
        zr1 = z11(1,1);
        zi1 = z11(1,2);
        if (zi1 < 0.05)
          zi1 = 0;
        end
%
      end
%
%     file input
%
      if (type == 3)
%
        z11 = fscanf(fid,'%f',2);
        zr1 = z11(1,1);
        zi1 = z11(2,1);
%
      end
%
%     place zero on axis
%
      if ((type == 1 & zi1 < 0.04) | n == zn)
        zi1 = 0;
      end
%
      dist = sqrt(zi1*zi1+zr1*zr1);
      if ((type == 1 & zi1 < 0.04) | n == zn)
        zi1 = 0;
      end
      plot(zr1,zi1,'go')
      zarr(n,1) = zr1;
      zarr(n,2) = zi1;
      if (zi1 ~= 0)
        parflg = 1;
      end
%
    end
%
  end
%
% if file input, close file
%
  if (type == 3)
    fclose(fid);
  end
%
% display zeroes
%
  zarr
  stpaus
  pause
  clc
%
% number of poles
%
  pn1 = input(['  Number of poles' stdef1]);
  if (pn1 > 0.5)
    pn = fix(pn1);
  end
%
% if file input, open file
%
  if (type == 3)
    fid = fopen('pole.dat','r');
  end
%  
% set all initial poles to the origin
%
  parr = zeros(pn,2);
%
% initialize flag
%
  parflg = 0;
%
% get the poles
%
  for n = 1:pn
%
    if (parflg == 1)
      plot(pr1,-pi1,'r*')
      parr(n,1) = pr1;
      parr(n,2) = -pi1;
      parflg = 0;
    else
%
      if (type == 1)
%
%       keyboard input
%
        stnumb = int2str(n);
        pr1 = input([streal stzero stnumb stdef0]);
        ansz = size(pr1);
        if (ansz(1,1) == 0)
          pr1 = 0;
        end
        pi1 = input([stimag stpole stnumb stdef0]);
        ansz = size(pi1);
        if (ansz(1,1) == 0)
          pi1 = 0;
        end
%
      end
%
%     mouse input
%
      if(type == 2)
%
        p11 = ginput(1);
        pr1 = p11(1,1);
        pi1 = p11(1,2);
%
      end
%
      if (type == 3)
%
        p11 = fscanf(fid,'%f',2);
        pr1 = p11(1,1);
        pi1 = p11(2,1);
%
      end
%
      dist = sqrt(pi1*pi1+pr1*pr1);
      if (dist > 1)
        error('Position out of unit circle')
      end
      if ((type == 1 & pi1 < 0.04) | n == pn)
        pi1 = 0;
      end
      plot(pr1,pi1,'r*')
      parr(n,1) = pr1;
      parr(n,2) = pi1;
      if (pi1 ~= 0)
        parflg = 1;
      end
%
    end
%
  end
%
% if file input, close ifle
%
  if (type == 3)
    fclose(fid);
  end
%
% display the poles
%
  parr
  stpaus
  pause
%
% determine sampling time
%
  clc
  disp(' ')
  timsmp = input('  Enter sampling time > ');
  timsmp2 = timsmp*2;
  tmaxis = [0:1/(timsmp2*180):1/timsmp2];
end
%%%%%%%%%%


%
% close the display window
%
  close
%
% set variable form for later
%
  zrl = zarr(:,1);
  zim = zarr(:,2);
  prl = parr(:,1);
  pim = parr(:,2);
%
% loop on poles and zeroes for frequency response
%
% zeroes
%
  for k = (1:zn)
    zd(:,k) = (cxy2(:,1)-zrl(k)).*(cxy2(:,1)-zrl(k));
    zd(:,k) = zd(:,k)+(cxy2(:,2)-zim(k)).*(cxy2(:,2)-zim(k));
    zd(:,k) = sqrt(zd(:,k));
    za(:,k) = atan((cxy2(:,2)-zim(k))./(cxy2(:,1)-zrl(k)));
    zind = find(cxy2(:,1) < zrl(k));
    za(zind) = za(zind)+pi;
  end
%
% poles
%
  for k = (1:pn)
    pd(:,k) = (cxy2(:,1)-prl(k)).*(cxy2(:,1)-prl(k));
    pd(:,k) = pd(:,k)+(cxy2(:,2)-pim(k)).*(cxy2(:,2)-pim(k));
    pd(:,k) = sqrt(pd(:,k));
    pa(:,k) = atan((cxy2(:,2)-pim(k))./(cxy2(:,1)-prl(k)));
    pind = find(cxy2(:,1) < prl(k));
    pa(pind) = pa(pind)+pi;
  end
%
% initialize fields
  zdist = 1;
  pdist = 1;
  zang = 0;
  pang = 0;
%
% generate magnitude and phase from vectors of individual zeroes
% and poles
%
  if ( zn ~= 1)
    zdist = prod(zd');
    zang = sum(za');
  else
    zdist = zd';
    zang = za';
  end
  if ( pn ~= 1)
    pdist = prod(pd');
    pang = sum(pa');
  else
    pdist = pd';
    pang = pa';
  end
%
% prepare output for display
%
  mag = 20*log10(zdist./pdist);
  ang = (zang-pang)*180/pi;
%
% max and min of each vector for plot
%
  magmax = max(mag);
  mag = mag - magmax;   % JHC 2/22/97
  magmax = 0;           % max gain is 0 dB
  magmin = min(mag);
  mag = mag - magmax;
  angmax = max(ang);
  angmin = min(ang);
%
% graphic window 
%
%  wind

  fign = figure;
  set(fign,'Position',[46 50 500 600])
  
%
% set unit circle graphic
%
  a1 = axes('Position',[0.10 0.1 0.45 0.33]);
  axis('equal')
  axis([-1.2 1.2 -1.2 1.2])
  grid
  hold;
%
% set phase angle graphic
%
  a2 = axes('Position',[0.1 0.5  0.8 0.17]);
  xlabel('Frequency (Hz)')
  ylabel('Phase Angle (deg)')
  axis([0 1/(timsmp2) angmin angmax])
  grid
  hold;
%
% set magnitude graphic
%
  a3 = axes('Position',[0.1 0.75 0.8 0.17]);
  xlabel('Frequency (Hz)')
  ylabel('Mag(dB)')
  axis([0 1/(timsmp2) magmin magmax])
  grid
  hold;
%
% set pole zero text window
%
  a4 = axes('Position',[0.60 0.1 0.3 0.33]);
  axis('ij')
  axis([1 15 1 20])

%
% unit circle, poles, and zeroes
%
  set(fign,'CurrentAxes',a1)
  plot(cxy(:,1),cxy(:,2),'m')
  plot(zrl,zim,'go')
  plot(prl,pim,'r*')
%
% listing of poles and zeros
%
  set(fign,'CurrentAxes',a4)
  set(a4,'Visible','off')
  text(8,1,'Poles')
  for n = 1:pn
    if (parr(n,1) < 0) 
      text(2.3,2+n,num2str(parr(n,1)))
    else
      text(3,2+n,num2str(parr(n,1)))
    end
    if (parr(n,2) < 0) 
      text(9.3,2+n,num2str(parr(n,2)))
    else
      text(10,2+n,num2str(parr(n,2)))
    end
  end
  text(8,4+pn,'Zeros')
  for n = 1:zn
    if (zarr(n,1) < 0)
      text(2.3,5+pn+n,num2str(zarr(n,1)))
    else
      text(3,5+pn+n,num2str(zarr(n,1)))
    end
    if (zarr(n,2) < 0)
      text(9.3,5+pn+n,num2str(zarr(n,2)))
    else
      text(10,5+pn+n,num2str(zarr(n,2)))
    end
  end
  text(3.5,7+pn+zn,'(Real)')
  text(10.5,7+pn+zn,'(Imag)')
  sam_mes = ['Sampling period = ',num2str(timsmp),' sec.'];
  text(3.5,10+pn+zn,sam_mes)
  disp(stpaus)
  pause
%
% run movie
%
%  trcmvi

% loop for movie type display
  for n = (1:60)
%
%     unit circle trace
%
    set(fign,'CurrentAxes',a1)
    for k = (1:zn)
      lnm(1,:) = zarr(k,:);
      lnm(2,:) = cxy(n*3,:);
      if (n == 1)
        lzn(k) = line(lnm(:,1),lnm(:,2),'erasemode','xor');
      else
        set(lzn(k),'xdata',lnm(:,1),'ydata',lnm(:,2));
      end
    end
    for k = (1:pn)
      pnm(1,:) = parr(k,:);
      pnm(2,:) = cxy(n*3,:);
      if (n == 1)
        lpn(k) = line(pnm(:,1),pnm(:,2),'erasemode','xor');
      else
        set(lpn(k),'xdata',pnm(:,1),'ydata',pnm(:,2));
      end
    end
%
%     magnitude and phase response to trace
%
    if (n == 1)
      set(fign,'CurrentAxes',a3)
      p1 = plot(tmaxis(1:n*3),mag(1:n*3),'r-','erasemode','xor');
      set(fign,'CurrentAxes',a2)
      p2 = plot(tmaxis(1:n*3),ang(1:n*3),'r-','erasemode','xor');
    else
%     xtmp = [1:n*3];
      set(fign,'CurrentAxes',a3)
      set(p1,'xdata',tmaxis(1:n*3),'ydata',mag(1:n*3))
      set(fign,'CurrentAxes',a2)
      set(p2,'xdata',tmaxis(1:n*3),'ydata',ang(1:n*3))
    end
    pause(0.1)
%
%   clear when done
%
    if (n == 60)
      for k = (1:zn)
        delete(lzn(k))
      end
      for k = (1:pn)
        delete(lpn(k))
      end
    end
%
  end
%
% check for repeat movie
%
  clc
  disp(' ')
  strrpt = input('Repeat movie display (y or n)> ','s');
  if (strrpt == 'y')
    delete(p1)
    delete(p2)
%    trcmvi

% loop for movie type display
  for n = (1:60)
%
%     unit circle trace
%
    set(fign,'CurrentAxes',a1)
    for k = (1:zn)
      lnm(1,:) = zarr(k,:);
      lnm(2,:) = cxy(n*3,:);
      if (n == 1)
        lzn(k) = line(lnm(:,1),lnm(:,2),'erasemode','xor');
      else
        set(lzn(k),'xdata',lnm(:,1),'ydata',lnm(:,2));
      end
    end
    for k = (1:pn)
      pnm(1,:) = parr(k,:);
      pnm(2,:) = cxy(n*3,:);
      if (n == 1)
        lpn(k) = line(pnm(:,1),pnm(:,2),'erasemode','xor');
      else
        set(lpn(k),'xdata',pnm(:,1),'ydata',pnm(:,2));
      end
    end
%
%     magnitude and phase response to trace
%
    if (n == 1)
      set(fign,'CurrentAxes',a3)
      p1 = plot(tmaxis(1:n*3),mag(1:n*3),'r-','erasemode','xor');
      set(fign,'CurrentAxes',a2)
      p2 = plot(tmaxis(1:n*3),ang(1:n*3),'r-','erasemode','xor');
    else
%     xtmp = [1:n*3];
      set(fign,'CurrentAxes',a3)
      set(p1,'xdata',tmaxis(1:n*3),'ydata',mag(1:n*3))
      set(fign,'CurrentAxes',a2)
      set(p2,'xdata',tmaxis(1:n*3),'ydata',ang(1:n*3))
    end
    pause(0.1)
%
%   clear when done
%
    if (n == 60)
      for k = (1:zn)
        delete(lzn(k))
      end
      for k = (1:pn)
        delete(lpn(k))
      end
    end
%
  end
%
% check for repeat movie
%
  clc
  disp(' ')
  strrpt = input('Repeat movie display (y or n)> ','s');
  if (strrpt == 'y')
    delete(p1)
    delete(p2)
    trcmvi
  end
  
  end

%
% title
%
  set(fign,'CurrentAxes',a3)
  titplt = input(sttitl,'s')
  title(titplt) 
%
% print window to a file
%
  clc
  orient tall
  echo on
%
% to save your frequency response plot to a file, type
%
%     print -dps xxxxxx
%     (where xxxxxx) is a filename)
%
% to print the frequency response immediately, type
%
%     print
%
echo off
%%%%%%%%%%