% This script for discrete-time filter design displays the frequency response for filter
%     \mathbf{H}(e^{j\Omega}) = \frac{ C (e^{j\Omega}-\mathbf{z}_1)(e^{j\Omega}-\mathbf{z}_2)\,\cdot\cdot\cdot\,(e^{j\Omega}-\mathbf{z}_M)}{(e^{j\Omega}-\mathbf{p}_1)(e^{j\Omega}-\mathbf{p}_2)\,\cdot\cdot\cdot\,(e^{j\Omega}-\mathbf{p}_N)}
% Specify the gain constant C and up to six poles \textbf{p}_n and/or zeros \textbf{z}_m via radius r and angle \theta.  Any \theta other than 0 or \pm 1 \pi automatically creates a complex conjugate pair.

% LS10_1_PoleZero_z.mlx       % pole zero plot and spectrum for discrete filter
%   with 1 to 6 poles and/or zeros
format compact; close all;
pz = ["--","--","--","--","--","--","--"]; % intitalize parameter arrays 
p = []; z = [];                             
r = ones(1,6); theta = zeros(1,6);         

C = 1; % filter H gain
pz(1) = "ze"; r(1) =1; theta(1) =1;
pz(2) = "ze"; r(2) =1; theta(2) =0.75;
pz(3) = "ze"; r(3) =1; theta(3) =0.5;
pz(4) = "po"; r(4) =0.8; theta(4) =0;
pz(5) = "po"; r(5) =0.8; theta(5) =0.25;
pz(6) = "po"; r(6) =0.8; theta(6) =0.5;

for i = 1:6              % append values to pole, zero arrays
  if (theta(i) == 0)     % positive real axis, no conjugates
      vals = 1;
  elseif (theta(i) == 1) % negative real axis, no conjugates
      vals = -1;
  else                   % complex conjugate pair of poles or zeros
      vals = [exp(j*pi*theta(i)), exp(-j*pi*theta(i))];
  end

  if pz(i) == "po"
    p = [p, r(i)*vals];
  elseif pz(i) == "ze"
    z = [z, r(i)*vals];
  end
end

if isempty(p) & isempty(z)     
  disp('**** At least one pole or zero must be chosen via dropdown. ****');
else
  NUM = C * poly(z); DEN = poly(p);
  % spectrum plot of H(e^{jw})
  w = linspace(0,pi,150);       % frequencies for spectrum plot    
  H = polyval(NUM,exp(j*w)) ./ polyval(DEN,exp(j*w));  % H(e^{jw}) = transfer function values

  figure(position = [0,0,800,300]), 
   subplot(1,2,1), zplane(z',p')  
   subplot(1,2,2), plot(w,abs(H),'r', 'LineWidth', 1.2), grid minor, title('|H(e^{j\Omega})|'), 
      xlabel('\Omega (rad/sample)'), xticks(round([0, pi/4,pi/2,3*pi/4, pi],2)) , axis tight
end