You are free to use these programs or any derivative of them for any
scientific purpose but please reference this book. Up-dated versions
of these programs and others can be found on our web page at:
http:

```
function p = psa(h,kk)
% p = psa(h,kk) calculates samples of the scaling function
% phi(t) = p by kk successive approximations from the
% scaling coefficients h. Initial iteration is a constant.
% phi_k(t) is plotted at each iteration. csb 5/19/93
%
if nargin==1, kk=11; end; % Default number of iterations
h2= h*2/sum(h); % normalize h(n)
K = length(h2)-1; S = 128; % Sets sample density
p = [ones(1,3*S*K),0]/(3*K); % Sets initial iteration
P = p(1:K*S); % Store for later plotting
axis([0 K*S+2 -.5 1.4]);
hu = upsam(h2,S); % upsample h(n) by S
for iter = 0:kk % Successive approx.
p = dnsample(conv(hu,p)); % convolve and down-sample
plot(p); pause; % plot each iteration
% P = [P;p(1:K*S)]; % store each iter. for plotting
end
p = p(1:K*S); % only the supported part
L = length(p);
x = ([1:L])/(S);
axis([0 3 -.5 1.4]);
plot(x,p); % Final plot
title('Scaling Function by Successive Approx.');
ylabel('Scaling Function');
xlabel('x');
```

```
function p = pdyad(h,kk)
% p = pdyad(h,kk) calculates approx. (L-1)*2^(kk+2) samples of the
% scaling function phi(t) = p by kk+3 dyadic expansions
% from the scaling coefficient vector h where L=length(h).
% Also plots phi_k(t) at each iteration. csb 5/19/93
%
if nargin==1, kk = 8; end % Default iterations
h2 = h*2/sum(h); % Normalize
N = length(h2); hr = h2(N:-1:1); hh = h2;
axis([0,N-1,-.5,1.4]);
MR = [hr,zeros(1,2*N-2)]; % Generater row for M0
MT = MR; M0 = [];
for k = 1:N-1 % Generate convolution and
MR = [0, 0, MR(1:3*N-4)]; % downsample matrix from h(n)
MT = [MT; MR];
end
M0 = MT(:,N:2*N-1); % M0*p = p if p samples of phi
MI = M0 - eye(N);
MJ = [MI(1:N-1,:);ones(1,N)];
pp = MJ\[zeros(N-1,1);1]; % Samples of phi at integers
p = pp(2:length(pp)-1).';
x = [0:length(p)+1]*(N-1)/(length(p)+1); plot(x,[0,p,0]); pause
p = conv(h2,p); % value on half integers
x = [0:length(p)+1]*(N-1)/(length(p)+1); plot(x,[0,p,0]); pause
y = conv(h2,dnsample(p)); % convolve and downsample
p = merge(y,p); % interleave values on Z and Z/2
x = [0:length(p)+1]*(N-1)/(length(p)+1); plot(x,[0,p,0]); pause
for k = 1:kk
hh = upsample(hh); % upsample coefficients
y = conv(hh,y); % calculate intermediate terms
p = merge(y,p); % insert new terms between old
x = [0:length(p)+1]*(N-1)/(length(p)+1); plot(x,[0,p,0]); pause;
end
title('Scaling Function by Dyadic Expansion');
ylabel('Scaling Function');
xlabel('x');
axis;
```

```
function [hf,ht] = pf(h,kk)
% [hf,ht] = pf(h,kk) computes and plots hf, the Fourier transform
% of the scaling function phi(t) using the freq domain
% infinite product formulation with kk iterations from the scaling
% function coefficients h. Also calculates and plots ht = phi(t)
% using the inverse FFT csb 5/19/93
if nargin==1, kk=8; end % Default iterations
L = 2^12; P = L; % Sets number of sample points
hp = fft(h,L); hf = hp; % Initializes iteration
plot(abs(hf));pause; % Plots first iteration
for k = 1:kk % Iterations
hp = [hp(1:2:L), hp(1:2:L)]; % Sample
hf = hf.*hp/sqrt(2); % Product
plot(abs(hf(1:P/2)));pause; % Plot Phi(omega) each iteration
P=P/2; % Scales axis for plot
end;
ht = real(ifft(hf)); % phi(t) from inverse FFT
ht = ht(1:8*2^kk); plot(ht(1:6*2^kk)); % Plot phi(t)
```

