www.gusucode.com > wavelet工具箱matlab源码程序 > wavelet/wavelet/imodwt.m

    function xrec = imodwt(w,varargin)
%IMODWT Inverse maximal overlap discrete wavelet transform.
%   XREC = IMODWT(W) returns the reconstructed signal based on the maximal
%   overlap discrete wavelet transform (MODWT) coefficients in W. W is a
%   LEV+1-by-N matrix which is the MODWT of an N-point input signal down to
%   level LEV. By default, IMODWT assumes that you used the 'sym4' wavelet
%   with periodic boundary handling to obtain the MODWT. If you do not
%   modify the coefficients, XREC is a perfect reconstruction of the
%   signal.
%
%   XREC = IMODWT(W,WNAME) reconstructs the signal using the wavelet WNAME.
%   WNAME must be the same wavelet used in the analysis of the signal with
%   MODWT.
%
%   XREC = IMODWT(W,Lo,Hi) reconstructs the signal using the scaling
%   filter Lo and the wavelet filter Hi. You cannot specify both WNAME and
%   Lo and Hi. Lo and Hi must be the same filters used in the analysis with
%   MODWT.
%
%   XREC = IMODWT(...,LEV) reconstructs the signal up to level LEV. XREC is
%   a projection onto the scaling space at level LEV. LEV is a nonnegative
%   integer between 0 and strictly less than size(W,1)-1. The default is
%   LEV = 0, which results in perfect reconstruction if you do not modify
%   the coefficients. 
%
%   XREC = IMODWT(...,'reflection') uses the 'reflection' boundary
%   condition in the reconstruction. If you specify 'reflection', IMODWT
%   assumes that the length of the original signal is 1/2 the number of
%   columns in the input coefficient matrix W. You must enter the entire
%   string 'reflection'. If you added a wavelet named 'reflection' using
%   the wavelet manager, you must rename that wavelet prior to using this
%   option. 'reflection' may be placed in any position in the input
%   argument list after W. By default both MODWT and IMODWT assume periodic
%   signal extension at the boundary.
%
%   %Example 1:
%   %   Demonstrate perfect reconstruction of ECG data with the MODWT.
%
%   load wecg;
%   wecg = wecg(1:end-1);
%   w = modwt(wecg,'sym4',10);
%   xrec = imodwt(w);
%   max(abs(xrec-wecg'))
%   subplot(2,1,1)
%   plot(wecg); title('Original Signal');
%   subplot(2,1,2)
%   plot(xrec); title('Reconstructed Signal');
%
%   %EXAMPLE 2:
%   %   Reconstruct a signal approximation based on the level-3 and
%   %   level-4 wavelet coefficients.
%
%   load wecg;
%   wecg = wecg(1:end-1);
%   w = modwt(wecg,'db2',10);
%   idx = 3:4;
%   wnew = zeros(size(w));
%   wnew = w(idx,:);
%   xrec = imodwt(wnew);
%   subplot(211)
%   plot(wecg); title('Original Signal');
%   subplot(212)
%   plot(xrec);
%   title('Reconstruction from Level-3 and Level-4 Wavelet Coefficients');
%
%   See also modwt, modwtmra, modwtvar, modwtcorr, modwtxcorr

% Check number of input arguments
narginchk(1,5);

% Input cannot be a row or column vector. IMODWT expects at least a two row
% matrix
if (isrow(w) || iscolumn(w))
    error(message('Wavelet:modwt:InvalidCFSSize'));
end

% Input must be real-value, finite, and double precision
validateattributes(w,{'double'},{'real','nonnan','finite'});

% Parse input arguments
params = parseinputs(varargin{:});

% Get the original input size
% Get the level of the MODWT
N = size(w,2);
Nrep = N;
J = size(w,1)-1;


boundary = params.boundary;
if (~isempty(boundary) && ~strcmpi(boundary,'reflection'))
    error(message('Wavelet:modwt:Invalid_Boundary'));
end

% Adjust final output length if MODWT obtained with 'reflection'
if strcmpi(boundary,'reflection')
    N = N/2;
end

% If the wavelet is specified as a string, obtain filters from wavemngr
if (isfield(params,'wname') && ~isfield(params,'Lo'))
    [~,~,Lo,Hi] = wfilters(params.wname);
    wtype = wavemngr('type',params.wname);
    if (wtype ~= 1)
        error(message('Wavelet:modwt:Orth_Filt'));
    end
end

%If scaling and wavelet filters specified as vectors, ensure they
%satisfy the orthogonality conditions

if (isfield(params,'Lo') && ~isfield(params,'wname'))
    filtercheck = CheckFilter(params.Lo,params.Hi);
    if ~filtercheck
        error(message('Wavelet:modwt:Orth_Filt'));
    end
    Lo = params.Lo;
    Hi = params.Hi;
end


% Scale scaling and wavelet filters for MODWT
Lo = Lo./sqrt(2);
Hi = Hi./sqrt(2);


% If the number of samples is less than the length of the scaling filter
% we have to periodize the data and then truncate.
if (Nrep<numel(Lo))
    w = [w repmat(w,1,ceil(numel(Lo)-Nrep))];
    Nrep = size(w,2);
end

% Get the DFTs of the scaling and wavelet filters
G = fft(Lo,Nrep);
H = fft(Hi,Nrep);

% Get the level of the reconstruction
lev = params.lev+1;

% Error
if (lev>J)
    error(message('Wavelet:modwt:Incorrect_ReconLevel'));
end

vin = w(J+1,:);

% IMODWT algorithm

for jj = J:-1:lev
    vout = imodwtrec(vin,w(jj,:),G,H,jj);
    vin = vout;
end

% Return proper output length
xrec = vout(1:N);




%-----------------------------------------------------------------
function Vout = imodwtrec(Vin,Win,G,H,J)
N = length(Vin);
Vhat = fft(Vin);
What = fft(Win);
upfactor = 2^(J-1);
Gup = conj(G(1+mod(upfactor*(0:N-1),N)));
Hup = conj(H(1+mod(upfactor*(0:N-1),N)));
Vout = ifft(Gup.*Vhat+Hup.*What);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function params = parseinputs(varargin)
% Parse varargin and check for valid inputs

% Assign default inputs
params.boundary = [];
params.lev = 0;
params.wname = 'sym4';


% Check for 'reflection' boundary
tfbound = strcmpi(varargin,'reflection');

% Determine if 'reflection' boundary is specified
if any(tfbound)
    params.boundary = varargin{tfbound>0};
    varargin(tfbound>0) = [];
end

% If boundary is the only input in addition to the data, return with
% defaults
if isempty(varargin)
    return;
end

% Only remaining char variable must be wavelet name
tfchar = cellfun(@ischar,varargin);

if (nnz(tfchar) == 1)
    params.wname = varargin{tfchar>0};
end

% Only scalar input must be the level
tfscalar = cellfun(@isscalar,varargin);

% Check for numeric inputs
tffilters = cellfun(@isnumeric,varargin);

% At most 3 numeric inputs are supported
if nnz(tffilters)>3
    error(message('Wavelet:modwt:Invalid_Numeric'));
end

% If there are at least two numeric inputs, the first two must be the
% scaling and wavelet filters

if (nnz(tffilters)>1)
    idxFilt = find(tffilters,2,'first');
    params.Lo = varargin{idxFilt(1)};
    params.Hi = varargin{idxFilt(2)};
    if (length(params.Lo) < 2 || length(params.Hi) < 2)
        error(message('Wavelet:modwt:Invalid_Filt_Length'));
    end
    params = rmfield(params,'wname');
end


% Any scalar input must be the level
if any(tfscalar)
    params.lev = varargin{tfscalar>0};
end

% If the user specifies a filter, use that instead of default wavelet
if (isfield(params,'Lo') && any(tfchar))
    error(message('Wavelet:FunctionInput:InvalidWavFilter'));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = CheckFilter(Lo,Hi)
% For a user-supplied scaling and wavelet filter, check that
% both correspond to an orthogonal wavelet
Lo = Lo(:);
Hi = Hi(:);
Lscaling = length(Lo);
Lwavelet = length(Hi);
evenlengthLo = 1-rem(Lscaling,2);
evenlengthHi = 1-rem(Lwavelet,2);

if all([evenlengthLo evenlengthHi])
    evenlength = 1;
else
    evenlength = 0;
end

if (Lscaling ~= Lwavelet)
    equallen = 0;
else
    equallen = 1;
end

normLo = norm(Lo,2);
sumLo = sum(Lo);
normHi = norm(Hi,2);
sumHi = sum(Hi);
tol = 1e-7;

if (abs(normLo-1) > tol && abs(normHi -1) > tol)
    unitnorm = 0;
else
    unitnorm = 1;
end

if (abs(sumLo-sqrt(2))> tol && abs(sumHi) > tol)
    sumfilters = 0;
else
    sumfilters = 1;
end


zeroevenlags = 1;

if (Lscaling > 2)
    L = Lscaling;
    xcorrHi = conv(Hi,flipud(Hi));
    xcorrLo = conv(Lo,flipud(Lo));
    xcorrLo = xcorrLo(L+2:2:end);
    xcorrHi = xcorrHi(L+2:2:end);
    zeroevenlagsLo = 1-any(abs(xcorrLo>tol));
    zeroevenlagsHi = 1-any(abs(xcorrHi>tol));
    zeroevenlags = max(zeroevenlagsLo,zeroevenlagsHi);
end
out = all([evenlength equallen unitnorm sumfilters zeroevenlags]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%










% [EOF] imodwt.m