fp = .1; Ap = .25;  % User's parameters for
fs = .2; As =  40;  % discrete-time lowpass filter

wp = 2*tan(pi*fp);  % Prewarp DTFT frequencies
ws = 2*tan(pi*fs);

dp = 10^(Ap/10)-1;  % Compute order n
ds = 10^(As/10)-1;
n = ceil(log(ds/dp)/log(ws/wp)/2);

twon = 2*n; xx = 1/twon;      % For omega_c, use the average
wc = (wp/dp^xx+ws/ds^xx)/2;   % of its upper and lower bounds

k = [0:twon-1];
svec = wc*exp(j*pi*(n+2*k+1)/twon); % poles of H(s)H(-s)
lambda = svec(1:n);                 % poles of H(s)

t = linspace(0,2*pi,200);     % plot circle
x = wc*cos(t); y = wc*sin(t);
figure(1); plot(x,y); grid on;
hold on
plot(svec,'x')                % plot all 2n roots
plot(lambda,'o');             % circle LHP roots
hold off

lambdahat = (1+lambda/2)./(1-lambda/2); % poles of z transform
a = real(poly(lambdahat));    % difference eq. coefficients
b = poly(-ones(1,n));
b = (sum(a)/sum(b))*b;        % normalize to make H(z)=1 for z=1

f = linspace(0,1/2,2000);     % Plot H(exp(j*2*pi*f)) on dB scale
zinv = exp(-j*2*pi*f);
Hf = polyval(fliplr(b),zinv)./polyval(fliplr(a),zinv);
ff = [fp fs]; Hff = -[Ap As]; % Put o's on graph at design points
figure(2)
plot(f,10*log10(Hf.*conj(Hf)),ff,Hff,'ro'); grid on