% Edwin Huffstutler, EE451, 9/17/94
% ---------------------------------------------------------
% Program to place arbitrary z-plane poles and zeros
% Input is in polar coordinates
% Conjugate poles/zeros automatically
% Symmetric poles/zeros (about unit circle) also created automatically
% ----------------------------------------------------------

% First draw a unit circle
    figure(1);clg;
    [x,y] = pol2cart(-pi:2*pi/360:pi,1);
    plot(x,y,'g:');
    axis([-2 2 -2 2]); axis('image');
    hold on;
    line([0 0],[-1.5 1.5]);line([-1.5 1.5],[0 0]);

% Now place poles
    fprintf(1,'Placing z-plane poles\n:');
    n = 1; clear p;
    while(1==1), 
        theta = input('Pole angle (degrees): ');
        if (length(theta)==0) break; end;
        r = input('Pole radius: ');
        [x,y] = pol2cart(pi*theta/180,r);
        [x2,y2] = pol2cart(pi*theta/180,1/r);
        p(n) = x+i*y; n=n+1;
        plot(x,y,'wx')
        if (r~=1)
           plot(x2,y2,'wx');
           p(n) = x2+i*y2; n=n+1;
        end;
        if ((theta~=0) & (theta~=180) )
            plot(x,-y,'wx');
            p(n) = x-i*y; n=n+1;
            if (r~=1)
                plot(x2,-y2,'wx');
                p(n) = x2-i*y2; n=n+1;
            end;
        end;
    end;                

% Now place zeros
    fprintf(1,'Placing z-plane zeros\n:');
    n = 1; clear z;
    while(1==1)
        theta = input('Zero angle (degrees): ');
        if (length(theta)==0) break; end;
        r = input('Zero radius: ');
        [x,y] = pol2cart(pi*theta/180,r);
        [x2,y2] = pol2cart(pi*theta/180,1/r);
        z(n) = x+i*y; n=n+1;
        plot(x,y,'wo');
        if (r~=1)
           plot(x2,y2,'wo');
           z(n) = x2+i*y2; n=n+1;
        end;
        if ((theta~=0) & (theta~=180) )
            plot(x,-y,'wo');
            z(n) = x-i*y; n=n+1;
            if (r~=1)
                plot(x2,-y2,'wo');
                z(n) = x2-i*y2; n=n+1;
            end;
        end;
    end;               

% Now get system function H(z)
    num = 1; den = 1;
    for n = 1:length(z)
        num = conv(num,[1 -z(n)])
    end;
    for n = 1:length(p)
        den = conv(den,[1 -p(n)])
    end;

% Now plot frequency response
    figure(2);
    h = zeros(201,1);
    [h,w] = freqz(num,den,200,'whole');
    mag=20*log10(abs(h));
    phase=(180/pi)*angle(h);
    mag(101:200)=mag(1:100);
    mag(1:100)=flipud(mag(1:100));
    phase(101:200)=phase(1:100);
    phase(1:100)=fliplr(phase(1:100));
    subplot(211);
    plot(w/pi-1,mag), title('Magnitude (dB)'); 
    subplot(212);
    plot(w/pi-1,phase), title('Phase (degrees)');
    xlabel('Normalized Frequency');