% bscheb2.m - bandstop Chebyshev type 2 digital filter design
%
% [A, B, P] = bscheb2(fs, fpa, fpb, fsa, fsb, Apass, Astop)
%
% design parameters:
% P = [Wpass, Wstop, Wsa, Wsb, c, epass, estop, Nex, N, f3a, f3b, a];
% A, B are Kx5 matrices of cascaded fourth-order sections

function [A, B, P] = bscheb2(fs, fpa, fpb, fsa, fsb, Apass, Astop)

c = sin(2*pi*(fpa + fpb)/fs)/(sin(2*pi*fpa/fs) + sin(2*pi*fpb/fs));
Wpass = abs(sin(2*pi*fpb/fs) / (cos(2*pi*fpb/fs) - c));
Wsa = sin(2*pi*fsa/fs) / (cos(2*pi*fsa/fs) - c);
Wsb = sin(2*pi*fsb/fs) / (cos(2*pi*fsb/fs) - c);
Wstop = min(abs(Wsa), abs(Wsb));

epass = sqrt(10^(Apass/10) - 1);
estop = sqrt(10^(Astop/10) - 1);

Nex = acosh(estop/epass) / acosh(Wstop/Wpass);
N = ceil(Nex);  r = rem(N,2);  K = (N - r) / 2;

a = asinh(estop) / N;
W3 = Wstop / cosh(acosh(estop)/N);
f3a = (fs/pi) * atan((sqrt(1 + W3^2 *(1 - c^2)) - 1)/(W3*(c+1)));
f3b = (fs/pi) * atan((sqrt(1 + W3^2 *(1 - c^2)) + 1)/(W3*(c+1)));
P = [Wpass, Wstop, Wsa, Wsb, c, epass, estop, Nex, N, f3a, f3b, a];
W0 = sinh(a) / Wstop;                                   % reciprocal of text

if r==1,
    G = 1 / (1 + W0);
    a1 = -2 * c / (1 + W0);
    a2 = (1 - W0) / (1 + W0);
    A(1,:) = [1, a1, a2, 0, 0];
    B(1,:) = G * [1, -2*c, 1, 0, 0];
else
    A(1,:) = [1, 0, 0, 0, 0];
    B(1,:) = [1, 0, 0, 0, 0];
end

for i=1:K,
    th = pi * (N - 1 + 2 * i) / (2 * N);
    Wi = sin(th) / Wstop;                               % reciprocal of text
    D = 1 - 2 * W0 * cos(th) + W0^2 + Wi^2;
    G = (1 + Wi^2) / D;
    b1 = -4 * c / (1 + Wi^2);
    b2 = 2 * ((2*c^2 + 1) - Wi^2) / (1 + Wi^2);
    a1 = -4 * c * (1 - W0 * cos(th)) / D;
    a2 = 2 * (2*c^2 + 1 - W0^2 - Wi^2) / D;
    a3 = -4 * c * (1 + W0 * cos(th)) / D;
    a4 = (1 + 2 * W0 * cos(th) + W0^2 + Wi^2) / D;
    A(i+1,:) = [1, a1, a2, a3, a4];
    B(i+1,:) = G * [1, b1, b2, b1, 1];
end