% Usage: [phi,psi,t] = dyadicortho(c, r)
% Given vector c of lowpass filter coefficients for an 
% orthogonal filter and r this computes the appropriate scaling 
% function at the dyadic % points k/2^r for 0<k<2^r(N-1), where 
% N is the filter length. The scaling function is returned at 
% points t = k/2^r on the range 0 <= k <= 2^r(N-1), in the 
% argument "phi", a vector of length 2^r(N-1)+1. The wavelet 
% is returned in "psi", a vector. The dyadic points
% at which the function were computed is returned in "t".

% Note, this routine will assume coefficients of c sum to 
% sqrt(2), as is implicit in our version of the dilation 
% equation. In fact, this routine rescales them! 

function [phi,psi,t] = dyadicortho(c,r)

r2=sqrt(2);

N = length(c);
R = 2^r;
n = (N-1)*R+1;
phi=zeros(n,1);

c = c/sum(c)*r2; %Make sure scaling is correct

%First set up the matrix to determine phi at 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) = r2*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(N-2); %Get position of largest eigenvalue
v = real(V(:,s)); %Get dominant eigenvector

%Now fill in all dyadic points k/2^r
phi(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+phi(R*p+1)*c(j+1);
            end
            phi(k*F+1)=r2*s;
        end
    end
end

%Now compute the wavelet
psi=zeros(n,1);

%Highpass coefficients
d = zeros(1,N);
for k=1:N
    d(k)=(-1)^(k+1)*c(N-k+1);
end

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

t=((0:R*(N-1))/R)';