sobolev_exponent.m

几个关于多小波的程序

function s = sobolev_exponent(H)

% SOBOLEV_EXPONENT - estimate Sobolev exponent of multiwavelet (slower, but better)
%
%         s = sobolev_exponent(H)
%
% This routine implements the algorithm given in the paper
%
%        Qingtang Jiang, On the Regularity of Matrix Refinable Functions
%        SIAM J. Math. Anal 29 # 5, Sept. 1998, pp. 1157-1176
%
% The multiscaling function phi given by the coefficients H is in the Sobolev
% space W^\sigma for any \sigma < s = - log_4(rho(T|H^0))
% (the log of the spectral radius of the transition operator T
% restricted to a certain subspace H^0).
%
% This is a better estimate than the one provided by routine SOBOLEV_ESTIMATE,
% but it runs slower. (Try multiwavelet CL2 for an example of
% where the two differ).

% Copyright (c) 2004 by Fritz Keinert ([e-mail]),
% Dept. of Mathematics, Iowa State University, Ames, IA 50011.
% This software may be freely used and distributed for non-commercial
% purposes, provided this copyright statement is preserved, and
% appropriate credit for its use is given.
%
% Last update: Feb 20, 2004

global MPOLY_TOLERANCE

if isempty(MPOLY_TOLERANCE)
    tol = 1.e-12;
else
    tol = MPOLY_TOLERANCE;
end

T = transition_matrix(H);
[L,l,r] = findLlr(H);
if (size(L,1) > 0)
    V = [L;l;r]';  % H^0 = orthogonal complement of columns of V
    [Q,R,P] = qr(V);
    r = rank(V);
    TQ = Q'*T*Q;
    n = size(T,1);
    discrepancy = max(max(abs(TQ(1:r,r+1:n))));
    if (discrepancy > tol)
	disp('T does not map H^0 into H^0');
    end
    Tr = TQ(r+1:n,r+1:n);
else
    Tr = T;
end
lambda = eig(Tr);
rho = max(abs(lambda));
s = -log(rho)/log(4);

function [L,l,r] = findLlr(H)

% FINDLLR - find left eigenvectors of transition matrix
%
%        [L,l,r] = findLlr(H)
%
% Given multiwavelet coefficients H, find vectors L, l, r which are
% left eigenvectors of the transition matrix, from the paper
%
%        Qingtang Jian, On the Regularity of Matrix Refinable Equations,
%        SIAM J. Math. Anal. vol. 29 #5, Sept. 1998, pp. 1157-1176.
%
% On output, L, l and r contain the vectors from Jiang's paper in their rows.

% Note: in this program, m is the multiplicity (usually called r),
% because the letter r is needed for the right eigenvectors.

H = symbol(H);
rr = H.r;
H = H(1:rr,:);
Hmin = H.min;
Hmax = H.max;
Hlen = Hmax - Hmin + 1;
N = Hlen - 1;
p = approximation_order(H);
[p,Y] = approximation_order(H,2*p);
if (p == 0)
    L = zeros(0,rr^2*(2*N-1));
    l = zeros(0,rr^2*(2*N-1));
    r = zeros(0,rr^2*(2*N-1));
    return
end

% l0 contains the y^(l)_0 = l^j_0 in Jiang's notation
% this is just Y reordered for easier access

l0 = cell(1,2*p);
for i = 0:2*p-1
    l0{i+1} = Y(:,i+1)';
end

% l1 contains the l^i_k from Jiang

l1 = cell(2*p,2*N+1);
for j = 0:2*p-1
    for k = -N:N
	l1{j+1,k+N+1} = zeros(1,rr);
	for i = 0:j
	    l1{j+1,k+N+1} = l1{j+1,k+N+1} + binomial(j,i) * k^(j-i) * l0{i+1};
	end
    end
end

%                  i
% l2 contains the l (k) from Jiang

l2 = cell(2*p,2*N+1);
for j = 0:2*p-1
    for k = -N:N
	l2{j+1,k+N+1} = zeros(1,rr^2);
	for i = 0:j
	    l2{j+1,k+N+1} = l2{j+1,k+N+1} + (-1)^i * binomial(j,i) * kron(l1{i+1,k+N+1},l0{j-i+1});
	end
    end
end

%                 j
% L contains the L  from Jiang

L = cell(2*p,1);
for j = 0:2*p-1
    L{j+1} = zeros(2*N+1,rr^2);
    for k = -N:N
	L{j+1}(k+N+1,:) = l2{j+1,k+N+1};
    end
    temp = L{j+1}';
    L{j+1} = temp(:)';
end

%                    i    i
% l, r contain the  l ,  r  from Jiang
%                  j    j

l = cell(p,rr);
r = cell(p,rr);
I = eye(rr);
for i = 0:p-1
    for j = 1:rr
	l{i+1,j} = zeros(2*N+1,rr^2);
	r{i+1,j} = zeros(2*N+1,rr^2);
	for k = -N:N
	    l{i+1,j}(k+N+1,:) = kron(I(j,:),l1{i+1,-k+N+1});
	    r{i+1,j}(k+N+1,:) = kron(l1{i+1,k+N+1},I(j,:));
	end
	temp = l{i+1,j}';
	l{i+1,j} = temp(:)';	
	temp = r{i+1,j}';
	r{i+1,j} = temp(:)';
    end
end

% convert L, l, r to regular arrays

L = cat(1,L{:});
l = cat(1,l{:});
r = cat(1,r{:});