function [h0,h1,h2] = wframe(K0,K1)
% [h0,h1,h2] = wframe(K0,K1);
% 
% % Example:
% K0 = 4; K1 = K2 = 2;
% [h0,h1,h2] = wframe(4,2);
%
% Ivan Selesnick
% selesi@taco.poly.edu


chckmode = 0;

[r0,p,q] = maxflatI(K0,K1);

v = 1;
for k = 1:K0
   v = conv(v,[1 1]);
end

h0 = conv(sfact(q),v);
h0 = sqrt(2)*h0/sum(h0);
% N = length(h0);
N = K0 + K1;         % [length(h0) is equal to K0+K1]
L = 2*ceil(N/2);
h0 = [h0, zeros(1,L-N)];

h00 = h0(1:2:end);
h01 = h0(2:2:end);

r00 = xcorr(h00);
r01 = xcorr(h01);
r02 = ((1-(L/2):(L/2)-1)==0) - r00 - r01;
h02 = sfactK(r02,K1,1);
if chckmode > 0
	xcorr(h02)-r02  
end

P = [h00; h01; h02];
[U,P0] = colfac(P);
if chckmode > 0
	P0 - 1/sqrt(2)*[1 1 0]'
end

r = 1/sqrt(2);
H0 = [
   r  r  0;
   r -r  0;
   0  0  1
];

H = mtxsyn(U,H0);
h0_orig = h0;
h00 = H{1}(1,:); h01 = H{1}(2,:); h02 = H{1}(3,:);
h10 = H{2}(1,:); h11 = H{2}(2,:); h12 = H{2}(3,:);
h20 = H{3}(1,:); h21 = H{3}(2,:); h22 = H{3}(3,:);

if chckmode > 0
	P00 = xcorr(h00,h00) + xcorr(h01,h01) + xcorr(h02,h02);
	P01 = xcorr(h00,h10) + xcorr(h01,h11) + xcorr(h02,h12);
	P02 = xcorr(h00,h20) + xcorr(h01,h21) + xcorr(h02,h22);
	P11 = xcorr(h10,h10) + xcorr(h11,h11) + xcorr(h12,h12);
	P12 = xcorr(h10,h20) + xcorr(h11,h21) + xcorr(h12,h22);
	P22 = xcorr(h20,h20) + xcorr(h21,h21) + xcorr(h22,h22);
	[P00' P11' P22'; P01' P02' P12']
end

h0 = [h00; h01]; h0 = h0(:)';
h1 = [h10; h11]; h1 = h1(:)';
h2 = [h20; h21]; h2 = h2(:)';

h0 = h0(1:N);
h1 = h1(1:N);
h2 = h2(1:N);
[h1,h2] = makez(h1,h2,N);

if chckmode > 0
	h0_orig(1:N)-h0
	chck(h0,h1,h2);
end


% ----------------------------------------------------
%    subroutines
% ----------------------------------------------------


function [U,P0] = colfac(P);
% [U,P0] = colfac(P);
% Paraunitary column factorization
% (after PPV)
%
% P  : para-unitary column vector
% U  : projection vectors
% P0 : P(z=1)

N = size(P,2);
M = size(P,1);
U = zeros(M,N-1);
for k = N:-1:2;
   u = P(:,k); u = u/norm(u);
   A1 = P-u*u'*P;
   A2 = u*u'*P;
   P = A1(:,1:k-1) + A2(:,2:k);
   U(:,k-1) = u;
end
P0 = P;


% ----------------------------------------------------


function H = mtxsyn(U,H0)
% H = mtxsyn(U,H0)
% Paraunitary matrix synthesis
% (after PPV)
%
% U  : projection vectors
% H0 : H(z=1)
% H  : para-unitary matrix

M = size(H0,2);
for k = 1:M
   H{k} = colsyn(U,H0(:,k));
end


% ----------------------------------------------------


function P = colsyn(U,P0)
% P = colsyn(U,P0)
% Paraunitary column synthesis
% (after PPV)
%
% U  : projection vectors
% P0 : P(z=1)
% P  : para-unitary column vector

N = size(U,2)+1;
M = size(U,1);
P = P0;
for k = 1:N-1;
   u = U(:,k);
   P = [P-u*u'*P, zeros(M,1)] + [zeros(M,1), u*u'*P];
end