projection_shift.m

几个关于多小波的程序

function [Pnew,Ptnew,F] = projection_shift(P,Pt,nshift,side,Qsize)

% PROJECTION_SHIFT -- relative shift of biorthogonal multiwavelets
%
%        [Pnew,Ptnew,F] = projection_shift(P,Pt,nshift,side,Qsize)
%
% Given a biorthogonal pair of polyphase matrices P, Pt, we can use 
% a projection factor to shift the position of Pt one step to the left or right 
% while keeping the position of P invariant. This is useful for decomposing
% a biorthogonal multiwavelet into projection factors and an atom,
% or for constructing a multiwavelet from an atom.
%
% We do this NSHIFT times (to the right if NSHIFT is positive, to
% the left if it is negative). Default is -1.
%
% SIDE is the side on which the projection factors are
% applied. This is 'left' or 'right' (default is 'right').
% (A projection factor on the right can shift Pt either left
% or right; likewise for a factor on the left).
%
% We may have some freedom in the choice of the projection Q. If
% QSIZE=0 (the default), we make the range as small as
% possible. Anything else makes the range as large as possible.
%
% If SIDE = 'right', F is a cell array {F_n, F_{n-1}, ..., F_1} of
% projection factors so that 
%
%        P  = P_new   * F_n * F_{n-1} * .. * F_1
%
% If SIDE = 'left', F is a cell array {F_1, F_2, ..., F_n} of
% projection factors so that 
%
%        P  = F_1 * F_2 * .. * F_n * P_new
%        Pt = F_1 * F_2 * .. * F_n * Pt_new
%        Pt = Pt_new  * F_n * F_{n-1} * .. * F_1

% 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

% supply default arguments, do error checking

if (nargin < 5)
    Qsize = 0;
end

if (nargin < 4)
    side = 'right';
end
side = lower(side(1));

if (nargin < 3)
    nshift = -1;
end

if (nshift == 0)
    Pnew = P;
    Ptnew = Pt;
    F = cell(1,0);
    return
end

if (nshift > 0)
    if (nshift > P.max - Pt.min)
	error('shift not possible');
    end
else
    if (-nshift > Pt.max - P.min)
	error('shift not possible');
    end
end

P = polyphase(P);
Pt = polyphase(Pt);

% do the shifts

Pnew = P;
Ptnew = Pt;

direction = sign(nshift);
nshift = abs(nshift);
F = cell(1,nshift);

for i = 1:nshift
    [Pnew,Ptnew,factor] = shift1(Pnew,Ptnew,direction,side,Qsize);
    if (side == 'l')
	F{i} = factor;
    else
	F{nshift-i+1} = factor;
    end
end

function [Pnew,Ptnew,factor] = shift1(P,Pt,direction,side,Qsize);

% This is the routine that does the actual shift

I = eye(size(P));

switch direction
 case 1 % right shift
  switch side
   case 'l' % factor on left side
    switch Qsize
     case 0 % smallest possible range
      Q = range(Pt{Pt.min});
     otherwise % largest possible range
      Q = nullspace(P{P.max}');
    end
   otherwise % factor on right side
    switch Qsize
     case 0 % smallest possible range
      Q = range(Pt{Pt.min}');
     otherwise % largest possible range
      Q = nullspace(P{P.max});
    end
  end
  Q = Q * Q';
  factor = mpoly({I-Q,Q},0);
 otherwise % left shift
  switch side
   case 'l' % factor on left side
    switch Qsize
     case 0 % smallest possible range
      Q = range(Pt{Pt.max});
     otherwise % largest possible range
      Q = nullspace(P{P.min}');
    end
   otherwise % factor on right side
    switch Qsize
     case 0 % smallest possible range
      Q = range(Pt{Pt.max}');
     otherwise % largest possible range
      Q = nullspace(P{P.min});
    end
  end
  Q = Q * Q';
  factor = mpoly({Q,I-Q},-1);
end

switch side
 case 'l'
  Pnew = factor * P;
  Ptnew = factor * Pt;
 otherwise
  Pnew = P * factor;
  Ptnew = Pt * factor;
end
factor = factor';