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

    function xrec = imodwt(w,varargin)
%MATLAB Code Generation Library Function

%   Copyright 1995-2016 The MathWorks, Inc.
%#codegen

narginchk(1,5);
validateattributes(w,{'double'},{'real','nonnan','finite'},mfilename);
coder.internal.assert(~isvector(w),'Wavelet:modwt:InvalidCFSSize');
ncw = coder.internal.indexInt(size(w,2));
N = coder.internal.indexInt(ncw);
Nrep = N;
J = coder.internal.indexInt(size(w,1)) - 1;
ZERO = coder.internal.indexInt(0);
ONE = coder.internal.indexInt(1);
%Parse input arguments
defaultLev = ZERO;
defaultWave = 'sym4';
[lev,Lo,Hi,reflection] = parseModwtInputs( ...
    defaultLev,defaultWave,varargin{:});
lev = lev + 1;
coder.internal.assert(lev <= J,'Wavelet:modwt:Incorrect_ReconLevel');
% Scale the scaling and wavelet filters for the MODWT
Lo = Lo./sqrt(2);
Hi = Hi./sqrt(2);
% Adjust final output length if MODWT obtained with 'reflection'
if reflection
    N = eml_rshift(N,ONE);
end
% If the number of samples is less than the length of the scaling filter
% we have to periodize the data and then truncate.
nLo = coder.internal.indexInt(numel(Lo));
if Nrep < nLo
    ncopies = nLo - Nrep;
    Nrep = (1 + ncopies)*N;
else
    ncopies = ZERO;
end
% Allocate storage for a possibly expanded data array.
ww = coder.nullcopy(zeros(size(w,1),Nrep));
% Copy one page of elements.
ww(:,1:ncw) = w;
if ncopies > 0
    for k = 1:ncopies
        offset = k*ncw;
        for j = 1:ncw
            ww(:,offset + j) = ww(:,j);
        end
    end
end
% Obtain the DFT of the filters
G = fft(Lo,double(Nrep));
H = fft(Hi,double(Nrep));
% IMODWT algorithm
vout = ww(J+1,:);
for jj = J:-1:lev
    vout = imodwtrec(vout,ww(jj,:),G,H,jj);
end
% Return proper output length
xrec = vout(1:N);

%----------------------------------------------------------------------

function V = imodwtrec(V,W,G,H,J)
N = coder.internal.indexInt(length(V));
ONE = coder.internal.indexInt(1);
upfactor = eml_lshift(ONE,J - 1); % 2^(J - 1);
Vhat = fft(V);
What = fft(W);
for k = 1:N
    idx = 1 + mod(upfactor*(k - 1),N);
    Gupk = conj(G(idx));
    Hupk = conj(H(idx));
    % Reuse storage for Vhat: Vhat = Gup.*Vhat + Hup.*What
    Vhat(k) = Gupk*Vhat(k) + Hupk*What(k);
end
V = real(ifft(Vhat));

%--------------------------------------------------------------------------