idwt.m

几个关于多小波的程序

function s = idwt(s,Ht,nlevel)

%  IDWT -- Inverse Discrete Wavelet Transform
%
%        s = idwt(sd,Ht,nlevel)
%
%  Performs a discrete wavelet transform of the signal S
%  through NLEVEL levels, using the wavelet described by Ht.
%  Ht 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 input 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 = Ht.m;
r = Ht.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 reconstruction steps');
end

% check that Ht contains scaling function and all wavelets
[n1,n2] = size(Ht);
if (n1 ~= m*r | n2 ~= r)
    error('size of Ht not compatible with M and R');
end

% convert Ht into coefficients alone
switch Ht.type
 case 'symbol'
  Ht.coef = Ht.coef * sqrt(m);
  Ht.type = '';
  
 case 'polyphase'
  Ht = symbol(Ht);
  Ht.coef = Ht.coef * sqrt(m);
  Ht.type = '';
end

% build the filterbank
newmin = floor(Ht.min/m)*m;
newmax = ceil((Ht.max+1)/m)*m-1;
Ht = trim(Ht,[newmin,newmax]);
Ft = mpoly(reshape(Ht.coef,m*r,m*r,prod(size(Ht.coef))/(m^2*r^2)), newmin/m, '', m, m*r);
Ft = Ft';

% reconstruct
ls = length(s)/m^(nlevel-1);
for i = 1:nlevel
    s(1:ls) = idwt_one_level(s(1:ls),Ft);
    ls = ls*m;
end

function s = idwt_one_level(sd,Ft)

% sort input from form (s0,s1,...,d0,d1,...) into form (s0,d0,s1,d1,...)
m = Ft.m;
mr = Ft.r;
r = mr/m;
n = prod(size(sd))/mr;
sd = reshape(sd,r,prod(size(sd))/r);

temp = zeros(mr,n);
if (isa(sd,'sym'))
    temp = sym(temp);
end

for i = 0:m-1
    temp(r*i+1:r*(i+1),:) = sd(:,n*i+1:n*(i+1));
end

% pad or cut at the beginning and/or end
if (Ft.min < 0)
    head1 = n+Ft.min+1;
    head2 = n;
    middle1 = 1;
else
    head1 = 1;
    head2 = 0;
    middle1 = Ft.min + 1;
end

if (Ft.max > 0)
    middle2 = n;
    tail1 = 1;
    tail2 = Ft.max;
else
    middle2 = n + Ft.max;
    tail1 = 1;
    tail2 = 0;
end

temp = [temp(:,head1:head2),temp(:,middle1:middle2),temp(:,tail1:tail2)];

% do filtering
s = zeros(mr,n);
for i = 1:Ft.length
    s = s + Ft.coef(:,:,i) * temp(:,i:i+n-1);
end

s = s(:);