function [Wout, Vout]=modwtj(Vin, ht, gt, j)
% USAGE
%   [Wt_j, Vout] = modwtj(Vin, ht, gt, j);
%
% INPUTS
%   Vin             - Initial time series, or scaling coefficients for j-1 level.
%   ht              - MODWT avelet filter coefficients.
%   gt              - MODWT Scaling filter coefficients.
%   j               - Level of decomposition.
%
% OUTPUTS
%   Wt_j            - MODWT wavelet coefficients for jth level.
%   Vout            - MODWT scaling coefficients (residuals) for jth level.
%
%   This is matlab code to replace the c routine modwtj.c in wmtsa
%   Joel Hasbrouck (January 8, 2013)
%
L=length(ht);
%fprintf('modwtj j=%d L=%d/%d\n',j,L, length(gt));
Wout=ht(1)*Vin;
Vout = gt(1)*Vin;
N = length(Vin);
kShift = 2^(j-1);
for n=2:L
    shift = (n-1)*kShift;
    k = circshift(1:N,[0, shift]);
    Wout = Wout + ht(n)*Vin(k);
    Vout = Vout + gt(n)*Vin(k);
end

% Replaced c code:
% 
% void modwtj(double *Vin, int N, int j, int L, double *ht, double *gt, 
% 	    double *Wout, double *Vout)
% {
% 
%   int k, n, t;
% 
%   for (t = 0; t < N; t++) {
%     k = t;
%     Wout[t] = ht[0] * Vin[k];
%     Vout[t] = gt[0] * Vin[k];
%     for (n = 1; n < L; n++) {
%       k -= (int) pow(2.0, (double) j - 1.0);
%       if (k < 0) {
% 	k += N;
%       }
%       Wout[t] += ht[n] * Vin[k];
%       Vout[t] += gt[n] * Vin[k];
%     }
%   }
% 
% }

