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