# Connexions

You are here: Home » Content » Appendix C

### Recently Viewed

This feature requires Javascript to be enabled.

# Appendix C

Module by: C. Sidney Burrus. E-mail the author

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:////www-dsp.rice.edu/

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
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


## Content actions

PDF | EPUB (?)

### What is an EPUB file?

EPUB is an electronic book format that can be read on a variety of mobile devices.

My Favorites (?)

'My Favorites' is a special kind of lens which you can use to bookmark modules and collections. 'My Favorites' can only be seen by you, and collections saved in 'My Favorites' can remember the last module you were on. You need an account to use 'My Favorites'.

| A lens I own (?)

#### Definition of a lens

##### Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

##### What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

##### Who can create a lens?

Any individual member, a community, or a respected organization.

##### What are tags?

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

| External bookmarks

### Reuse / Edit:

Reuse or edit module (?)

#### Check out and edit

If you have permission to edit this content, using the "Reuse / Edit" action will allow you to check the content out into your Personal Workspace or a shared Workgroup and then make your edits.

#### Derive a copy

If you don't have permission to edit the content, you can still use "Reuse / Edit" to adapt the content by creating a derived copy of it and then editing and publishing the copy.