% Usage: [phi,phi2,psi,psi2,t,t2,t3] = dyadicbiortho(c,d,r)
% Given vectors c, d of lowpass and highpass filter 
% coefficients for the analysis filter for a biorthogonal 
% filter bank, with c of length M, d of length N (and the 
% synthesis filters are assumed to be lowpass filter 
% (-1)^k*d(k), 0<=k<=N-1, highpass filter (-1)^(k+1)*c(k), 
% 0<=k<=M-1)), and r this computes both scaling functions at 
% the dyadic points t=k/2^r for 0<=k<=2^r(M-1) (analysis scaling) 
% or 0<=k<=2^r(N-1) (synthesis scaling), or 0<=k<\2^r((M+N)/2-1)
% (both wavelets). The analysis scaling function is returned in 
% "phi", wavelet in  "psi", synthesis scaling and wavelet in 
% "phi2", "psi2". The arrays t, t2, t3 contain the times at which 
% phi, phi2 and the wavelets were computed, respectively.

% Note, this routine will assume coefficients of c sum to 1 and the
% coefficients of the highpass filter d have alternating sum 1, i.e.,
% sum((-1)^k*d(k), k=0..N-1)=1.

function [phi,phi2,psi,psi2,t,t2,t3] = dyadicbiortho(c,d,r)

M = length(c);
N = length(d);
Q = (M+N)/2;
R = 2^r;
nM = (M-1)*R+1;
nN = (N-1)*R+1;
nQ = (Q-1)*R+1;
phi=zeros(nM,1);
phi2=zeros(nN,1);

c = c/sum(c); %Make sure scaling is correct: coefs sum to 1.

%Set up the matrix which determines phi the integers 1 to N-2.

A = zeros(M-2,M-2);
for p=1:(M-2)
    for k=1:M
        m = 2*p-k+1;
        if m>0 && m<M-1
            A(p,m) = 2*c(k);
        end
    end
end

%Get the eigenvector with eigenvalue 1
[V,D] = eig(A); %V holds vectors, D eigenvalues
[C,I] = sort(diag(D)); %sort eigenvalues
s = I(M-2); %Get position of largest eigenvalue
v = real(V(:,s)); %Get dominant eigenvector

%Now fill in all dyadic points k/2^r for phi
phi(1+R*(1:M-2)) = v;
for q=1:r
    Q=2^q; F=R/Q;
    for k=1:2:(Q*(M-1))
        s=0;
        for j=0:(M-1)
            p=2*k/Q-j;
            if p>0 && p<(M-1)
                s = s+phi(R*p+1)*c(j+1);
            end
            phi(k*F+1)=2*s;
        end
    end
end

%Repeat the above process for phi2, the synthesis low-pass 
%scaling filter. Get synthesis low-pass from analysis high-pass.
c2 = zeros(1,N);
for k=1:N
    c2(k)=(-1)^(k+1)*d(k);
end

c2 = c2/sum(c2); %Make sure scaling is correct: coefs sum to 1.
d = d/sum(c2);

%Set up the matrix which determines phi at the integers 1 to N-2.

A = zeros(N-2,N-2);
for p=1:(N-2)
    for k=1:N
        m = 2*p-k+1;
        if m>0 && m<N-1
            A(p,m) = 2*c2(k);
        end
    end
end

%Get the eigenvector with eigenvalue 1
[V,D] = eig(A); %V holds vectors, D eigenvalues
[C,I] = sort(diag(D)); %sort eigenvalues
s = I(N-2); %Get position of largest eigenvalue
v = real(V(:,s)); %Get dominant eigenvector

%Now fill in all dyadic points k/2^r for phi2
phi2(1+R*(1:N-2)) = v;
for q=1:r
    Q=2^q; F=R/Q;
    for k=1:2:(Q*(N-1))
        s=0;
        for j=0:(N-1)
            p=2*k/Q-j;
            if p>0 && p<(N-1)
                s = s+phi2(R*p+1)*c2(j+1);
            end
            phi2(k*F+1)=2*s;
        end
    end
end

%Now compute the wavelet psi corresponding to phi.
psi=zeros(nQ,1);

%Wavelet construction.
for m=1:nQ
    s=0;
    for k=1:N
        q=2*(m-1)-(k-1)*R+1;
        if q>0 && q<=nM
            s=s+d(k)*phi(q);
        end
    end
    psi(m)=2*s;
end

%Now compute the wavelet psi2 corresponding to phi2.
psi2=zeros(nQ,1);

%Get synthesis high-pass from analysis low-pass.
d2 = zeros(1,M);
for k=1:M
    d2(k)=(-1)^k*c(k);
end

%Wavelet construction for phi2.
for m=1:nQ
    s=0;
    for k=1:M
        q=2*(m-1)-(k-1)*R+1;
        if q>0 && q<=nN
            s=s+d2(k)*phi2(q);
        end
    end
    psi2(m)=2*s;
end

% Times at which each scaling function was computed.
t=((0:R*(M-1))/R)';
t2=((0:R*(N-1))/R)';
t3=((0:R*((M+N)/2-1))/R)';