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

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

narginchk(1,5);
validateattributes(wtmp,{'double','single'},{'nonnan','finite'},mfilename);
coder.internal.assert(~isvector(wtmp),'Wavelet:modwt:InvalidCFSSize');
isSingle = isUnderlyingType(wtmp,'single');
filterSingle = false;
ncw = coder.internal.indexInt(size(wtmp,2));
Nsig = coder.internal.indexInt(size(wtmp,3));
N = coder.internal.indexInt(ncw);
Nrep = N;
J = coder.internal.indexInt(size(wtmp,1)) - 1;
ZERO = coder.internal.indexInt(0);
ONE = coder.internal.indexInt(1);
%Parse input arguments
defaultLev = ZERO;
defaultWave = 'sym4';
[lev,Lotmp,Hitmp,reflection] = parseModwtInputs( ...
    defaultLev,defaultWave,varargin{:});
if isUnderlyingType(Lotmp,'single') || isUnderlyingType(Hitmp,'single')
    filterSingle = true;
end
lev = lev + 1;
coder.internal.assert(lev <= J,'Wavelet:modwt:Incorrect_ReconLevel');
% Scale the scaling and wavelet filters for the MODWT
LoDouble = double(Lotmp);
HiDouble = double(Hitmp);
LoScal = LoDouble./sqrt(2);
HiScal = HiDouble./sqrt(2);
if isSingle || filterSingle
    Lo = cast(LoScal,'single');
    Hi = cast(HiScal,'single');
else
    Lo = LoScal;
    Hi = HiScal;
end

if filterSingle && ~isSingle
    w = cast(wtmp,'single');
else
    w = wtmp;
end
% 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,Nsig,'like',w));
% 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 = squeeze(vout(:,1:N,:));

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

function V = imodwtrec(V,W,G,H,J)
N = int64(size(V,2));
ONE = int64(1);
upfactor = eml_lshift(ONE,J - 1); % 2^(J - 1);
isReal = isreal(V);
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,:) = bsxfun(@times,Gupk,Vhat(:,k,:)) + bsxfun(@times,Hupk,What(:,k,:));
end
if isReal
    V = real(ifft(Vhat));
else
    V = ifft(Vhat);
end

        

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