%% Copyright (C) 1999 Paul Kienzle <pkienzle@users.sf.net>
%%
%% 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; see the file COPYING. If not, see
%% <https://www.gnu.org/licenses/>.

%% -*- texinfo -*-
%% @deftypefn  {Function File} {[@var{numd}, @var{dend}] =} bilinear2 (@var{num}, @var{den}, @var{fs})
%% @deftypefnx {Function File} {[@var{zd}, @var{pd}, @var{kd}] =} bilinear2 (@var{z}, @var{p}, @var{k}, @var{fs})
%% @deftypefnx {Function File} {[@var{Ad}, @var{Bd}, @var{Cd}, @var{Dd}] =} bilinear2 (@var{A}, @var{B}, @var{C}, ,@var{D}, @var{fs})
%% @deftypefnx {Function File} {[__] =} bilinear2 (__, @var{fp})
%%
%% Transform a s-plane filter specification into a z-plane specification. Filters 
%% can be specified in either transfer function, zero-pole-gain or
%% state-space form. fs is the sampling frequency. The
%% fp parameter does prewarping to match fp in the s and z planes
%%
%% blinear2 differs from the bilinear function in the signal processing
%% toolbox, which uses T as an input argument rather than fs (1/T).  bilinear2
%% is compatible with the MATLAB function.  The bilinear function does the
%% transformation on the poles and zeros while bilinear2 does the transformation
%% on the coefficients in the transfer function representation.
%%
%% Demos are included, try "demo bilinear2"
%%
%% @end deftypefn

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

function [Ad, Bd, Cd, Dd] = bilinear2(Ass, Bss, Css, Dss, xarg1, xarg2)
  % Transform all representations to transfer function (coefficient form)

  prewarp = false;  
  if nargout<=2
    if nargin==3
      B = Ass; A = Bss; fs = Css;
    elseif nargin==4
      B = Ass; A = Bss; fs = Css; fp = Dss; prewarp = true;
    else
      print_usage;
    endif
  elseif nargout==3  
    if nargin==4
      [B, A] = zp2tf(Ass, Bss, Css);
      fs = Dss;
    elseif nargin==5
      [B, A] = zp2tf(Ass, Bss, Css);
      fs = Dss; fp = xarg1; prewarp = true;
    else
      print_usage;
    endif
  elseif nargout==4
    if nargin==5
      [B, A] = ss2tf(Ass, Bss, Css, Dss);
      fs = xarg1;
    elseif nargin==6
      [B, A] = ss2tf(Ass, Bss, Css, Dss);
      fs = xarg1; fp = xarg2; prewarp = true;
    else
      print_usage;
    endif
  else
    print_usage;
  endif  
  
  % Initialization
  B = B(:)'; A = A(:)'; 
  if prewarp
    C = 2*pi*fp/tan(pi*fp/fs);
  else
    C = 2*fs;
  endif
  M = length(B);  N = length(A);
  % Numerator and denominator of BLT
  biltn = [1  -1];  biltd = [1  1];
  
  % Transform numerator (B) coefficients
  if M > 1
    Bt = zeros(1, M);
    for n = 1:M
      zp = 1;
      for k = 1:(n-1)
        zp = conv(zp, biltd);
      end
      for k = 1:(M-n)
        zp = conv(zp, C*biltn);
      end
      Bt = Bt + B(n)*zp;
    end
  else
    Bt = 1;
  end
  for k = 1:(N-M)
    Bt = conv(Bt, biltd);
  end
 
  % Transform denominator (A) coefficients
  if N > 1
    At = zeros(1, N);
    for n = 1:N
      zp = 1;
      for k = 1:(n-1)
        zp = conv(zp, biltd);
      end
      for k = 1:(N-n)
        zp = conv(zp, C*biltn);
      end
      At = At + A(n)*zp;
    end
  else
    At = 1;
  end 
  for k = 1:(M-N)
    At = conv(At, biltd);
  end
  
  %Normalize the coefficient vectors
  Bt = Bt/At(1); At = At/At(1);

  % Convert transfer function to desired form
  if nargout==3
    % Zero-pole
    [Ad, Bd, Cd] = tf2zp(Bt, At);
  elseif nargout==4
    % State-space
    [Ad, Bd, Cd, Dd] = tf2ss(Bt, At);
  else
    % Transfer function
    Ad = Bt; Bd = At;
  end
end

%!demo
%!
%! [z,p,k] = cheb1ap(10,3);
%! [A,B,C,D] = zp2ss(z,p,real(k));
%!
%! Fs = 2e3;
%! u1 = 2*Fs*tan(100*(2*pi/Fs)/2);   
%! u2 = 2*Fs*tan(500*(2*pi/Fs)/2); 
%! Bw = u2 - u1;                     
%! Wo = sqrt(u1*u2);                 
%! [At,Bt,Ct,Dt] = lp2bp(A,B,C,D,Wo,Bw);
%! [b,a] = ss2tf(At,Bt,Ct,Dt);
%!
%! figure(1)
%! clf
%! w = 0:50:10000;
%! [h] = freqs(b,a,w);     
%! plot(w,mag2db(abs(h)))
%! hold on
%! ylim([-165 5])
%! [U1,U2] = meshgrid([u1 u2],ylim);
%! plot(U1,U2)
%! legend('Magnitude response','Lower Passband Edge','Upper Passband Edge')
%! hold off
%! xlabel('Angular Frequency (rad/s)')
%! ylabel('Magnitude (dB)')
%! grid
%!
%! [Ad,Bd,Cd,Dd] = bilinear2(At,Bt,Ct,Dt,Fs);
%! [bz,az] = ss2tf(Ad,Bd,Cd,Dd);
%! [hz, fz] = freqz(bz, az, 512, Fs);
%!
%! figure(2)
%! clf
%! plot(fz, mag2db(abs(hz)))
%! ylim([-260 10])
%! grid on
%! xlabel('Frequency (Hz)')
%! ylabel('Magnitude (dB)')
%! title('Magnitude Response (dB)')

%!demo
%!
%! Fc = 20;
%! Fs = 200;                             
%! [z,p,k] = ellip(6,3,90,2*pi*Fc,'s');
%! [num,den] = zp2tf(z,p,k);
%!
%! f = 0:0.5:50;
%! [h] = freqs(num,den,2*pi*f);
%! figure(1)
%! clf
%! plot(f,mag2db(abs(h)))
%! hold on
%! xlim([0 50])
%! [l1,l2] = meshgrid(Fc,[-120 0]);
%! plot(l1,l2)
%! grid
%! legend('Magnitude response','Passband Edge')
%! xlabel('Frequency (Hz)')
%! ylabel('Magnitude (dB)')
%!
%! [numd,dend] = bilinear2(num,den,Fs,20);
%! [hz, fz] = freqz(numd, dend, 512, Fs);
%!
%! figure(2)
%! clf
%! plot(fz, mag2db(abs(hz)))
%! ylim([-110 5])
%! grid on
%! xlabel('Frequency (Hz)')
%! ylabel('Magnitude (dB)')
%! title('Magnitude Response (dB)')

