projection_factorization.m

几个关于多小波的程序

function [factors,factorst] = projection_factorization(P,Pt)

%  PROJECTION_FACTORIZATION -- factorize a poplyphase matrix
%
%        F = projection_factorization(P)
%        [F,Ft] = projection_factorization(P,Pt)
%
%  The first form is for orthogonal wavelets, the second
%  form is for biorthogonal wavelets
%
% In the orthogonal case, F is a cell array of the form
%            {Q, F_1, F_2, ..., F_n, z^k}
%
% with
%
%        P(z) = Q * F_1(z) * ... * F_n(z) * z^k
%
% where Q is a constant orthogonal matrix and the F_k are projection factors.
% The term z^k is omitted if k=0.
%
% In the biorthogonal case, F is a cell array of the form
%
%            {Q, A, F_1, ..., F_n, z^k}
%
%  with
%
%        P(z) = Q * A(z) * F_1(z) * ... * F_n(z) * z^k
%
% and likewise for Pt and Ft. Here Q is a constant matrix, 
% A is an atom, and the F_k are projection factors.
% The term z^k is omitted if k=0.

% 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

P = polyphase(P);
r = P.r;
m = P.m;
z = mpoly(1,1);
I = eye(m*r);

if (nargin == 1) % orthogonal case

% error check
    if (~isidentity(P*P'))
	error('this is not an orthogonal polyphase matrix');
    end

% pull out power of z to make P start at 0
    if (P.min ~= 0)
	factors = {mpoly(1,P.min)};
	P.min = 0;
    else
	factors = {};
    end

% factor
    while (P.length > 1)
	U = range(P{end}');
	F = projection_factor(U,m,r);
	Finv = reverse(F);
	factors = {F, factors{:}};
% prevent infinite loops
	oldlength = P.length;
	P = P * Finv;
	if (P.length >= oldlength)
	    error('P is not getting shorter; this should not happen');
	end
    end
% what is left is a constant matrix
    factors = {P.coef, factors{:}};

else % biorthogonal case
    
% error check
    Pt = polyphase(Pt);
    if (~isidentity(P*Pt'))
	error('this is not a biorthogonal polyphase matrix pair');
    end

% pull out power of z to make P start at 0
    if (P.min ~= 0)
	factors  = {mpoly(1,P.min)};
	factorst = {mpoly(1,P.min)};
	Pt.min = Pt.min - P.min;
	P.min = 0;
    else
	factors  = {};
	factorst = {};
    end
    
% factor
   while (Pt.max > P.min)
       oldPtmax = Pt.max;
       [P,Pt,F] = projection_shift(P,Pt);
       factors = {F{1}, factors{:}};
       factorst = {F{1}.', factorst{:}};
% prevent infinite loops
	if (Pt.max >= oldPtmax | P.min ~= 0)
	    error('Pt is not shifting correctly; this should not happen');
	end
   end
% we have reached either an atom, or the last projection factor
   Q = moment(P);
   P = Q\P;
   Pt = Q'*Pt;
   factors = {Q,P,factors{:}};
   factorst = {inv(Q)',Pt,factorst{:}};
end