wavelet.m

小波变换所有常见例程

  g0 = [0,a;0,sqrt(2)*a];
  g1 = [b,c;sqrt(2)*b,sqrt(2)*c];
  g2 = [-1,c;0,-sqrt(2)*c];
  g3 = [b,a;-sqrt(2)*b,-sqrt(2)*a];
  H = mpoly({[h0;g0],[h1;g1],[h2;g2],[h3;g3]},-2,'symbol',2,2); 
  H = H/2;

  ht0 = [0,at;0,0];
  ht1 = [bt,ct;0,0];
  ht2 = [1,ct;0,dt];
  ht3 = [bt,at;et,dt];
  gt0 = [0,at;0,sqrt(2)*at];
  gt1 = [bt,ct;sqrt(2)*bt,sqrt(2)*ct];
  gt2 = [-1,ct;0,-sqrt(2)*ct];
  gt3 = [bt,at;-sqrt(2)*bt,-sqrt(2)*at];
  Ht = mpoly({[ht0;gt0],[ht1;gt1],[ht2;gt2],[ht3;gt3]},-2,'symbol',2,2); 
  Ht = Ht/2;

 case 'jrzb'
%  Full form: 'jrzb', s, t, lambda, mu
%  Jia-Riemenschneider-Zhou biorthogonal, |2 lambda + mu| < 2.
%  scaling function only; wavelet functions / dual not known
%  symmetric, m=2, r=2, l=3, support [0,2], p=1
%
%  R.-Q. Jia, S. D. Riemenschneider, and D.-X. Zhou, Vector
%        Subdivision Schemes and Multiple Wavelets, Math. Comp. 67 (1998),
%        pp. 1533 - 1563

  s = varargin{1};
  t = varargin{2};
  lambda = varargin{3};
  mu = varargin{4};

  h0 = [1/4,s/4;t/2,lambda/2];
  h1 = [1/2,0;0,mu/2];
  h2 = [1/4,-s/4;-t/2,lambda/2];
  H = mpoly({h0,h1,h2},0,'symbol',2,2); 

 otherwise
  error('unknown type of wavelet');
end

if symbolic
    return
end

H = double(H);
if (nargout >= 2)
    Ht = double(Ht);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                       %
%  End of main program                                                  %
%  The rest are routines that calculate some special types of wavelets  % 
%  (Daubechies, Cohen)                                                  %
%                                                                       %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function S = daubechies(p,type)

% DAUBECHIES --  produce symbol of Daubechies wavelet
%
%        H = daubechies(p,type)
%
% Calculates the coefficients for the Daubechies wavelet with p
% vanishing moments, and returns them as a symbol. If type = 1 (or
% omitted), this routine calculates the original Daubechies wavelets
% as given in
%  
%        Ingrid Daubechies
%        Orthonormal Bases of Compactly Supported Wavelets
%        Comm. Pure Appl. Math. 41 (1988), p. 909 - 996
%
% Otherwise, it calculates the least asymmetric Daubechies
% wavelets, as defined in section 2 of
%
%        Ingrid Daubechies
%        Orthonormal Bases of Compactly Supported Wavelets II:
%            Variations on a Theme
%        Siam J. Math. Anal., vol. 24, no. 2, March 1993
%        p. 499 - 519
%        
% For p <= 3, the two types are identical.      
%
% The algorithm is that given in chapter 6 of Daubechies' book "Ten
% Lectures on Wavelets" (TLW).
%  
% In all cases, there are really two solutions, one of them being the
% other one in reverse. I have not been able to find any standard
% choice in the literature or in other people's programs. My choice
% is to select the ordering so that the largest coefficient occurs in
% the first half. The peak and center of mass of the scaling function
% are therefore also to the left of center. If you are trying to
% duplicate some results from the literature, make sure to check that
% you are using the same wavelet. Likewise, someone else's
% coefficients may have a minus sign attached to the high-pass
% coefficients.
%
% For p = 1, 2, 3, the coefficients are returned in symbolic form.
% For larger p, the coefficients are always numeric.

if (round(p) ~= p | p < 1)
    error('p must be an integer >=1');
end

if (nargin < 2)
    type = 1;
end

