function [B,A] = f_low2bss (b,a,F_0,F_1)

%F_LOW2
%-------------------------------------------------------------------------------
 % Usage:       [B,A] = f_low2bss (b,a,F_0,F_1)
%
% Description: Convert normalized analog lowpass filter H(s) = b(s)/a(s) 
%              with cutoff frequency 1/(2*pi) to an analog bandstop filter 
%              with cutoff frequencies F_0 and F_1.
%
% Inputs:      b   = vector of length m+1 containing coefficients of analog
%                    numerator polynomial.
%              a   = vector of length n+1 containing coefficients of analog
%                    denominator polynomial (n >= m).
%              F_0 = lower stopband cutoff frequency in Hz
%              F_1 = upper stopband cutoff frequency in Hz
%
% Outputs:     B = (m+1) by 1 vector containing coefficients of highpass 
%                  numerator polynomial.
%              A = (n+1) by 1 vector containing coefficients of highpass
%                  denominator polynomial.
%
% See also:    f_low2lows, f_low2highs, f_low2bps
%-------------------------------------------------------------------------------

% Initialize

F_0 = f_clip (F_0,0,F_0);
F_1 = f_clip (F_1,0,F_1);
Omega_0 = 2*pi*F_0; 
Omega_1 = 2*pi*F_1; 
Omega = Omega_0*Omega_1;
Delta = Omega_1 - Omega_0;
m = length(b) - 1;
n = length(a) - 1;

% Compute poles and zeros

z = roots(b);
r1 = (Delta./z + sqrt(Delta^2./(z.^2) - 4*Omega))/2;
r2 = (Delta./z - sqrt(Delta^2./(z.^2) - 4*Omega))/2;
r = [f_torow(r1) f_torow(r2)];
if n > m
   r(2*m+1:2*m+n-m) = j*sqrt(Omega);
   r(m+n+1:m+n+n-m) = -j*sqrt(Omega);
end
p = roots(a);
q1 = (Delta./p + sqrt(Delta^2./(p.^2) - 4*Omega))/2;
q2 = (Delta./p - sqrt(Delta^2./(p.^2) - 4*Omega))/2;
q = [f_torow(q1) f_torow(q2)];

% Compute gain

b0 = b(1);
beta = b0*((-1)^(m-n))*prod(z)/prod(p);

% Compute coefficients

B = beta*poly(r);
A = poly(q);

%---------------------------------------------------------------------------------