fs = 1e6;
fpass = 0.10*fs;
fstop = 0.13*fs;
rpass = 0.1;
rstop = 60;

f = linspace(0, fs/2, 1000);
w = 2*pi*f;

% Elliptic IIR
[Ne, We] = ellipord(fpass*2/fs, fstop*2/fs, rpass, rstop);
[be,ae] = ellip(Ne, rpass, rstop, We);
He = tf(be, ae, 1/fs);

% Remez FIR
dev = [(10^(rpass/20)-1)/(10^(rpass/20)+1)  10^(-rstop/20)]; 
[Nr,fo,ao,Wr] = remezord([fpass fstop], [1 0],dev, fs);
br = remez(Nr, fo, ao, Wr);
Hr = tf(br, 1, 1/fs);
Hr = Hr / 10^(rpass/40);

% output
Ne
Nr

Me = freqresp(He, w);
Me = Me(1,:);

Mr = freqresp(Hr, w);
Mr = Mr(1,:);

plot(f, 20*log10(abs(Mr)), 'r');
hold on;
plot(f, 20*log10(abs(Me)), 'b');
plot([0 fpass fpass], [-rpass -rpass -rstop], 'g:');
plot([0 fstop fstop fs/2], [0 0 -rstop -rstop], 'g--');
fixfig;
xlabel('Frequency  0...f_s/2  [Hz]');
ylabel('Magnitude  [dB]');
legend('FIR', 'IIR');
axis([min(f) max(f) -rstop-10 5]);
hold off;

% phase
plot(f, unwrap(angle(Mr))*180/pi);
xlabel('Frequency  [Hz]');
ylabel('Phase  [degrees]');
fixfig
return

[Gd, Fd] = grpdelay(br, 1, 1000, fs);
plot(Fd, Gd);