if (p == 1)
    S = mpoly({sym('[1/2;1/2]'),sym('[1/2;-1/2]')},0,'symbol',2,1);
    return;
end

if (p == 2)
    S = mpoly({sym('[1+sqrt(3);1-sqrt(3)]'),...
	       sym('[3+sqrt(3);-3+sqrt(3)]'),...
	       sym('[3-sqrt(3);3+sqrt(3)]'),...
	       sym('[1-sqrt(3);-1-sqrt(3)]')},0,'symbol',2,1);
    S = S/8;
    return
end

if (p == 3)
    S = mpoly({sym('[1+sqrt(10)+sqrt(5+2*sqrt(10));1+sqrt(10)-sqrt(5+2*sqrt(10))]'),...
	       sym('[5+sqrt(10)+3*sqrt(5+2*sqrt(10));-5-sqrt(10)+3*sqrt(5+2*sqrt(10))]'),...
	       sym('[10-2*sqrt(10)+2*sqrt(5+2*sqrt(10));10-2*sqrt(10)-2*sqrt(5+2*sqrt(10))]'),...
	       sym('[10-2*sqrt(10)-2*sqrt(5+2*sqrt(10));-10+2*sqrt(10)-2*sqrt(5+2*sqrt(10))]'),...
	       sym('[5+sqrt(10)-3*sqrt(5+2*sqrt(10));5+sqrt(10)+3*sqrt(5+2*sqrt(10))]'),...
	       sym('[1+sqrt(10)-sqrt(5+2*sqrt(10));-1-sqrt(10)-sqrt(5+2*sqrt(10))]')},...
	       0,'symbol',2,1);
    S = S/32;
    return
end

%  Compute polynomial P = P_p(y) (eq. 6.1.12 in TLW)

k = 0:p-1;

P = mpoly(reshape(binomial(p-1+k,k),1,1,p));

%  Convert P_p(y), y = sin^2(w/2) to form Q(z), z = exp(iw)

y = mpoly({-1/4,1/2,-1/4},-1);
Q = P(y);

%  Find the roots of Q

r = roots(Q.coef(:));

%  Now we need to select one root out of each real pair, or one pair
%  out of each complex quadruple.
%
%  If type == 1, we choose the roots outside the unit disk.
%  If type ~= 1, we have to generate all possible choices, and
%  calculate the deviation from symmetry for each one.

if (type == 1)
    r = r(find(abs(r)>1));
else
    %  choose one root from each group (real pair or complex quadruple)
    r = r(find(abs(r)>1 & imag(r) >= 0));
    %  generate all possibilities for inverting any element of r except
    %  the first
    R = r;
    l = length(r);
    for k = 2:l
	s = R;
	s(k,:) = 1 ./ s(k,:);
	R = [R s];
    end
    %  calculate max(abs(psi_tot)) for all of them, select the smallest
    P = zeros(1,2^(l-1));
    for k = 1:2^(l-1)
	P(k) = psi_tot(R(:,k));
    end
    [minP,mini] = min(P);
    r = R(:,mini)';
    %  put back the second representative of each complex quadruple
    l = length(r);
    for i = 1:l
	if (imag(r(i)) ~= 0)
	    r = [r conj(r(i))];
	end
    end
end

%  construct q = sqrt(Q) from selected roots
%  normalize phase so that q(0) is real

q = poly(r);
q0 = q(length(q));
q = q*(conj(q0)/abs(q0));

% Everything should be real, so drop the imaginary parts

q = real(q);

% Now compute h by multiplying with (0.5(1+z))^p

for k = 1:p,
  q = conv(q,[1/2 1/2]);
end

% Normalize

h = q/sum(q);

%  choose order so that the largest coefficient shows up in the first half

[hmax,imax] = max(abs(h));
if (imax > length(h)/2)
    h = h(length(h):-1:1);
end

% find g, assemble the result

g = h(length(h):-1:1);
g(2:2:end) = - g(2:2:end);
S = mpoly(reshape([h;g],2,1,length(h)),0,'symbol',2,1);
return

function P = psi_tot(r)

