% USAGE:
% cplane (SYS, RES, AX, FUNC, RESP)
% Plot the magnitude or phase of a transfer function in a complex plane
% Inputs:
% SYS
% LTI model, must be a single input and single output model
% RES
% (optional) Real and imaginary resolution of complex plane
% AX
% (optional) Axis limits of real and imaginary axes.
% FUNC
% (optional) define mesh plot function. Options are: 'mesh', 'meshc','meshz', 'surf', 'surfc', 'pcolor'
% RESP
% (optional) specify type of response, either "mag" (default) or "phase".
% Output:
% [oxx,oyy,ozz]
% (optional) The coordinatates to call plot functions, or export to other tools
%
% See also pzmap, bode, bodemag
%
% examples:
% % 4th order butterworth low-pass filter, corner frequency=1000 [rad/s]
% [num, den]=butter(4,1000,'s');
% h=tf(num,den);
% cplane(h);
%
% %lowpass elliptic/Cauer filter, 5th-order, 0.1/20dB's pass/stop ripple, fc=20 [rad/s]
% [num, den]=ellip(5,0.1,20,20,'s');
% h=tf(num,den);
% cplane(h,150,[-20 0 -30 30]);
%
% %bessel 9th order low-pass filter, fc=100 [rad/s]
% [num, den]=besself(9,100,'s');
% h=tf(num,den);
% cplane(h,[],[], 'surfc'); % use the contour surface function
% view(0,90); shading interp; colormap gray; %create a log density plot view
%
% Copyright (C) 2012 Ronald van der Meer
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2,
% or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; If not, see <http://www.gnu.org/licenses/>.
function [oxx, oyy, ozz] = cplane(sys,res,ax,func,resp)
%narginchk(1, 5) %check the correct number of arguments
if (nargin < 5) || (exist('resp','var')==0) || (length(resp)==0);
resp='mag'; %default response type
end
if (nargin < 4) || (exist('func','var')==0) || (length(func)==0);
func='mesh'; %default plot function
end
if (nargin < 3) || (exist('ax','var')==0) || (length(ax)==0); %determine if auto-scaling is needed
pz=[pole(sys) ;zero(sys)];
rmin=min(real(pz));
rmax=max(real(pz));
imin=min(imag(pz));
imax=max(imag(pz));
a=0.2; %This is the scaling factor e.g. 0.2 = 20% beyond the min/max of the pole/zero's
ax=[rmin-(rmax-rmin)*a rmax+(rmax-rmin)*a imin-(imax-imin)*a imax+(imax-imin)*a];
end
if (nargin < 2) || (exist('res','var')==0) || (length(res)==0);
res=80; %the default xy resolution of the plot
end;
x=linspace(ax(1),ax(2),res);
y=linspace(ax(3),ax(4),res);
[xx,yy]=meshgrid(x,y);
z=polyval(sys.num{1},xx+i*yy)./polyval(sys.den{1},xx+i*yy);
if (strcmp(resp,'mag'))
zz=20*log10(abs(z)); %convert to [dB] magnitude
else
zz=atan2(imag(z),real(z)); %do phase instead of magnitide
end
if (nargout==0) %when there are no output arguments
eval(strcat(func,'(x,y,zz)'));
axis('tight');
set(gca,'TickDir','out') %let the tick-lines be outside the plot. This will position the tick-text a little better
xlabel('Real Axis');
ylabel('Imaginary Axis');
if (strcmp(resp,'mag'))
zlabel('Magnitude [dB]');
else
zlabel('Phase [rad]');
end
%change the view so we look "into" the negetive real half plane
[az,el]=view;
az=az+90;
view(az,el);
else %don't plot but return coordinates
oxx=x;
oyy=y;
ozz=zz;
end
end
%!demo
%! [num, den]=ellip(5,0.1,20,20,'s');
%! h=tf(num,den);
%! cplane(h);
Here is some example octave code that produces a 4th order Butterworth filter complex plane plot:
graphics_toolkit fltk
figure
[num, den]=butter(4,1000,'s');
h=tf(num,den);
cplane(h);
print -dpng butter4example.png -FArialMT:8 "-S800,600"
system('convert butter4example.png -trim butter4example.png')
Here is some example octave code that produces an elliptic 5-th order filter complex plane plot:
graphics_toolkit fltk
figure
[num, den]=ellip(5,0.1,20,10,'s');
h=tf(num,den);
cplane(h,150,[],'meshc');
print -dpng cauer5example.png -FArialMT:8 "-S1200,1000"
system('convert cauer5example.png -trim cauer5example.png')
Here is some example octave code for a pole-zero map, with matching log density plot:
graphics_toolkit fltk
figure
[num, den]=ellip(5,0.1,20,10,'s');
h=tf(num,den);
subplot(1,2,1)
pzmap(h);
legend('poles','zeros');
axislimits=[get(gca,'xlim') , get(gca,'ylim')]
axis('equal');
subplot(1,2,2);
cplane(h,80,axislimits);
set(gca,'TickDir','in') %let the tick-lines be inside to match pzmap
title('Log density of h')
axis('equal');
view(0,90); shading interp; colormap gray;
print -dpng logdensityexample.png -FArialMT:8 "-S1000,1100"
system('convert logdensityexample.png -trim logdensityexample.png')
Here is some example octave code for a notch filter with matching frequency magnitude and phase plots:
%graphics_toolkit fltk
figure(1)
wc=10; %center freq.
bw=10; %-3dB bandwidth=10, Q-factor=1
num = [1 0 wc^2]; %notch numerator polynomials
den = [1 bw wc^2]; %notch denominator polynomials
h=tf(num,den);
cplane(h,[],[-6 0 0 16 ]);
print -dpng notchexample_tmp1.png -FArialMT:6 "-S540,350"
system('convert notchexample_tmp1.png -trim notchexample_cmag.png')
figure(2)
cplane(h,[],[-6 0 0 16 ],[],'phase')
print -dpng notchexample_tmp2.png -FArialMT:6 "-S540,350"
system('convert notchexample_tmp2.png -trim notchexample_cphase.png')
w=linspace(0,16,100);
[mag,phase,w]=bode(h,w);
figure(3)
plot(w,20*log(mag))
xlabel('freq in [rad/s]');
ylabel('Magnitude [dB]');
grid;
axis tight;
print -dpng notchexample_tmp3.png -FArialMT:6 "-S490,200"
system('convert notchexample_tmp3.png -trim notchexample_mag.png')
figure(4)
plot(w,phase/180*pi)
xlabel('freq in [rad/s]');
ylabel('Phase [rad]');
grid;
axis tight;
print -dpng notchexample_tmp4.png -FArialMT:6 "-S490,200"
system('convert notchexample_tmp4.png -trim notchexample_phase.png')
|
|
|
|
|
|