function [Wdwt,Vdwt] = wavmodwt(X,J,filtertype);

% USAGE: X is the time series to be transformed, J is the level the transform
%        is being computed to, filetype is a .txt or .dat file which contains
%        a vector of the scaling filter coefficients.  The '.txt' is to be 
%        explicitly indluded in the function call.  Examples of scaling filters:
%        
%        haar = [1/sqrt(2) 1/sqrt(2)]';
%        LA8_scaling = [-0.0758,-0.0296,0.4976,0.8037,0.2979,-0.0992,...
%                       -0.0126,0.0322]';
%        C6_scaling = [-0.0157,-0.0727,0.3849,0.8526,0.3379,-0.0727]';
%        
%        Each output matrix, Wdwt, and Vdwt, are JxN matrices of the Jth
%        level wavelet and scaling coefficient vectors, respectively. The
%        'lower' level coefficients are in the top rows of Wdwt and Vdwt.
%        
%        Example: 
%        
%        [Wdwt,Vdwt] = modwt(X,3,'C6_scaling.txt');
%
%        Computes the 3rd level MODWT wavelet and scaling coefficients.
%        
%        See also 
%        MODWTPLOT
%
% Created by Josh Carmichael
% 09.May.2007

%% note: you will see a lot of "+1" because MatLab does not
%% deal with indicies that are less than 1, so adjustments for
%% this must be made.

N = length(X);

% if(floor(log2(N)) ~= log2(N))
%     disp('Make the time series a power of 2 in length');
% end

% the circular shift matrix
T = spdiags(ones(N,1),-1,N,N);
T(1,end) = 1;

% the wavelet scaling filter
g=load(filtertype);
g=g./(sqrt(2));
L=length(g);

% calculate the wavelet filter 
for l=0:L-1
    h(l+1)=(-1)^(l)*g(L-1-l+1);
end

V0=X;

% calculate the MODWT for j=1

j = 1;
for (t=0:N-1);
    k = t;
    W(j,t+1)=h(0+1)*V0(k+1);
    V(j,t+1)=g(0+1)*V0(k+1);
    for n=1:(L-1)
        k=k-2^(j-1);
        if k<0; k=mod(k,N); end
        W(j,t+1)=W(j,t+1)+h(n+1)*V0(k+1);
        V(j,t+1)=V(j,t+1)+g(n+1)*V0(k+1);
    end
end

% calculate MODWT for j=2:J
for j=2:J
    for t=0:(N-1)
        k = t;
        W(j,t+1)=h(0+1)*V(j-1,k+1);
        V(j,t+1)=g(0+1)*V(j-1,k+1);
        for n=1:(L-1)
            k=k-2^(j-1);
            if k<0; k=mod(k,N); end
            W(j,t+1)=W(j,t+1)+h(n+1)*V(j-1,k+1);
            V(j,t+1)=V(j,t+1)+g(n+1)*V(j-1,k+1);
        end
    end
end

Wdwt=W;
Vdwt=V;
return;

% V(1:size(V,1)-1,:)=0;

V(1:J-1,:) = 0;

% calculate inverse MODWT for j=2:J
for j=J:-1:2
    for t=0:(N-1)
        k=t;
	    V(j-1,t+1)=h(0+1)*W(j,k+1)+g(0+1)*V(j,k+1);
        for n=1:(L-1)
            k=k+2^(j-1);
            if k>=N; k=mod(k,N); end
            V(j-1,t+1)=V(j-1,t+1)+h(n+1)*W(j,k+1)+g(n+1)*V(j,k+1);
        end
    end
end


% calculate inverse MODWT for j=1
j=1;
for t=0:(N-1)
    k=t;
    V0(t+1)=h(0+1)*W(j,k+1)+g(0+1)*V(j,k+1);
    for n=1:(L-1)
        k=k+2^(j-1);
        if k>=N; k=mod(k,N); end
        V0(t+1)=V0(t+1)+h(n+1)*W(j,k+1)+g(n+1)*V(j,k+1);
    end
end




Xout=V0;

% plot the input time series, and time series created from inverse MODWT
a=[0 N-1];

L=J;


subplot(L+4,1,L+4)
plot([0:N-1],X,'k')
xlim(a)
%legend('      X_t')
xlabel('time index (t)')

i=1;
    subplot(L+4,1,L+4-i)
    plot([0:N-1],(inv(T)^2)*Wdwt(i,:)','b')
    xlim(a)
    y=max(abs(Wdwt(i,:)))
    ylim([-y y]);
    %legend(['{\itT}^{  -2}W_',num2str(i)])
    
    i=2;
    subplot(L+4,1,L+4-i)
    plot([0:N-1],(inv(T)^5)*Wdwt(i,:)','b')
    xlim(a)
    y=max(abs(Wdwt(i,:)));
    ylim([-y y]);
    %legend(['{\itT}^{  -5}W_',num2str(i)])
    
    i=3;
    subplot(L+4,1,L+4-i)
    plot([0:N-1],(inv(T)^10)*Wdwt(i,:)','b')
    xlim(a)
    y=max(abs(Wdwt(i,:)));
    %ylim([-y y]);
    %legend(['{\itT}^{ -10}W_',num2str(i)])

subplot(L+4,1,1)
plot([0:N-1],(inv(T)^7)*Vdwt(3,:)','r')
ylim([-1 1]);
xlim(a)
%legend(['V_',num2str(3)])

subplot(L+4,1,2)
plot([0:N-1],(inv(T)^3)*Vdwt(2,:)','r')
ylim([-1 1]);
xlim(a)
%legend(['V_',num2str(2)])

subplot(L+4,1,3)
plot([0:N-1],(inv(T))*Vdwt(1,:)','r')
ylim([-1 1]);
xlim(a)
%legend(['V_',num2str(1)])

figure
plot([0:N-1],X)
hold on
plot([0:N-1],Xout,'rx')
legend('Inputted Time Series','Time Series computed from inverse DWT',2)
xlabel('time index (t)')

figure
subplot(3,1,1)
plot([0:N-1],Vdwt(1,:))
hold on
plot([0:N-1],V(1,:),'rx')
xlim(a)
legend('V_1 from modwt','V_1 from invmodwt',-1)
subplot(3,1,2)
plot([0:N-1],Vdwt(2,:))
hold on
plot([0:N-1],V(2,:),'rx')
xlim(a)
legend('V_2 from modwt','V_2 from invmodwt',-1)
subplot(3,1,3)
plot([0:N-1],Vdwt(3,:))
hold on
plot([0:N-1],V(3,:),'rx')
xlim(a)
legend('V_3 from modwt','V_3 from invmodwt',-1)
xlabel('time index (t)')