%% 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{AT BT CT DT} =} sstrans (@var{A}, @var{B}, @var{C},@var{D}, @var{W}, @var{STOP})
%%
%% Transform band edges of a generic lowpass filter (cutoff at W=1)
%% represented in state-space form.  W is the edge of the
%% target filter (or edges if band pass or band stop). Stop is true for
%% high pass and band stop filters or false for low pass and band pass
%% filters. Filter edges are specified in radians/second.
%%
%% @seealso{lp2lp, lp2hp, lp2bp, lp2bs}
%% @end deftypefn

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

function [AT BT CT DT] = sstrans (A, B, C, D, W, stop)
  if (nargin != 6)
    print_usage;
  endif

  if length(W)==2
    Fl = W(1); Fh = W(2); BW = Fh - Fl;
    Wo = sqrt(Fl*Fh); Q = Wo/BW;
    [ma, m] = size(A); n = size(B, 2); mc = size(C, 1);
    if stop
      %% Band Stop       
      %%        S*BW
      %% S ->   -------- 
      %%        S^2+FhFl 
      AT = [Wo/Q*inv(A) Wo*eye(ma); -Wo*eye(ma) zeros(ma)];
      BT = -[Wo/Q*(A\B); zeros(ma,n)];
      CT = [C/A zeros(mc,ma)];
      DT = D - C/A*B;
    else
      %% Band Pass       
      %%        S^2+FhFl 
      %% S ->   -------- 
      %%        S*BW
      AT = Wo*[A/Q eye(ma,m); -eye(ma,m) zeros(ma,m)];
      BT = Wo*[B/Q; zeros(ma,n)];
      CT = [C zeros(mc,ma)];
      DT = D;
    endif  
  else
    Fc = W;
    if stop
      %% High Pass
      %% S ->   Fc/S     
      AT = W*inv(A); BT = -W*(A\B); CT = C/A; DT = D-C/A*B;
    else
      %% Low Pass
      %% S ->   S/Fc     
      AT = W*A; BT = W*B; CT = C; DT = D;
    endif
  endif
endfunction
