function [w,G,H] = modwt(x,varargin)
%MATLAB Code Generation Library Function

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

narginchk(1,5);
validateattributes(x,{'double','single'},{'nonnan','finite','2d'},mfilename);
isSingle = isUnderlyingType(x,'single');
isReal = isreal(x);
% Create column vector if input is row vector
if isrow(x)
    xtmp = x(:);
else
    xtmp = x;
end
datalength = coder.internal.indexInt(size(xtmp,1));
coder.internal.assert(datalength >= 2,'Wavelet:modwt:LenTwo');
%Parse input arguments
defaultLev = coder.internal.indexInt(floor(log2(double(datalength))));
defaultWave = 'sym4';
[lev,Lotmp,Hitmp,reflection] = parseModwtInputs( ...
    defaultLev,defaultWave,varargin{:});

Nsig = size(xtmp,2);
%Check that the level of the transform does not exceed floor(log2(numel(x))
maxlev = floor(log2(double(datalength)));
coder.internal.assert(lev > 0 && lev <= maxlev,'Wavelet:modwt:MRALevel');
% Scale the scaling and wavelet filters for the MODWT
if isSingle
    Lo = cast(Lotmp./sqrt(2),'single');
    Hi = cast(Hitmp./sqrt(2),'single');
else
    Lo = Lotmp./sqrt(2);
    Hi = Hitmp./sqrt(2);
end

% increase signal length if 'reflection' is specified
nLo = coder.internal.indexInt(numel(Lo));
if reflection
    siglen = 2*datalength;
else
    siglen = datalength;
end
% If the signal length is less than the filter length, need to
% periodize the signal in order to use the DFT algorithm
if siglen < nLo
    ncopies = nLo - siglen;
    Nrep = (1 + ncopies)*siglen;
else
    ncopies = coder.internal.indexInt(0);
    Nrep = siglen;
end
% Allocate storage for a possibly expanded data array.
xx = coder.nullcopy(zeros(Nrep,Nsig,'like',xtmp));

% Copy one siglen worth of elements.
if reflection
    xx(1:datalength,:) = xtmp;
    xx(datalength+1:siglen,:) = xtmp(end:-1:1,:);
else
    xx(1:datalength,:) = xtmp;
end
if ncopies > 0
    % x = [x,repmat(x,1,nLo - siglen)];
    for k = 1:ncopies
        offset = k*siglen;
        for j = 1:siglen
            xx(offset + j,:) = xx(j,:);
        end
    end
end
% Allocate coefficient array
w = zeros(lev+1,siglen,Nsig,'like',x);
% Obtain the DFT of the filters
G = fft(Lo,double(Nrep));
H = fft(Hi,double(Nrep));

%Obtain the DFT of the data
Vhat = fft(xx);
% Main MODWT algorithm
for jj = 1:lev
    [Vhat,What] = modwtdec(Vhat,G,H,jj);
    if isReal
        wtmp = real(ifft(What,size(What,1)));
    else
        wtmp = ifft(What,size(What,1));
    end
    % Truncate data to length of boundary condition
    w(jj,:,:) = wtmp(1:siglen,:,:);
end
if isReal
    wtmp = real(ifft(Vhat,size(Vhat,1)));
else
    wtmp = ifft(Vhat,size(Vhat,1));
end
% Truncate data to length of boundary condition
w(lev+1,:,:) = wtmp(1:siglen,:);

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

function [Vhat,What] = modwtdec(X,G,H,J)
% [Vhat,What] = modwtfft(X,G,H,J)
ONE = int64(1);
N = int64(size(X,1));
upfactor = eml_lshift(ONE,J - 1); % 2^(J - 1);
Vhat = coder.nullcopy(X);
What = coder.nullcopy(X);
for k = 1:N
    idx = 1 + mod(upfactor*(k - 1),N);
    Vhat(k,:) = bsxfun(@times,G(idx),X(k,:));
    What(k,:) = bsxfun(@times,H(idx),X(k,:));
end

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