Skip to content Skip to navigation

OpenStax_CNX

You are here: Home » Content » Appendix C

Navigation

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

Content actions

Download module as:

PDF | EPUB (?)

What is an EPUB file?

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

Downloading to a reading device

For detailed instructions on how to download this content's EPUB to your specific device, click the "(?)" link.

| More downloads ...

Add module to:

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? tag icon

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