continuous_moment.m

几个关于多小波的程序

function Mu = continuous_moment(H,k)

% CONTINUOUS_MOMENT -- compute multiwavelet continuous moments
%
%        Mu = continuous_moment(H,k)
%
% The kth continuous moment of the multiscaling function is defined as
%
%        \mu_k  = \int x^k  \phi(x) dx
%
% and likewise for the moments \nu_k of the multiwavelets.
% k must be an integer >= 0 (default is k=0).
%
% Continuous moments are only defined up to some constant
% multiplier, which is fixed by choosing some particular \mu_0. We
% normalize \mu_0 to have 2-norm 1, with the largest entry positive.
%
% The matrix Mu contains all continuous moments up to order k in
% its columns (since we have to compute them, anyway).
%
% H can be a symbol or a polyphase matrix.

% 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)
    k = 0;
end

r = H.r;
m = H.m;
nwavelet = (size(H.coef,1)/r) - 1;

% compute discrete moments up to order k
M = discrete_moment(H,0:k);

% compute \mu_0
mu_0 = nullspace(M(1:r,:,1) - eye(r));
if (size(mu_0,2) > 1)
    mu_0
    error('mu_0 is not unique');
elseif (size(mu_0,2) == 0)
    disp('eigenvalues of H(0) are');
    eig(M(1:r,:,1));
    error('mu_0 not found');
end
mu_0 = normalize(mu_0);

% compute zeroth moment of wavelet functions
for i = 1:nwavelet
    mu_0(i*r+1:(i+1)*r) = M(i*r+1:(i+1)*r,:,1) * mu_0(1:r);
end

% compute the rest
Mu = zeros(size(H,1),k+1);
if (isa(H.coef,'sym'))
    Mu = sym(Mu);
end
Mu(:,1) = mu_0;
for j = 1:k
    % scaling function
    rhs = zeros(r,1);
    for s = 0:j-1
	rhs = rhs + binomial(j,s)*M(1:r,:,j-s+1)*Mu(1:r,s+1);
    end
    Mu(1:r,j+1) = (m^j * eye(r) - M(1:r,:,1)) \ rhs;
    for i = 1:nwavelet
    % wavelets
	for s = 0:j
	    Mu(i*r+1:(i+1)*r,j+1) = Mu(i*r+1:(i+1)*r,j+1) + ...
		binomial(j,s)*M(i*r+1:(i+1)*r,:,j-s+1)*Mu(1:r,s+1);
	end
	Mu(i*r+1:(i+1)*r,j+1) = Mu(i*r+1:(i+1)*r,j+1) / m^j;
    end
end
if (isa(H.coef,'sym'))
    Mu = simplify(Mu);
end


function y = normalize(x)
% normalize a vector to have length 1
% with the largest entry positive

n = size(x);
x = x(:);
nrm = sqrt(x'*x);
y = x / nrm;
[ymax,imax] = max(abs(y));
y = y / sign(y(imax));
y = reshape(y,n);