```
function hn = daub(N2)
% hn = daub(N2)
% Function to compute the Daubechies scaling coefficients from
% her development in the paper, "Orthonormal bases of compactly
% supported wavelets", CPAM, Nov. 1988 page 977, or in her book
% "Ten Lectures on Wavelets", SIAM, 1992 pages 168, 216.
% The polynomial R in the reference is set to zero and the
% minimum phase factorization is used.
% Not accruate for N > 20. Check results for long h(n).
% Input: N2 = N/2, where N is the length of the filter.
% Output: hn = h(n) length-N min phase scaling fn coefficients
% by rag 10/10/88, csb 3/23/93
a = 1; p = 1; q = 1; % Initialization of variables
hn = [1 1]; % Initialize factors of zeros at -1
for j = 1:N2-1,
hn = conv(hn,[1,1]); % Generate polynomial for zeros at -1
a = -a*0.25*(j+N2-1)/j; % Generate the binomial coeff. of L
p = conv(p,[1,-2,1]); % Generate variable values for L
q = [0 q 0] + a*p; % Combine terms for L
end;
q = sort(roots(q)); % Factor L
hn = conv(hn,real(poly(q(1:N2-1)))); % Combine zeros at -1 and L
hn = hn*sqrt(2)/(sum(hn)); % Normalize
```

```
function h = h246(a,b)
% h = h246(a,b) generates orthogonal scaling function
% coefficients h(n) for lengths 2, 4, and 6 using
% Resnikoff's parameterization with angles a and b.
% csb. 4/4/93
if a==b, h = [1,1]/sqrt(2); % Length-2
elseif b==0
h0 = (1 - cos(a) + sin(a))/2; % Length-4
h1 = (1 + cos(a) + sin(a))/2;
h2 = (1 + cos(a) - sin(a))/2;
h3 = (1 - cos(a) - sin(a))/2;
h = [h0 h1 h2 h3]/sqrt(2);
else % Length-6
h0 = ((1+cos(a)+sin(a))*(1-cos(b)-sin(b))+2*sin(b)*cos(a))/4;
h1 = ((1-cos(a)+sin(a))*(1+cos(b)-sin(b))-2*sin(b)*cos(a))/4;
h2 = (1+cos(a-b)+sin(a-b))/2;
h3 = (1+cos(a-b)-sin(a-b))/2;
h4 = (1-h0-h2);
h5 = (1-h1-h3);
h = [h0 h1 h2 h3 h4 h5]/sqrt(2);
end
```

```
function [a,b] = ab(h)
% [a,b] = ab(h) calculates the parameters a and b from the
% scaling function coefficient vector h for orthogonal
% systems of length 2, 4, or 6 only. csb. 5/19/93.
%
h = h*2/sum(h); x=0; % normalization
if length(h)==2, h = [0 0 h 0 0]; x=2; end;
if length(h)==4, h = [0 h 0]; x=4; end;
a = atan2((2*(h(1)^2+h(2)^2-1)+h(3)+h(4)),(2*h(2)*(h(3)-1)-2*h(1)*(h(4)-1)));
b = a - atan2((h(3)-h(4)),(h(3)+h(4)-1));
if x==2, a = 1; b = 1; end;
if x==4, b = 0; end;
```

```
function y = upsample(x)
% y = upsample(x) inserts zeros between each term in the row vector x.
% for example: [1 0 2 0 3 0] = upsample([1 2 3]). csb 3/1/93.
L = length(x);
y(:) = [x;zeros(1,L)]; y = y.';
y = y(1:2*L-1);
```

```
function y = upsam(x,S)
% y = upsam(x,S) inserts S-1 zeros between each term in the row vector x.
% for example: [1 0 2 0 3 0] = upsample([1 2 3]). csb 3/1/93.
L = length(x);
y(:) = [x;zeros(S-1,L)]; y = y.';
y = y(1:S*L-1);
```

```
function y = dnsample(x)
% y = dnsample(x) samples x by removing the even terms in x.
% for example: [1 3] = dnsample([1 2 3 4]). csb 3/1/93.
L = length(x);
y = x(1:2:L);
```

```
function z = merge(x,y)
% z = merge(x,y) interleaves the two vectors x and y.
% Example [1 2 3 4 5] = merge([1 3 5],[2 4]).
% csb 3/1/93.
%
z = [x;y,0];
z = z(:);
z = z(1:length(z)-1).';
```

```
function w = wave(p,h)
% w = wave(p,h) calculates and plots the wavelet psi(t)
% from the scaling function p and the scaling function
% coefficient vector h.
% It uses the definition of the wavelet. csb. 5/19/93.
%
h2 = h*2/sum(h);
NN = length(h2); LL = length(p); KK = round((LL)/(NN-1));
h1u = upsam(h2(NN:-1:1).*cos(pi*[0:NN-1]),KK);
w = dnsample(conv(h1u,p)); w = w(1:LL);
xx = [0:LL-1]*(NN-1)/(LL-1);
axis([1 2 3 4]); axis;
plot(xx,w);
```

