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