% David Ferrone, January 2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Produces the matrices involved in the periodized biorthogonal transform.
% Make sure that s_dual contains more non-trivial coefficients,
% s_prim & s_dual are the same length (add 0's to the end of s_prim),
% s_prim & s_dual should be biorthogonal.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The idea is to produce banded Q = UTV where T is unbanded periodic-form of
% the discrete wavelet transform.
% This constructs the matrices V_prim and V_dual.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If not using BiMRAsetup, set the coefficients and signal length here.
%n_f=11; % Matrices will be constructed for signals of length 2n_f.
%s_dual=[0,2,0,0]; % Primal scaling coefficients
%s_prim=[0,0,1/2,1,1/2,0]; % Dual scaling coefficients
%s_dual=[0,0,0,0,.25, .75, .75, .25,0,0,0,0];
%s_prim=[-.01953125,.05859375,.07421875,-.37890625,-.1015625,1.3671875, 1.3671875,-.1015625,-.37890625, .07421875, .05859375,-.01953125];
%s_prim=[0,0,0,1/2,1,1/2,0,0,0,0];
%s_dual=1/64*[3,-6,-16,38,90,38,-16,-6,3,0];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
w_prim = zeros(length(s_prim)); w_dual = zeros(length(s_prim));
n_s=(length(s_prim)-2)/2;
% There are 2n_s + 1 total scaling coefficients: s_0, ..., s_{2n_s+1}
for k=1:length(s_prim)
w_prim(k)=(-1)^k*s_dual(length(s_prim)-k+1);
w_dual(k)=(-1)^k*s_prim(length(s_prim)-k+1);
end;
S_Prim_a=Create_block_matrix_A(s_prim);
S_Prim_b=Create_block_matrix_B(s_prim);
S_Dual_a=Create_block_matrix_A(s_dual);
S_Dual_b=Create_block_matrix_B(s_dual);
W_Prim_a=Create_block_matrix_A(w_prim);
W_Prim_b=Create_block_matrix_B(w_prim);
W_Dual_a=Create_block_matrix_A(w_dual);
W_Dual_b=Create_block_matrix_B(w_dual);
Ident=eye(2*(n_f-n_s));
Z_top=zeros(2*n_s,2*(n_f-n_s));
Z_side=zeros(2*(n_f-n_s),n_s);
R_a = inv(S_Dual_a*W_Prim_a');
R_b = inv(S_Dual_b*W_Prim_b');
% inv() might not be the best approach, but at least this is readable.
V_dual = [horzcat((R_a*S_Dual_a)', Z_top, (R_b*S_Dual_b)'); horzcat(Z_side, Ident, Z_side)];
V_prim = [horzcat(W_Prim_a', Z_top, W_Prim_b'); horzcat(Z_side, Ident, Z_side)];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% T_prim, T_dual - the two periodized transformation matrices.
T_prim=Create_T_Matrix(n_f,n_s,s_prim,w_prim);
T_dual=Create_T_Matrix(n_f,n_s,s_dual,w_dual);
Q_prim=(1/sqrt(2))*T_prim*V_dual;
Q_dual=(1/sqrt(2))*T_dual*V_prim;