function [A,B,C,D,p]= butter_synth_ss( N, fc, fs, lh);
% [A,B,C,D]= butter_synth_ss( N, fc, fs, lh);
% INPUT     DESCRIPTION
%   N        filter order
%   fc       -3 dB frequency in Hz
%   fs        sample frequency in Hz   (optional)
%   lh       'l' for low-pass filter (default) and 'h' for high-pass filter
%
% OUTPUT    DESCRIPTION
% A,B,C,D    state-space realization 
%    p       poles of the continuous time system
%    f  fs  is included,  (A,B,C,D)  is in discrete time with first order hold
%    otherwise,  (A,B,C,D)  is in companion matrix continuous time
 
% https://www.dsprelated.com/showarticle/1119.php

if nargin < 4, lh = 'l'; end
 
% poles of the Butterworth filter in continuous time
p = 2*pi*fc * exp(1i*pi*(0.5 + (2*[1:N]-1)/(2*N)));

a = flip(real(poly(p)));         % coeff's of denominator polynomial

A = [ zeros(N-1,1) , eye(N-1) ; -a(1:N) ];    % dynamics matrix

switch lh

 case 'h' % continuous-time state-space realization of high-pass filter
   B = [ zeros(N-1,1) ;  1 ];                         %       input matrix
   C = [ -a(1:N) ];                                   %      output matrix
   D = [  1 ];                                        % feedthrough matrix

 otherwise  % continuous-time state-space realization of  low-pass filter
   B = [ zeros(N-1,1) ;  1 ];                         %       input matrix
   C = [ a(1) , zeros(1,N-1) ];                       %      output matrix
   D = [ 0 ];                                         % feedthrough matrix

end

% discrete-time state-space realization of high-pass Butterworth filter
if size(fs,1) > 0
  if fc >= fs/2, error('fc must be less than fs/2'); end
  [ A , B , C , D ] = con2dis(A,B,C,D,1/fs,'foh');
end

% w = 2*pi*logspace(-3,2,50); 
% mimoBode(A,B,C,D,w,[],100,'b');
% mimoBode(A,B,C,D,w,1/fs,200,'b');
 
% ===========================================================================
% HP Gavin, 2021-07-19
