function g = er_fconv(f1,f2)
% function g = er_fconv(f1,f2)
%   Function to convolve given transfer function polynomials in the frequency
%   domain.  Takes two input structs F1 and F2 and returns result in struct G.
%   Uses algorithm described in Sec. 7, pp. 1-6 in ECEN5632 notes by C.T.
%   Mullis.  This program works for causal transfer functions based on given 
%   parameterization. 
%   
%   Author:  Evan Ruzanski, CU-Boulder, ECEN5632 MATLAB assignment, FA2004

% Unpack input struct
num1 = f1(1).tf_complete;
den1 = f1(2).tf_complete;
num2 = f2(1).tf_complete;
den2 = f2(2).tf_complete;
L = (length(den1)-1)*(length(den2)-1);
x = deconvolve(den1,num1,L); % Perform deconvolution on a1, b1 s.t. x(k) = f(k)
y = deconvolve(den2,num2,L); % Perform deconvolution on a2, b2 s.t. y(k) = g(k)
x = x.*y; % x(k) = h(k)
ahat = join(den1,den2); % Join denominators a1, a2 using Newton moments
for k = 1:L + 1 % Create numerator
    summ = 0;
    for m = 1:k
        summ = summ + ahat(m)*x(k + 1 - m);
    end
    bhat(k) = summ;
end
numcconv = bhat;
dencconv = ahat;
% Repack struct for output
f1(1).tf_complete = numcconv;
f1(2).tf_complete = dencconv; 
g = f1;

%%%%% DECLARE LOCAL FUNCTIONS %%%%%
function y = deconvolve(a, b, L)
% DECONVOLVE  Deconvolve polynomials a and b
y(1:L + 1) = 0; % Initialize output
for k = 1:L + 1
    if (k == 1)
        y(k) = b(k);
    elseif ((k >= 2) & (k <= length(b)))
        sum1 = 0;
        for m = 2:k
            sum1 = sum1 + a(m)*y(k + 1 - m);
        end
        y(k) = b(k) - sum1;
    elseif (k > length(b))
        sum2 = 0;
        for m = 2:length(b)
            sum2 = sum2 + a(m)*y(k + 1 - m);
        end
        y(k) = -sum2;
    end
end

function ah = join(a, b)
% JOIN  Joins given polynomials using Newton moments
L = (length(a)-1)*(length(b)-1);
etaf = f2mom(a,L);
etag = f2mom(b,L);
etah = etaf.*etag;
ah = mom2f(etah);

function eta = f2mom(a, L)
% F2MOM  Convert polynomial to Newton moments
n = length(a) - 1; % Order of a
eta(1:L + 1) = 0;
for k = 1:L + 1
    if (k == 1)
        eta(k) = n;
    elseif (k <= n + 1)
        sum1 = 0;
        for m = 2:k
            sum1 = sum1 + a(m)*eta(k + 1 - m);
        end
        y(k) = -k*a(k) - sum1;
    elseif (k > n + 1)
        sum2 = 0;
        for m = 2:n + 1
            sum1 = sum1 + a(m)*eta(k + 1 - m);
        end
        eta(k) = -sum2;
    end
end

function ahat = mom2f(etah)
% MOM2F  Convert Newton moments to polynomial
ahat(1:(etah(1)+1)) = 0;
for k = 1:length(ahat) - 1
    if k == 1
        ahat(k) = 1;
    else
        sum3 = 0;
        for m = 1:k
            sum3 = sum3 + ahat(m)*etah(k + 1 - m);
        end
        ahat(k) = -1/k*sum3;
    end
end