几个关于多小波的程序

				function s = dwt(s,H,nlevel)
				
				%  DWT -- Discrete Wavelet Transform
				%
				%        sd = dwt(s,H,nlevel)
				%
				%  Performs a discrete wavelet transform of the signal S
				%  through NLEVEL levels, using the wavelet described by H.
				%  H can be an untyped matrix polynomial, in which case we 
				%  use the coefficients the way they are, or a symbol or 
				%  polyphase matrix. The default for NLEVEL is 1.
				%
				%  The output SD contains the S and D coefficients, in the form
				%
				%          (s , d , d   , ..., d   )
				%            n   n   n+1        m-1
				%
				%  where M = finest level, N = M - NLEVEL.
				
				% Copyright (c) 2004 by Fritz Keinert (keinert@iastate.edu),
				% 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
				
				m = H.m;
				r = H.r;
				s = s(:);
				
				if (nargin < 3)
				    nlevel = 1;
				end
				
				% check that M^NLEVEL*R divides the length of S
				if (round(length(s)/(m^nlevel*r))*m^nlevel*r ~= length(s))
				    error('length of S is incompatible with number of decomposition steps');
				end
				
				% check that H contains scaling function and all wavelets
				[n1,n2] = size(H);
				if (n1 ~= m*r | n2 ~= r)
				    error('size of H not compatible with M and R');
				end
				
				% convert H into coefficients alone
				switch H.type
				 
				 case 'symbol'
				  H.coef = H.coef * sqrt(m);
				  H.type = '';
				  
				 case 'polyphase'
				  H = symbol(H);
				  H.coef = H.coef * sqrt(m);
				  H.type = '';
				
				end
				
				% build the filterbank
				newmin = floor(H.min/m)*m;
				newmax = ceil((H.max+1)/m)*m-1;
				H = trim(H,[newmin,newmax]);
				F = mpoly(reshape(H.coef,m*r,m*r,prod(size(H.coef))/(m^2*r^2)), newmin/m, '', m, m*r);
				
				% decompose
				ls = length(s);
				for i = 1:nlevel
				    s(1:ls) = dwt_one_level(s(1:ls),F);
				    ls = ls/m;
				end
				
				function sd = dwt_one_level(s,F)
				
				% divide S into segments of length mr
				m = F.m;
				mr = F.r;
				r = mr/m;
				n = prod(size(s))/mr;
				s = reshape(s,mr,prod(size(s))/mr);
				
				% pad or cut at the beginning and/or end
				if (F.min < 0)
				    head1 = n+F.min+1;
				    head2 = n;
				    middle1 = 1;
				else
				    head1 = 1;
				    head2 = 0;
				    middle1 = F.min + 1;
				end
				
				if (F.max > 0)
				    middle2 = n;
				    tail1 = 1;
				    tail2 = F.max;
				else
				    middle2 = n + F.max;
				    tail1 = 1;
				    tail2 = 0;
				end
				
				s = [s(:,head1:head2),s(:,middle1:middle2),s(:,tail1:tail2)];
				
				% do filtering
				temp = zeros(mr,n);
				for i = 1:F.length
				    temp = temp + F.coef(:,:,i) * s(:,i:i+n-1);
				end
				
				% sort output from form (s0,d0,s1,d1,...) into form (s0,s1,...,d0,d1,...)
				sd = zeros(r,m*n);
				if (isa(temp,'sym'))
				    sd = sym(sd);
				end
				
				for i = 0:m-1
				    sd(:,n*i+1:n*(i+1)) = temp(r*i+1:r*(i+1),:);
				end
				sd = sd(:);