% PSI_TOT -- maximum deviation from linear phase for given choice of roots
%
%        P = psi_tot(r)
%        
% This is used to find the most symmetric Daubechies wavelets.
%
% See page 505 in
%
%        Ingrid Daubechies
%        Orthonormal Bases of Compactly Supported Wavelets II: Variations on a Theme
%        Siam J. Math. Anal. vol. 24, no. 2, March 1993
%        p. 499 - 519

npoints = 65;
xi = (0:npoints-1)/(npoints-1) * 2 * pi;
Psi = zeros(1,npoints);

for j = 1:length(r)
    if (imag(r(j)) == 0)
	temp = atan((1+r(j))/(1-r(j))*tan(xi/2));
    else
	Rl = abs(r(j));
	alphal = angle(r(j));
	temp = atan((1 - Rl^2)*sin(xi)./((1+Rl^2)*cos(xi)-2*Rl*cos(alphal)));
    end	
    %  get rid of jumps in arctangent
    for k = 2:npoints
	if (temp(k) - temp(k-1) > pi/2)
	    temp(k:npoints) = temp(k:npoints) - pi;
	end
	if (temp(k) - temp(k-1) < -pi/2)
	    temp(k:npoints) = temp(k:npoints) + pi;
	end
    end
    temp = temp - xi/(2*pi) * temp(npoints);
    Psi = Psi + temp;
end

P = max(abs(Psi));
return

function [H,Ht] = bidaubechies(m1,m2,roots1)

% BIDAUBECHIES --  coefficients of symmetric biorthogonal Daubechies wavelet
%
%        [H,Ht] = bidaubechies(m1,m2,roots)
%
% Calculates the coefficients for the biorthogonal symmetric
% Daubechies wavelets with m1 coefficients for H, m2 coefficients for
% Ht (m1 + m2 = divisible by 4), as described on p. 278 of "Ten
% Lectures on Wavelets" (TLW). The total number of vanishing moments
% between F and Ft is (m1+m2)/2.
%
% Algorithm: Form the polynomial Q(z) described in "Ten Lectures on
% Wavelets" (TLW). It is called p_A there (page 172). This polynomial
% has roots that come in real pairs (r, 1/r) and complex quadruples
% (z, z_bar, 1/z, 1/z_bar).
%
% For orthogonal wavelets (m1 = m2 = 2m), we pick one real root
% from each pair, and one complex conjugate pair from each complex
% quadruple to get the factor polynomials, and add one factor of
% (1+z)^m to each.
%
% In the present setting, we need to assign each real pair and complex
% quadruple entirely to one polynomial or the other, and we can split
% the total factor (1+z)^(m1+m2) between them any way we want.
%
% If m=(m1+m2)/4, then at least up to m=10, we have 0 real pairs if m
% is odd, 1 real pair if m is even, the rest are complex
% quadruples. We always assign the real pair to H (otherwise, just
% switch H, Ht). The complex pairs with subscripts given in the
% subscript vector ROOTS also get put into F. Then we adjust the (1+z)
% factors to get the correct lengths.
%
% I cannot guarantee that this routine will work for m > 10.

% error checking on arguments

if (nargin < 3)
    error('usage: [S,St] = bidaubechies(m1,m2,roots)');
end

m = (m1 + m2)/4;

if (round(m) ~= m)
    error('m1 + m2 must be divisible by 4');
end

if (m1 <= 0 | m2 <= 0)
    error('m1, m2 must be >= 1');
end

ncq = floor((m-1)/2);  % number of complex quadruples
if (max(roots1) > ncq)
    error(['you cannot pick more than ',num2str(ncq),' complex roots']);
end

l = (ncq - length(roots1)); % number of quadruples used in q2
n2 = m2 - (l * 4 + 1);
n1 = 2*m - n2;
if (n1 < 0 | n2 < 0)
    error(['for selected roots, m2 must be between ',num2str(l*4+1),...
	' and ',num2str(2*m+l*4+1)]);
end

%  Compute polynomial P = P_m(y) (eq. 6.1.12 in TLW)

k = 0:m-1;
P = mpoly(reshape(binomial(m-1+k,k),1,1,m));

%  Convert P_m(y), y = sin^2(w/2) to form Q(z), z = exp(iw)

y = mpoly({-1/4,1/2,-1/4},-1);
Q = P(y);

