% 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 && m0 && 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 && m0 && 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)';