
%% Compute Coefficients for Rational Z-Domain Transfer Function
%
% Computes the numerator/denominator coefficients for the fourth order 
% transfer function:
% 
% Causal Portion (+):
% F(z^-1) = numer/denom
% numer   = n0 + n1*z^-1 + n2*z^-2 + n3*z^-3
% denom   = 1  + d1*z^-1 + d2*z^-2 + d3*z^-3 + d4*z^-4
%
% Non-Causal Portion (-):
% F2(z^-1) = numer/denom
% numer    = n1*z + n2*z^2 + n3*z^3 + n4^z^4
% denom    = 1  + d1*z^1 + d2*z^2 + d3*z^3 + d4*z^4

% given the following impulse response:
%
% h[n] =   (a0*cos(w0*n/s) + a1*sin(w0*n/s))*exp(-b0*n/s)
%        + (c0*cos(w1*n/s) + c1*sin(w1*n/s))*exp(-b1*n/s)
%
%References:
%1] Deriche, Rachid. "Recursively implementating the Gaussian and 
%   its derivatives." (1993): 24.


%% Input Arguments

% param - Impulse response parameters (see description above).
%         1x8 vector a0,a1,c0,c1,w0,w1,b0,b1 of scalars (real).

% sigma - Sigma (standard deviation of gaussian)
%         Scalar (>0).

% causality - Optional specification of causal vs. noncausal sequence.
%             If 'noncausal' is chosen then computes coefficients for
%             non-causal portion of impulse response. It's assumed the
%             input parameters correspond to fitting to causal portion of
%             (symmetric ) gaussian or (antisymmetric) derivative of
%             gaussian.
%             String {default:'causal', 'noncausal'}

% symmetry - Optional specify whether function to approximate is symmetric
%            antisymmetric.
%            String {default: 'symmetric', 'antisymmetric'}

%% Output Arguments

% cNum - Causal numerator coefficients.
%        1x4 vector [n0,n1,n2,n3].

% cDem - Causal denominator coefficients.
%        1x4 vector [d1,d2,d3,d4].

% cNum - Non-Causal numerator coefficients.
%        1x4 vector [n1,n2,n3,n4].

% cDem - Non-Causal denominator coefficients.
%        1x4 vector [d1,d2,d3,d4].

% scale - Normalization factors for:
%         1)smoothing filter 2)derivative filter
%         1x2 vector Scalar.

%% Compute Coefficients

function [cNum, cDem, ncNum, ncDem, scale] = deriche_compute_coefs(param, sigma, varargin)

%parse input parameters
[a0,a1,c0,c1,w0,w1,b0,b1,sigma,isCausal,isSymmetric] = parse_inputs(param,sigma, varargin{:});

%temporary variables
eb0   = exp(-b0/sigma);
eb1   = exp(-b1/sigma);
sinw1 = sin(w1/sigma);
sinw0 = sin(w0/sigma);
cosw1 = cos(w1/sigma);
cosw0 = cos(w0/sigma);

%numerator coefficients
n0 = a0 + c0;
n1 = eb1*( c1*sinw1 - (c0 + 2*a0)*cosw1 ) + eb0*( a1*sinw0 - (2*c0 + a0)*cosw0);
n2 = 2*eb0*eb1*( (a0 + c0)*cosw1*cosw0 - a1*cosw1*sinw0 - c1*cosw0*sinw1 ) + c0*eb0^2 + a0*eb1^2;
n3 = eb1*eb0^2*( c1*sinw1 - c0*cosw1 ) + eb0*eb1^2*( a1*sinw0 - a0*cosw0 );

cNum = [n0,n1,n2,n3];

%denominator coefficients
d1 = -2*eb1*cosw1 - 2*eb0*cosw0;
d2 = 4*cosw1*cosw0*eb0*eb1 + eb1^2 + eb0^2;
d3 = -2*cosw0*eb0*eb1^2 - 2*cosw1*eb1*eb0^2;
d4 = eb0^2*eb1^2;

cDem = [d1,d2,d3,d4];

%determine non-causal coefficient if necessary
if isCausal
    
    ncNum = [];
    ncDem = [];
    
else
   
    ncDem      = cDem;
    ncNum(1:3) = cNum(2:4) - n0*cDem(1:3);
    ncNum(4)   = -n0*cDem(4);
        
    if ~isSymmetric        
       
        ncNum = -ncNum;
        
    end   
    
end

%compute normalization factors
scale = zeros(1,2);
scale(1) = -a0*( eb0^(-2) - 1 )/( 2*cosw0/eb0 - 1 - eb0^(-2) );
scale(2) = 2*a1*eb0^(-1)*( cosw0^2*eb0^(-2) - eb0^(-2) - cosw0^2 + 1 )/...
           (sinw0*( 4*cosw0*eb0^(-1) - 1 - 2*eb0^(-2) - 4*cosw0^2*eb0^(-2) + 4*cosw0*eb0^(-3) - eb0^(-4) )); 

end%end of deriche compute coefficients


%% Input Parser

function [a0, a1, c0, c1, w0, w1, b0, b1, sigma, isCausal, isSymmetric] = parse_inputs(param, sigma, varargin)

% parse input arguments
p = inputParser;
p.addRequired('param', @(x) ...
				validateattributes(x,{'numeric'},{'real','vector','size',[1,8],'nonnan'})); 
p.addRequired('sigma', @(x) ...
				validateattributes(x,{'numeric'},{'real','scalar','positive'}));   
p.addParameter('causality', 'causal', @(x) ischar(x) && (strcmp(x,'causal') || strcmp(x,'noncausal')));
p.addParameter('symmetry', 'symmetric', @(x) ischar(x) && (strcmp(x,'symmetric') || strcmp(x,'antisymmetric')));				       
p.parse(param,sigma, varargin{:});
inputs    = p.Results;
param     = inputs.param;
sigma     = inputs.sigma;
causality = inputs.causality;
symmetry  = inputs.symmetry;

a0 = param(1);
a1 = param(2);
c0 = param(3);
c1 = param(4);
w0 = param(5);
w1 = param(6);
b0 = param(7);
b1 = param(8);

if strcmp(causality, 'causal')
    isCausal = true;
else
    isCausal = false;
end    

if strcmp(symmetry, 'symmetric')
    isSymmetric = true;
else
    isSymmetric = false;
end    


end %end of parse inputs



