%% Copyright (C) 2021 Tony Richardson
%%
%% This program is free software: you can redistribute it and/or modify
%% it under the terms of the GNU General Public License as published by
%% the Free Software Foundation, either version 3 of the License, or
%% (at your option) any later version.
%%
%% This program is distributed in the hope that it will be useful,
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%% GNU General Public License for more details.
%%
%% You should have received a copy of the GNU General Public License
%% along with this program.  If not, see <https://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @deftypefn {} {@var{[B, A]} =} lp2bp (@var{Bp}, @var{Ap}, @var{w0}, @var{W})
%% Transform a lowpass s-plane filter into an s-plane bandpass filter.
%% Bp and Ap are numerator and denominator coefficients of the lowpass
%% prototype filter.  w0 and W are the center frequency and bandwidth
%% (rad/s) of the desired bandpass filter.  The numerator and denominator
%% coefficients of the desired bandpass filter are returned.
%%
%% @seealso{lp2lp, lp2hp, lp2bs}
%% @end deftypefn

%% Author: Tony Richardson <richardson.tony@gmail.com>
%% Created: 2021-04-09

function [B, A] = lp2bp(Bp, Ap, w0, W)
  N = length(Ap);
  M = length(Bp);
  NBP = 2*(N - 1);
  MBP = 2*(M - 1);
  v1 = [1 0  w0^2]; v2 = [ W 0 ]; 
  A = zeros(1, NBP + 1);
  for n = 1:N
    vT = 1;
    for k = 1:(N-n)
      vT = conv(vT, v1);      
    endfor
    for k = 1:(n-1)
      vT = conv(vT, v2);      
    endfor
    A = A + [zeros(1,NBP+1-length(vT)) Ap(n)*vT];
  endfor
  B = zeros(1, MBP + 1);
  for n = 1:M
    vT = 1;
    for k = 1:(M-n)
      vT = conv(vT, v1);      
    endfor
    for k = 1:(n-1)
      vT = conv(vT, v2);      
    endfor
    B = B + [zeros(1,MBP+1-length(vT)) Bp(n)*vT];
  endfor
  while (M < N)
    B = conv(B, v2);
    M = M + 1;
  end
  NORM = A(1);
  A = A/NORM;
  B = B/NORM;
end

%!demo
%! n = 14;
%! [z p k] = buttap(n);
%! [b a] = zp2tf(z, p, real(k));
%!
%! w = logspace(-1, 1, 501);
%! H = freqs(b, a, w);
%! figure(1)
%! clf
%! subplot(2, 1, 1)
%! loglog(w, abs(H))
%! ylim([1e-15 1]);
%! yticks(logspace(-15, 0, 4))
%! grid on
%! subplot(2, 1, 2)
%! semilogx(w, rad2deg(angle(H)))
%! grid on
%!
%! fl = 30; fh = 100;
%! Wo = 2*pi*sqrt(fl*fh);
%! Bw = 2*pi*(fh - fl);
%! [bt, at] = lp2bp(b, a, Wo, Bw);
%!
%! w = logspace(0, 4, 501);
%! H = freqs(bt, at, w);
%!
%! figure(2)
%! clf
%! subplot(2, 1, 1)
%! loglog(w, abs(H))
%! ylim([1e-20 1]);
%! yticks(logspace(-20, 0, 3))
%! grid on
%! subplot(2, 1, 2)
%! semilogx(w, rad2deg(angle(H)))
%! grid on
