itst.m

几个关于多小波的程序

function [Hn,Htn] = itst(H,Ht,C)

% ITST -- perform inverse two-scale similarity transform
%
%        Hn = itst(H,C)
%        [Hn,Htn] = itst(H,Ht,C)
%
% Perform one step of inverse two-scale similarity transform on a 
% multiwavelet symbol H, or on a pair of dual symbols H, Ht.
%
% The argument C is either
%
%  -  a TST matrix; C must be a linear matrix polynomial
%
%  -  the string 'lr', in which case we use
%
%        C(z) = I - (r l^*) / (l^* r) z
%
%     where l, r are the left and right eigenvectors of H(0) to
%     eigenvalue 1.
%
%  -  the string 'll', in which case we use r = l as above.

% 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

if (nargin == 2)
% arguments are H, C
    C = Ht;
end

% make sure H, Ht are symbols
H  = symbol(H);
if (nargin > 2)
    Ht = symbol(Ht);
end
r = H.r;
m = H.m;

M0 = moment(H,0);
M0 = M0(1:r,:);
z = mpoly(1,1);

% find TST matrix C
if isa(C,'char')
    
% find r, l
    R = null(M0 - eye(r));
    if (size(R,2) == 0)
	error('r not found');
    elseif (size(R,2) > 1)
	error('r is not unique')
    end
    
    L = null(M0' - eye(r));
    if (size(L,2) == 0)
	error('l not found');
    elseif (size(L,2) > 1)
	error('l is not unique')
    end
    
    switch C
     case 'lr'
      A = (R * L') / (L' * R);

     case 'll'
      A = (L * L') / (L' * L);
      
     otherwise
      error('third argument must be ''lr'' or ''ll''');
    end

    C = eye(r) - A * z;

elseif (isa(C,'mpoly'))
    if (C.min ~= 0 | C.max ~= 1)
	error('C must be a linear matrix polynomial');
    end
else
    error('C must be a string (''lr'' or ''rr'') or a linear matrix polynomial');
end

% decompose C = A + B(1-z)

B = - C{1};
A = C{0} - B;
DE = inv([A,B;B,A]);
n = size(C,1);
D = DE(1:n,1:n);
E = DE(n+1:2*n,1:n);
C0 = E + D*(1-z);

Cm = C;
Cm{m} = C{1};
Cm{1} = 0;
Cm0 = C0;
Cm0{m} = C0{1};
Cm0{1} = 0;

% apply inverse TST
[n1,n2] = size(H);
if (n1 == r) % H = scaling function only
    Hn = m * ((1-z^m) \ (Cm0 * H * C));
    
else
    Hn = m * ([(1-z^m) \ (Cm0 * H(1:r,:) * C);H(r+1:n1,:) * C]);
end

% apply TST to dual
if (nargin > 2)
    if (n1 == r) % Ht = scaling function only
	Htn = (1/m) * (Cm' * Ht * C0') / (1-z^(-1));
    else
	Htn = (1/m) * ([Cm' * Ht(1:r,:);Ht(r+1:n1,:)] * C0' )/ (1-z^(-1));
    end
end