```
function g = dwt(f,h,NJ)
% function g = dwt(f,h,NJ); Calculates the DWT of periodic g
% with scaling filter h and NJ scales. rag & csb 3/17/94.
%
N = length(h); L = length(f);
c = f; t = [];
if nargin==2, NJ = round(log10(L)/log10(2)); end; % Number of scales
h0 = fliplr(h); % Scaling filter
h1 = h; h1(1:2:N) = -h1(1:2:N); % Wavelet filter
for j = 1:NJ % Mallat's algorithm
L = length(c);
c = [c(mod((-(N-1):-1),L)+1) c]; % Make periodic
d = conv(c,h1); d = d(N:2:(N+L-2)); % Convolve & d-sample
c = conv(c,h0); c = c(N:2:(N+L-2)); % Convolve & d-sample
t = [d,t]; % Concatenate wlet coeffs.
end;
g = [c,t]; % The DWT
```

```
function f = idwt(g,h,NJ)
% function f = idwt(g,h,NJ); Calculates the IDWT of periodic g
% with scaling filter h and NJ scales. rag & csb 3/17/94.
%
L = length(g); N = length(h);
if nargin==2, NJ = round(log10(L)/log10(2)); end; % Number of scales
h0 = h; % Scaling filter
h1 = fliplr(h); h1(2:2:N) = -h1(2:2:N); % Wavelet filter
LJ = L/(2^NJ); % Number of SF coeffs.
c = g(1:LJ); % Scaling coeffs.
for j = 1:NJ % Mallat's algorithm
w = mod(0:N/2-1,LJ)+1; % Make periodic
d = g(LJ+1:2*LJ); % Wavelet coeffs.
cu(1:2:2*LJ+N) = [c c(1,w)]; % Up-sample & periodic
du(1:2:2*LJ+N) = [d d(1,w)]; % Up-sample & periodic
c = conv(cu,h0) + conv(du,h1); % Convolve & combine
c = c(N:N+2*LJ-1); % Periodic part
LJ = 2*LJ;
end;
f = c; % The inverse DWT
```

```
function r = mod(m,n)
% r = mod(m,n) calculates r = m modulo n
%
r = m - n*floor(m/n); % Matrix modulo n
```

```
function g = dwt5(f,h,NJ)
% function g = dwt5(f,h,NJ)
% Program to calculate the DWT from the L samples of f(t) in
% the vector f using the scaling filter h(n).
% csb 3/20/94
%
N = length(h);
c = f; t = [];
if nargin==2
NJ = round(log10(L)/log10(2)); % Number of scales
end;
h1 = h; h1(1:2:N) = -h1(1:2:N); % Wavelet filter
h0 = fliplr(h); % Scaling filter
for j = 1:NJ % Mallat's algorithm
L = length(c);
d = conv(c,h1); % Convolve
c = conv(c,h0); % Convolve
Lc = length(c);
while Lc > 2*L % Multi-wrap?
d = [(d(1:L) + d(L+1:2*L)), d(2*L+1:Lc)]; % Wrap output
c = [(c(1:L) + c(L+1:2*L)), c(2*L+1:Lc)]; % Wrap output
Lc = length(c);
end
d = [(d(1:N-1) + d(L+1:Lc)), d(N:L)]; % Wrap output
d = d(1:2:L); % Down-sample wlets coeffs.
c = [(c(1:N-1) + c(L+1:Lc)), c(N:L)]; % Wrap output
c = c(1:2:L); % Down-sample scaling fn c.
t = [d,t]; % Concatenate wlet coeffs.
end % Finish wavelet part
g = [c,t]; % Add scaling fn coeff.
```

```
function a = choose(n,k)
% a = choose(n,k)
% BINOMIAL COEFFICIENTS
% allowable inputs:
% n : integer, k : integer
% n : integer vector, k : integer
% n : integer, k : integer vector
% n : integer vector, k : integer vector (of equal dimension)
nv = n;
kv = k;
if (length(nv) == 1) & (length(kv) > 1)
nv = nv * ones(size(kv));
elseif (length(nv) > 1) & (length(kv) == 1)
kv = kv * ones(size(nv));
end
a = nv;
for i = 1:length(nv)
n = nv(i);
k = kv(i);
if n >= 0
if k >= 0
if n >= k
c = prod(1:n)/(prod(1:k)*prod(1:n-k));
else
c = 0;
end
else
c = 0;
end
else
if k >= 0
c = (-1)^k * prod(1:k-n-1)/(prod(1:k)*prod(1:-n-1));
else
if n >= k
c = (-1)^(n-k)*prod(1:-k-1)/(prod(1:n-k)*prod(1:-n-1));
else
c = 0;
end
end
end
a(i) = c;
end
```