Archive
Archive for April, 2012
plot_z_surface
April 7, 2012
1 comment
This function plots the z surface given the positions of poles and zeros. Also used by the zpgui function
% plot_z_surface.
% By David Dorran (david.dorran@dit.ie)
%
%
function plot_z_surface(pole_positions_orig, zero_positions_orig, varargin)
surface_limit = 1.5;
min_val = -surface_limit;
max_val = surface_limit;
CameraPos=[3.7719 -15.7111 275.3541];
CameraUpVec=[-0.1237 0.5153 23.1286];
if(nargin >2 )
CameraPos= varargin{1};
CameraUpVec=varargin{2};
end
surface_display_opts = 0;
if(nargin == 5)
surface_display_opts = varargin{3};
end
%check if any of the poles or zeros fall outside the maximum values
%specified
if length(find(abs(real([pole_positions_orig zero_positions_orig])) > max_val)) || length(find(abs(imag([pole_positions_orig zero_positions_orig])) > max_val))
error(['the real part of the complex numbers indicating the the positions of the poles and zeros must be less than ' num2str(max_val) ' and greater than -' num2str(max_val) '. Also, the imaginary part of the complex numbers indicating the the positions of the poles and zeros must be less than ' num2str(max_val) 'j and greater than -' num2str(max_val) 'j.']);
end
%set up the s_plane grid from -5j to 5j along the imaginary axis and -5 to
%5 along the real axis. A resolution of 0.02 provides a visually pleasing
%plot
surface_resolution = 0.02;
a_vals = min_val:surface_resolution:max_val;
b_vals = a_vals;
grid_len = length(a_vals);
%create a grid of values along the z surface indicating position on
%the z-domain
ones_vector = ones(1, grid_len);
z_grid = ones_vector'*a_vals + (fliplr(b_vals*j))'*ones_vector;
%round all the poles and zeros to 1 decimal place because of resolution limitation
%Also, add .01+0.01j to each pole and zero position so that they are positioned in the
%center of a grid square of the z-surface created above. This ensures that
%the all the poles and zeros appear the same on the z-surface given the
%resolution constraints
pole_positions = round(pole_positions_orig*10)/10 + surface_resolution/2+(surface_resolution/2)*j;
zero_positions = round(zero_positions_orig*10)/10 +surface_resolution/2 +(surface_resolution/2)*j;
z_surface = zeros(length(a_vals)); % initialise the s_surface
% Since zeros are numerator terms obtain an log s_plane surface for each zero and add each individual plane
% since a logarithmic addition is equivalent to a multiply i.e.
% log(a.b) = log(a)+log(b)
for k = 1 : length(zero_positions)
z_surface = z_surface + 20*log10(abs(z_grid - zero_positions(k)));
% transfer_function_numerator = [transfer_function_numerator '(s-' num2str(zero_positions_orig(k)),')'];
end
% Since poles are denominator terms obtain an log s_plane surface for each pole and subtract each individual plane
% since a logarithmic subtraction is equivalent to a divide i.e.
% log(a/b) = log(a)-log(b)
transfer_function_denominator = [];
for k = 1 :length(pole_positions)
z_surface = z_surface - 20*log10(abs(z_grid - pole_positions(k) ));
transfer_function_denominator = [transfer_function_denominator '(z-' num2str(pole_positions_orig(k)),')'];
end
%
% for k = 1 :length(pole_positions)
% round(imag(pole_positions(k))/surface_resolution)
% z_surface(round(-imag(pole_positions(k))/surface_resolution)+surface_limit/surface_resolution, round(real(pole_positions(k))/surface_resolution)+surface_limit/surface_resolution) = 40;
% end
% for k = 1 :length(zero_positions)
% z_surface(round(imag(zero_positions(k))/surface_resolution)+surface_limit/surface_resolution, round(real(zero_positions(k))/surface_resolution)+surface_limit/surface_resolution) = -40;
% end
%half_z_surface = s_surface;
%half_s_surface(:,1:round(length(a_vals)/2)) = 20*log10(0.001);
% transfer_function = num2str(gain);
% if(length(transfer_function_numerator))
% if(gain ~= 1)
% transfer_function = [num2str(gain) '(' transfer_function_numerator ')'];
% else
% transfer_function = transfer_function_numerator;
% end
% end
% if(length(transfer_function_denominator))
% transfer_function = [ transfer_function '/(' transfer_function_denominator ')'];
% end
total_indices = 1:numel(z_grid);
temp_zero = 20*log10(abs(z_grid+.000001));
indices_outside_unit_circle = find(temp_zero > -0.2);
indices_inside_unit_circle = find(temp_zero < 0.2);
indices_on_unit_circle = intersect(indices_outside_unit_circle, indices_inside_unit_circle);
indices_not_on_unit_circle = setdiff( total_indices, indices_on_unit_circle );
positive_imaginary_indices = find(imag(z_grid) > 0 );
negative_inaginary_indices = find(imag(z_grid) < 0 );
positive_indices_outside_unit_circle = intersect(positive_imaginary_indices, indices_outside_unit_circle);
mask = ones(size(z_surface));
if(surface_display_opts == 1)
mask_indices = indices_not_on_unit_circle;
else
mask_indices = [];
end
mask(mask_indices) = NaN;
colors = z_surface;
colors(indices_on_unit_circle) = 40;
mesh(a_vals,b_vals,z_surface.*mask, colors)
set(gca,'YTick',[min_val:max_val])
set(gca,'YTickLabel',[min_val:max_val]')
set(gca,'XTick',[min_val:max_val])
set(gca,'XTickLabel',[min_val:max_val]')
set(gca,'CameraPosition',CameraPos);
set(gca,'CameraUpVector',CameraUpVec);
xlh = ylabel('Im');
ylh =xlabel('Re');
zlh = zlabel('Mag(db)');
axis([-1.1 1.1 -1.1 1.1 -30 30])
axis square
box on
Categories: matlab code
zpgui
April 7, 2012
15 comments
The zpgui function was developed by Tom Krauss in Perdue – I’ve adapted it to allow the z surface be displayed. It relies on another function I’ve written called plot_z_surface
Here is a link to updated files from Johan Peltenburg which work with the latest version of matlab (as at 2016) download zip file here
function varargout = zpgui(varargin)
% ZPGUI Zero Pole dragging Graphic User Interface
% Allows you to add, remove, and move the zeros and
% poles of a filter interactively.
% To begin, type:
% zpgui
% Author: Tom Krauss, 9/1/98
% Adapted by David Dorran (2011) to force poles and zeros to be conjugate
% pairs (maybe not the most flexible feature but it saves me time). Updated in 2012 to plot the z-surface.
global zh ph Lresp Nfft b a z h Y
global ax1 ax2 ax3 z_surface_CameraPos z_surface_CameraUpVec surface_display_opts
if nargin == 0
action = 'init';
z = [0.29-0.41j 0.29+0.41j]' ;
p = [-0.6+0.73j -0.6-0.73j]' ;
elseif nargin >= 2
p = varargin{1}';
z = varargin{2}';
action = 'init';
else
action = varargin{1};
end
if nargin == 3
jpg_filename = varargin{3};
end
switch action
case 'init'
z_surface_CameraPos=[3.7719 -15.7111 275.3541];
z_surface_CameraUpVec=[-0.1237 0.5153 23.1286];
surface_display_opts = 0;;
set(0,'defaultaxesfontsize',10)
subplot(2,2,1)
[zh,ph,cruff]=zplaneplot(z,p);
%set(cruff,'hittest','off')
ax1 = gca;
ylabel('Imaginary Part')
xlabel('Real Part')
uicontrol('style','pushbutton',...
'string','Add Zeros',...
'fontsize',8,...
'units','normalized',...
'position',[.15 .46 .18 .04],...
'callback','zpgui(''addzero'')');
uicontrol('style','pushbutton',...
'string','Remove Zeros',...
'fontsize',8,...
'units','normalized',...
'position',[.55 .46 .18 .04],...
'callback','zpgui(''removezero'')');
uicontrol('style','pushbutton',...
'units','normalized',...
'fontsize',8,...
'position',[.35 .46 .18 .04],'string','Add Poles',...
'callback','zpgui(''addpole'')');
uicontrol('style','pushbutton',...
'units','normalized',...
'fontsize',8,...
'position',[.75 .46 .18 .04],'string','Remove Poles',...
'callback','zpgui(''removepole'')');
uicontrol('style','pushbutton',...
'units','normalized',...
'fontsize',8,...
'position',[.5 .96 .18 .04],'string','Toggle surface display',...
'callback','zpgui(''toggle_surface_display'')');
subplot(2,1,2)
[b,a]=zp2tf(z,p,1);
Nfft = 512;
Y = fft(b,Nfft)./fft(a,Nfft);
Lresp = plot((0:Nfft-1)/Nfft*2-1, 20*log10(fftshift(abs(Y))),'linewidth',2);
ax2 = gca;
set(ax2,'xlim',[0 1])
%set(ax2,'xtick', [0:1/16:0.5])
%set(ax2,'xticklabel', '0 \pi/8 \pi/8 \pi/8 \pi/8 \pi/8 \pi/8 \pi/8 ')
%get(ax2,'xticklabel')
grid on
xlabel('Normalised Frequency (x\pi rads/sample)')
ylabel('Magnitude (db) ')
set(Lresp,'erasemode','xor')
set(zh,'buttondownfcn','zpgui(''zeroclick'')',...
'markersize',8,'linewidth',1)
set(ph,'buttondownfcn','zpgui(''poleclick'')',...
'markersize',8,'linewidth',1)
subplot(2,2,2)
plot_z_surface(p,z);
ax3 = gca;
set(gcf, 'ToolBar', 'figure');
if(nargin == 3) % create some figures ans save them to a file
figure
set(gcf,'position', [ 232 387 417 279]);
plot_z_surface(p,z);
saveas(gcf, strcat('z_surface_', jpg_filename, '.bmp'), 'bmp');
figure
set(gcf,'position', [ 232 387 417 279]);
plot((0:Nfft-1)/Nfft - .5, 20*log10(fftshift(abs(Y))),'linewidth',2);
set(gca,'xlim',[0 .5])
grid on
xlabel('Normalised Frequency')
ylabel('Magnitude (db) ')
saveas(gcf, strcat('z_response_', jpg_filename, '.bmp'), 'bmp');
figure
set(gcf,'position', [ 232 387 417 279]);
zplaneplot(z,p)
ylabel('Imaginary Part')
xlabel('Real Part')
saveas(gcf, strcat('z_pz_', jpg_filename, '.bmp'), 'bmp');
figure
set(gcf,'position', [232 605 392 61])
b_coeff_vals = poly(z)
a_coeff_vals = poly(p)
text(0.05,0.1 ,strcat('a = [',num2str(a_coeff_vals) ,']'))
text(0.05,0.8 ,strcat('b = [',num2str(b_coeff_vals) ,']'))
set(gca, 'visible','off')
pause(0.1)
saveas(gcf, strcat('z_ba_', jpg_filename, '.bmp'),'bmp');
close all
end
case 'toggle_surface_display'
surface_display_opts = surface_display_opts + 1;
if(surface_display_opts >= 2)
surface_display_opts = 0;
end
zpgui('recompute')
case 'addzero'
if length(zh)>0
zh(end+1) = line(.5,0,'parent',ax1,'buttondownfcn','zpgui(''zeroclick'')',...
'markersize',8,'linewidth',1,'marker','o','linestyle','none');
else
zh = line(.5,0,'parent',ax1,'buttondownfcn','zpgui(''zeroclick'')',...
'markersize',8,'linewidth',1,'marker','o','linestyle','none');
set(findobj('string','Remove Zeros'),'enable','on')
end
zpgui('recompute')
case 'removezero'
delete(zh(end))
zh(end)=[];
if length(zh)==0
set(gco,'enable','off')
end
zpgui('recompute')
case 'addpole'
if length(ph)>0
ph(end+1) = line(.5,0,'parent',ax1,'buttondownfcn','zpgui(''poleclick'')',...
'markersize',8,'linewidth',1,'marker','x','linestyle','none');
else
ph = line(.5,0,'parent',ax1,'buttondownfcn','zpgui(''poleclick'')',...
'markersize',8,'linewidth',1,'marker','x','linestyle','none');
set(findobj('string','Remove Poles'),'enable','on')
end
zpgui('recompute')
case 'removepole'
delete(ph(end))
ph(end)=[];
if length(ph)==0
set(gco,'enable','off')
end
zpgui('recompute')
case 'zeroclick'
set(gcf,'userdata','')
set(gcf,'windowbuttonmotionfcn','set(gcf,''userdata'',''motion'')')
set(gcf,'windowbuttonupfcn','set(gcf,''userdata'',''up'')')
ind = find(zh==gco);
%set(zh(ind),'erasemode','xor')
% set(Lresp,'erasemode','xor')
pair = length(get(zh(ind),'xdata'))==2;
done = 0;
pt = get(ax1,'currentpoint');
pt = pt(1,1:2);
title(['selected position: ' num2str(pt) 'j'])
while ~done
waitfor(gcf,'userdata')
switch get(gcf,'userdata')
case 'motion'
pt = get(ax1,'currentpoint');
pt = pt(1,1:2);
title(['selected position: ' num2str(pt) 'j'])
%if pair
if 1
set(zh(ind),'xdata',[pt(1) pt(1)],'ydata',[pt(2) -pt(2)])
else
set(zh(ind),'xdata',pt(1),'ydata',pt(2))
end
zpgui('recompute')
case 'up'
done = 1;
end
set(gcf,'userdata','')
end
set(gcf,'windowbuttonmotionfcn','')
set(gcf,'windowbuttonupfcn','')
set(zh(ind),'erasemode','normal')
set(Lresp,'erasemode','normal')
set(ax2,'ylimmode','auto')
ylim = get(ax2,'ylim');
Y = get(Lresp,'ydata');
if ylim(1)>min(Y)-3 | ylim(2)>max(Y)+3
set(ax2,'ylim',[min(Y)-3 max(Y)+3])
end
case 'poleclick'
set(gcf,'userdata','')
set(gcf,'windowbuttonmotionfcn','set(gcf,''userdata'',''motion'')')
set(gcf,'windowbuttonupfcn','set(gcf,''userdata'',''up'')')
ind = find(ph==gco);
% set(ph(ind),'erasemode','xor')
%set(Lresp,'erasemode','xor')
pair = length(get(ph(ind),'xdata'))==2;
done = 0;
pt = get(ax1,'currentpoint');
pt = pt(1,1:2);
title(['selected position: ' num2str(pt) 'j'])
while ~done
waitfor(gcf,'userdata')
switch get(gcf,'userdata')
case 'motion'
pt = get(ax1,'currentpoint');
pt = pt(1,1:2);
title(['selected position: ' num2str(pt) 'j'])
if 1
set(ph(ind),'xdata',[pt(1) pt(1)],'ydata',[pt(2) -pt(2)])
else
set(ph(ind),'xdata',pt(1),'ydata',pt(2))
end
zpgui('recompute')
case 'up'
done = 1;
end
set(gcf,'userdata','')
end
set(gcf,'windowbuttonmotionfcn','')
set(gcf,'windowbuttonupfcn','')
set(ph(ind),'erasemode','normal')
set(Lresp,'erasemode','normal')
set(ax2,'ylimmode','auto')
Y = get(Lresp,'ydata');
ylim = get(ax2,'ylim');
if ylim(1)>min(Y)-3 | ylim(2)<max(Y)+3
set(ax2,'ylim',[min(Y)-3 max(Y)+3])
end
case 'recompute'
z = [];
p = [];
b = 1;
a = 1;
for i=1:length(zh)
zx = get(zh(i),'xdata');
zy = get(zh(i),'ydata');
if length(zx)==1
b = conv(b,[1 -(zx+sqrt(-1)*zy)]);
else
b = conv(b,[1 -2*zx(1) zx(1).^2+zy(1).^2]);
end
z = [z zx+sqrt(-1)*zy];
end
for i=1:length(ph)
px = get(ph(i),'xdata');
py = get(ph(i),'ydata');
if length(px)==1
a = conv(a,[1 -(px+sqrt(-1)*py)]);
else
a = conv(a,[1 -2*px(1) px(1).^2+py(1).^2]);
end
p = [p px+sqrt(-1)*py];
end
Y = fft(b,Nfft)./fft(a,Nfft);
%Y = Y/max(abs(Y));
set(Lresp,'ydata',20*log10(fftshift(abs(Y))))
subplot(2,2,2)
z_surface_CameraPos = get(ax3, 'CameraPosition');
z_surface_CameraUpVec = get(ax3, 'CameraUpVector');
plot_z_surface(p,z, z_surface_CameraPos, z_surface_CameraUpVec, surface_display_opts);
end
Categories: matlab code
dadorran 