%  Find the roots of Q, sort into real pairs and complex quadruples

r = roots(Q.coef(:));
rr = r(find(abs(r) > 1 & imag(r) == 0));
rc = r(find(abs(r) > 1 & imag(r) > 0));

% select the desired complex quadruples

rc1 = rc(roots1);
roots2 = 1:length(rc);
roots2 = roots2(~ismember(roots2,roots1));
rc2 = rc(roots2);

% put back the other representatives of each real pair or complex quadruple

rr = [rr; 1./rr];
rc1 = [rc1; 1./rc1; conj(rc1); 1./conj(rc1)];
rc2 = [rc2; 1./rc2; conj(rc2); 1./conj(rc2)];

% construct q1, q2 from selected roots
% by construction, they should both be real
% zero any complex components introduced by roundoff

q1 = real(poly([rr;rc1]));
q2 = real(poly(rc2));

% Multiply with powers of (0.5(1+z))

for k = 1:n1
  q1 = conv(q1,[1/2 1/2]);
end

for k = 1:n2
  q2 = conv(q2,[1/2 1/2]);
end

% Normalize

h  = q1 / sum(q1);
ht = q2 / sum(q2);

% find duals

g = ht(length(ht):-1:1);
g(2:2:end) = - g(2:2:end);
gt = h(length(h):-1:1);
gt(1:2:end) = - gt(1:2:end);

% adjust h, ht to common center

H = mpoly(reshape(h,1,1,length(h)),-(length(h)-1)/2,'symbol',2,1);
Ht = mpoly(reshape(ht,1,1,length(ht)),-(length(ht) - 1)/2,'symbol',2,1);
% H = mpoly(reshape(h,1,1,length(h)),0,'symbol',2,1);
% Ht = mpoly(reshape(ht,1,1,length(ht)),(length(h) - length(ht))/2,'symbol',2,1);

% put result together

P = polyphase(H);
Pt = polyphase(Ht);
P = [P;Pt(1,2)',-Pt(1,1)'];
Pt = [Pt;P(1,2)',-P(1,1)'];
H = symbol(P);
Ht = symbol(Pt);

function [H,Ht] = cohen(n,ntilde)

% COHEN -  coefficients for Cohen biorthogonal wavelets
%
%        [H,Ht] = cohen(n,ntilde)
%  
% Returns the wavelet coefficients for the wavelets introduced on
% page 540 in 
%  
%       Cohen - Daubechies - Feauveau
%       Biorthogonal Bases of Compactly Supported Wavelete
%       Comm. Pure Appl., Math 45 (1992), p. 485 - 560
%       
% n and ntilde are integers, with ntilde >= 1, n + ntilde even

%  check arguments

if (ntilde < 1)
    error('cohen: ntilde must be >= 1');
end

if (~ iseven(n + ntilde))
    error('cohen: n + ntilde must be even');
end

if (iseven(n))
    k = 0;
else
    k = 1;
end

%  h coefficients
%  these are just binomials, i.e. phi is a B-spline

h = binomial(n,0:n) / 2^n;
h = reshape(h,1,1,length(h));
hmin = (k - n)/2;
h = mpoly(sym(h),hmin);

%  g coefficients

z = mpoly(sym(1),1);
term1 = ((1-z)/2)^ntilde;

L = (n + ntilde)/2;
PL = cell(1,L);
for i = 0:L-1
    PL{i+1} = binomial(L-1+i,i);
end
PL = mpoly(PL);

term2 = PL(mpoly({1/4,1/2,1/4},-1));
g = term1 * term2;

l = (2 - k - ntilde)/2;
g = (-1)^(l+1) * g;
g.min = g.min + l; 

H = [h;g];
H.type = 'symbol';
H.r = 1;
H.m = 2;

%  htilde, gtilde coefficients

P = polyphase(H);
Pt = longinv(P');
Ht = symbol(Pt);
    
function l = iseven(x)

% ISEVEN -- check whether an integer is even or odd
%  
%        l = iseven(x)
% 
% returns 1 = true   if x is an even integer,
%         0 = false  otherwise

if (x == 2 * round(x/2))
    l = 1;
else
    l = 0;
end