approximation_order.m

几个关于多小波的程序

function [p,Y] = approximation_order(H,n)

% APPROXIMATION_ORDER -- compute approximation order and approximation vectors
%
%        [p,Y] = approximation_order(H)
%
% This routine determines the approximation order of a refinable 
% function vector, as well as the approximation vectors.
% The vectors are returned in the columns of a matrix Y.
%
% H can be a symbol (of the multiscaling function alone, or of the
% entire multiwavelet) or a polyphase matrix.

% Algorithm:
%
% This routine is based on the sum rules in terms of partial moments.
% Instead of solving for y_0, y_1, ... in turn, we solve for all
% y_k, k=0,...,p-1 simultaneously, until it doesn't work any more
% for some p.
%
% So, when we calculate y_2, we recompute y_0 and y_1. What's the point?
% In numerical experiments years ago I ran into cases where you can
% get multiple solutions for y_1, for example, but when you look
% for y_2, only one of the y_1 works.
%
% I never completely figured out why this happens. Most likely, these
% scaling functions are not stable. Anyway, this routine will continue
% when multiple solutions are found, and hope they will go away at
% a later point. It does print warning messages when this happens.

% 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

% convert H into the symbol of the refinable function vector alone
H = symbol(H);
r = H.r;
m = H.m;
H = H(1:r,:);
H = trim(H);

% check for approximation order 1
M = mpoly(partial_moment(H,0)');
I = kron(ones(m,1),eye(r));

A = M{0} - I/m;
y = nullspace(A);

if (size(y,2) > 1)
    disp('multiple y_0 found');
end

if (size(y,2) == 0)
    p = 0;
    Y = zeros(r,0);
    return
end


% check for higher orders

for p = 1:1000 
    M{p} = partial_moment(H,p)';
    for s = 0:p
	A(p*r*m+1:(p+1)*r*m,s*r+1:(s+1)*r) = binomial(p,s)*(-1)^(p-s)*m^s*M{p-s};
    end
    A(p*r*m+1:(p+1)*r*m,p*r+1:(p+1)*r) = A(p*r*m+1:(p+1)*r*m,p*r+1:(p+1)*r) - I/m;
    ynew = nullspace(A);

    if (size(ynew,2) > 1)
	disp(['multiple y_',num2str(p),' found']);
    end

    if (size(ynew,2) == 0)
	break
    end
    
    y = ynew;
end

% try to find more vectors if requested

if (nargin < 2)
    n = p;
end

if (n > p)
    M = mpoly(moment(H,0)');
    for j = 1:p-1
	Mp = moment(H,j);
	M{j} = Mp';
    end
    A = A(1:p*r*m,1:p*r);
end

for j = p:n-1
    Mp = moment(H,j);
    M{j} = Mp';
    for s = 0:j
	B((j-p)*r+1:(j-p+1)*r,s*r+1:(s+1)*r) = binomial(j,s) * (-1)^(j-s)*m^s*M{j-s};
    end
    B((j-p)*r+1:(j-p+1)*r,j*r+1:(j+1)*r) = B((j-p)*r+1:(j-p+1)*r,j*r+1:(j+1)*r) - eye(r);
    ynew = nullspace([A,zeros(p*m*r,(j-p+1)*r);B]);
    if (size(ynew,2) > 1)
	disp(['multiple y_',num2str(j),' found']);
    end

    if (size(ynew,2) == 0)
	break
    end

    y = ynew;
end    

% normalize y_0 so that (y_0^*)(mu_0) = 1;
Y = reshape(y,r,prod(size(y))/r);
mu_0 = continuous_moment(H,0);
Y = Y / (Y(:,1)' * mu_0);

function M = partial_moment(H,n)

% compute the nth partial moments of H, which are
% defined as
%                         _
%                   1    \          n
%        M     = ------  /_  (mk+nu)  h     . 
%         n,nu   sqrt(m)  k            mk+nu      
%
% n must be an integer >= 0 (default is n=0).
%
% The result is returned as a single matrix
%
% [M   , M   , ..., M     ]
%   n,0   n,1        n,m-1


if (nargin < 2)
    n = 0;
end

m = H.m;
[n1,n2] =size(H);
M = zeros(n1,m*n2);
if isa(H.coef,'sym')
    M = sym(M);
end
for k = H.min:H.max
    nu = mod(k,m);
    M(:,nu*n2+1:(nu+1)*n2) = M(:,nu*n2+1:(nu+1)*n2) + k^n * H{k};
end

% we don't need to divide by sqrt(m) because the symbol H
% already has the coefficients divided by sqrt(m)