dwt.m

几个关于多小波的程序

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 ([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

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(:);