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