correlation.m

几个关于多小波的程序

function H = correlation(H1,H2,type)

% CORRELATION -- find symbol for correlation of multiwavelets H1, H2
%
%        H = correlation(H1,H2,'type')
%
% Assume \phi_1, \phi_2 are two refinable vectors
%
%        \phi_j(x) = \sqrt(m) \sum_k h_{jk} \phi_j (mx-k),    j = 1,2
%
% Their correlation is the function
%
%        C(x) = \int \phi_1(y) \phi_2(y+x)^*  dy
% 
% This is a matrix-valued refinable function. When stretched out into
% a vector c = vec(C), it satisfies a recursion relation of the form
%
%        c(x) = \sqrt(m) \sum_k h_k  c(mx-k),
%
% where
%
%        h_k  =  1/\sqrt(m) \sum_j  \overline{h_{2,j+k}} \otimes h_{1,j}
%
% with \otimes = Kronecker product. (That is a complex conjugate on h_2, 
% but not a complex conjugate transpose).
%
% The optional third argument is a string of two characters, each '+' or '-'.
% If it is present, calculate
%
%        C(x) = \int phi_1(y) phi_2(+-y +- x)^*  dy
% 
% instead, which leads to
%
%        h_k =  1/\sqrt(m) \sum_j \overline{h_{2,+-j+-k}} \otimes h_{1,j}
%
%  The first +- corresponds to y and j, the second +- corresponds to x and k.
%
% If H2 is omitted, use H2 = H1. Default for TYPE is '+-'. For H2 = H1,
% that is the autocorrelation.

% 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)
    H2 = H1;
    type = '+-';
end

if (nargin == 2)
    if isa(H2,'char')
	type = H2;
	H2 = H1;
    else
	type = '+-';
    end
end

if ((type(1) ~= '-' & type(1) ~= '+') | (type(2) ~= '-' & type(2) ~= '+'))
    error('TYPE must be two plus or minus signs');
end

% reduce H1, H2 to symbols of scaling functions alone

H1 = symbol(H1);
r1 = H1.r;
m1 = H1.m;
H1 = trim(H1(1:r1,:));
h1min = H1.min;
h1max = H1.max;
h1len = H1.length;

H2 = symbol(H2);
r2 = H2.r;
m2 = H2.m;
H2 = trim(H2(1:r2,:));
h2min = H2.min;
h2max = H2.max;
h2len = H2.length;

if (m1 ~= m2)
    error('incompatible dilation factors');
end


llen = h2len + h1len - 1;
Hcoef = zeros(r1*r2,r1*r2,llen);
if (isa(H1.coef,'sym') & isa(H2.coef,'sym'))
    Hcoef = sym(Hcoef);
end

if (type(1) == '+')
    lmin = h2min - h1max;
    lmax = h2max - h1min;
    for l = lmin:lmax
	for j = h1min:h1max
	    Hcoef(:,:,l-lmin+1) = Hcoef(:,:,l-lmin+1) + kron(conj(H2{j+l}),H1{j});
	end
    end
    H = mpoly(Hcoef,lmin,'symbol',m1,r1*r2);
else
    lmin = h2min + h1min;
    lmax = h2max + h1max;
    for l = lmin:lmax
	for j = h1min:h1max
	    Hcoef(:,:,l-lmin+1) = Hcoef(:,:,l-lmin+1) + kron(conj(H2{-j+l}),H1{j});
	end
    end
    H = mpoly(Hcoef,lmin,'symbol',m1,r1*r2);
end

% reverse order if type(2) is '-'
if (type(2) == '-')
% H = mpoly(H,coef(:,:,llen:-1:1),-lmax,'symbol',m1,r1*r2);
    H = (H').';
end

% scaling is not needed: the elements of the symbols H1, H2 are
% actually h1_j/\sqrt(m), h2_j/\sqrt(m), and since H is again a symbol,
% the elements are actually H_j/\sqrt(m). When we consider that, plus
% the factor 1/sqrt(m) from the definition of H_j, all the factors
% cancel.