function g = er_blaschke(f, s, phi)
% function g = er_blaschke(f, s, phi)
%   Performs desired frequency transform on prototype digital LPF.  Takes filter
%   coefficient vectors in struct F, sign specification in S (S = +/- 1), and frequency 
%   point vector in PHI.  Returns new transfer function coefficient vectors in struct G.
%
%   EXAMPLE: >> g = er_blaschke(f, 1,[0.2 0.4 0.6 0.8]), yields multi-band
%   filter from original prototype digital LPF with passband scetions between 
%   (0,0.2*pi), (0.4*pi,0.6*pi), (0.8*pi,pi).
%
%   Ref: R.A. Roberts and C.T. Mullis, Digital Signal Processing, Massachusetts: Adison-Wesley, 1987. Chapter 6.
%
%   Author: Evan Ruzanski, CU-Boulder, ECEN5632 MATLAB Assignment, FA2004

%%%%%%%%%% UNPACK INPUT STRUCT %%%%%%%%%%
num1 = f(1).tf_complete;
den1 = f(2).tf_complete;

%%%%%%%%%% FIND SPECTRAL TRANSFORMATION PARAMETERS (Ref. Roberts, Mullis, p.207) %%%%%%%%%%
phi = phi*pi; % Normalize input vector to pi
n = length(phi);
% Initialize algorithm
v = 0.5;
p(1) = 1;
% Iterate algorithm as in Roberts, Mullis, p. 207
for k = 1:n
    v = -v;
    phip = (phi(k) - pi)*v;
        for j = 0:k
        alpha = 0;
        beta = 0;
        if j > 0
            alpha = alpha + p(j);
            beta = beta - p(k - j + 1);
        end
        if j < k
            alpha = alpha + p(j + 1);
            beta = beta + p(k - j);
        end
        q(j + 1) = alpha*cos(phip) + beta*sin(phip);
    end
    for j = 0:k
        p(j + 1) = q(j + 1);
    end
end

%%%%%%%%%% PERFORM SPECTRAL TRANSFORMATION %%%%%%%%%%
m = 0; % No shift
A = p;
n = length(num1) - 1; % Order of numerator, from given filter, assume order of numerator = order of denominator
for k = 1:n + 1
    A = s*A;
    % Step 1: Handle exponents, first allpass function in equation
    if (k == n + 1)
        t1 = 1;
    elseif (k == n)
        t1 = A;
    else
        t1 = A;
        for q = 1:(n - k)
            tt1 = conv(t1,A);
            t1 = tt1;
        end
    end
    A = A/s;
    % Step 2: Handle exponents, second reversed allpass function in
    % equation
    if (k == 1)
        t2 = 1;
    elseif (k == 2)
        t2 = fliplr(A);
    else
        t2 = fliplr(A);
        for p = 1:(k - 2)
            tt2 = conv(t2,fliplr(A));
            t2 = tt2;
        end
    end
    % Step 3: Combine convolution results
    AAA = conv(t1,t2);
    % Step 4: Shift results and zeropad accordingly
    BBB = [zeros(1,m*((n + 1) - k)) AAA];
    AAA = [zeros(1,m*((n + 1) - k)) AAA];
    % Step 5: Multiply by appropriate coefficient
    F(k,1:length(BBB)) = num1(k)*BBB;
    G(k,1:length(AAA)) = den1(k)*AAA;
end
% Step 6: Create polynomials
ss = size(G);
for k = 1:ss(2)
    bspec(k) = sum(F(:,k));
    aspec(k) = sum(G(:,k));
end
bspec = bspec/aspec(1);
aspec = aspec/aspec(1);

%%%%%%%%%% REPACK OUTPUT STRUCT %%%%%%%%%%
f1(1).tf_complete = bspec;
f1(2).tf_complete = aspec; 
g = f1;