function [x,Tx]=waverecfix(W,Mr,wvlt)
%
%	A.R. Butz, 8 June 2005
%
%   This program performs the inverse wavelet transform by fixing
%	the W and Mr vectors so that waverec can be used with zero
%   padding.
%
%   x is the signal constituting the inverse transform and Tx < 2
%   indicates the time that x(1) is to be associated with. Since
%   the synthesis filters are anti-causal, Tx will often be negative.
%
%	wvlt is the wavelet being used, e.g.
%	'db3' is the Daubechies wavelet of length 6
%
%	W and Mr give the decomposition into L detail levels,
%   elsewhere designated X(m,l) with l < L+1, and one approximation
%   level, elsewhere designated X(m,L+1). L is not supplied as a
%   parameter; rather, its value is inferred from W and Mr.
%
%	W: detail and approximation coefficients, elsewhere designated
%   X(m,l). A row vector consisting of L+1 segments to represent
%   the L detail segments and the one approximation segment. Its length
%   is sum(Mr(1:L+1)).
%
%	Mr: an integer row vector of length at least L+1 giving the
%	lengths of the L+1 segments of W. W(1:Mr(1)) is the approximation
%	segment, i.e. the values of X(m,L+1), m >= 1. W(Mr(1)+1:Mr(1)+Mr(2))
%   is the values of the lowest detail segment, i.e. the values of X(m,L),
%   m >= 1, etc.. The minimum allowable value of any Mr(k), for k <= L+1,
%   is 1. Values of Mr(k) for k > L+1, if present, are ignored.
%
%   The first values in each segment of W are thought of as occurring
%   at the times 2^L, 2^L, 2^(L-1), 2^(L-2), . . . . , 4, 2 . The time
%   increments associated with each segment are the same, e.g. 2^L for
%   the first two segments. Thus X(m,l) is thought of as occurring at
%   the time m2^l if l <= L, at the time m2^L if l = L+1.
%
%   x(1) returned will usually = 0, but not always. For example, if
%   the wavelet filters have unpadded length 6, and L=1, then
%   x(1) ~= 0. The following session shows this, and also shows that
%   this program cannot be used with symmetric padding:
%
%   W=[0 1];
%   Mr=[1 1];
%   [x,Tx]=waverecfix(W,Mr,'db3')
%   (x(1) ~= 0 and Tx = -3, indicating a time shift of 4)
%   [Wp,Mrp]=wavedec(x,1,'db3')
%   (you get back the starting values in Mr, allowing for the time shift)
%   (now do same with symmetric padding)
%   (comment out the line "dwtmode('zpd') below and then)
%   dwtmode('sym')
%   [x,Tx]=waverecfix(W,Mr,'db3')
%   [Wp,Mrp]=wavedec(x,1,'db3')
%   (you do not get back the starting W)
%
%   That test also fails if waverec is used instead of waverecfix,
%   regardless of whether zero or symmetric padding is used:
%
%   W=[0 1];
%   Mr=[1 1 10];
%   x=waverec(W,Mr,'db3')
%   [Wp,Mrp]=wavedec(x,1,'db3')
%
%   Another example:
%
%  W=sign(randn(1,9))
%  W =
%     1     1    -1     1     1    -1    -1     1     1
%  Mr=[2 3 4];
%  [x,Tx]=waverecfix(W,Mr,'bior2.2)
%
% inverse of a 2 level DWT being computed
%  x =
%         0         0         0    0.0625    0.1250    0.1875    0.2500
%    0.8125    1.3750    0.7607    0.1464    0.8839   -1.2071    1.0607
%    0.5000   -1.1232    0.0821   -1.0714    0.6036    0.3643    0.1250
%    0.0625         0         0
%  Tx =
%   -11
%  [W,Mr]=wavedec(x,2,'bior2.2')
%   W =
%         0         0    0.0000    1.0000    1.0000    0.0000    0.0000
%    0.0000         0         0         0   -0.0000   -1.0000    1.0000
%    1.0000   -0.0000         0         0         0         0         0
%   -0.0000    0.0000   -0.0000   -1.0000   -1.0000    1.0000    1.0000
%         0         0         0         0
%   Mr =
%     9     9    14    24
%
%   Matlab assumes that x(1) in the arguments of wavedec is at time 1,
%   but waverecfix said it should be considered to be at time -11, a
%   difference of 12, resulting in a shift of 6 positions at level 1
%   3 positions at levels 2 and 3.
%
%                     PERIODIC SIGNALS
%
%   This program can be used to find inverse transforms of periodic
%   transforms as follows. First, let X(m,l)=0 for l < L. If that is
%   not true, then the lower values of l can be handled with
%   superposition.
%
%   The four filters of the wavelet set wvlt are assumed to have Matlab
%   length N, which is the lowest even number that is >= all four
%   lengths in the normal sense (e.g. N = 6 for bior2.2, since the filters
%   are of length 3 and 5).
%
%   Let X(m,L) and X(m,L+1) have minimum common period Q = P/2^L,
%   so that {x} has period P. Specify X(m,L) and X(m,L+1) for 1 <=
%   m <= M, where M >= Q + (N-1)*(1-2^(-L)). Specify X(1,l)=0 for
%   l < L. Apply waverecfix in the customary way and then compute
%
%   xperiod=x(2-Tx:P+1-Tx);
%
%   as one period of {x}, with subscript 1 representing time 1.
%
%   If two periods are wanted, take M >= 2*Q + (N-1)*(1-2^(-L)) and say
%
%   xperiod=x(2-Tx:2*P+1-Tx);
%
%   Example. The wavelet set is bior2.2, so N = 6. X(m,l) = 1 for
%   l = 3 and m-2 a multiple of 5, = 0 otherwise. Thus L = 3, P = 40,
%   Q = 5. For two periods of {x} we require M >= 14.375, so M = 15.
%   Then
%
%   W=[0 1 0 0 0];
%   W=[0 W W W 0 0];
%   Mr=[1 15 1 1];
%   [x,Tx]=waverecfix(W,Mr,'bior2.2');
%   x=x(2-Tx:81-Tx);
%
%   yields two periods of {x}.
%
dwtmode('zpd')
%
%   check arguments
%
[rows,lenW]=size(W);
if rows ~= 1
    fprintf('\nW must be a row vector in call to waverecfix\n')
    return
end
[rows,lenMr]=size(Mr);
if rows ~= 1
    fprintf('\nMr must be a row vector in call to waverexfix\n')
    return
end
Lp1=0;
lW=0;
while (lW < lenW) & (Lp1 < lenMr)
    Lp1=Lp1+1;
    if Mr(Lp1)<1
        fprintf('\ndisallowed 0 value in Mr in call to waverecfix\n')
        return
    end
    lW=lW+Mr(Lp1);
end
if (lW ~= lenW) | (Lp1 < 2)
    fprintf('\nW or Mr improperly specified in call to waverecfix\n')
    return
end
L=Lp1-1;
Mr=Mr(1:Lp1);
fprintf('\ninverse of a %i level DWT being computed\n',L)
if W == zeros(size(W))
    x=0;
    Tx=1;
    fprintf('\nthe wavelet transform supplied was zero\n')
    return
end
%
%   Determine the support interval of the inverse transform.
%
nmin=Inf;
nmax=-Inf;
[ga,ha,gs,hs]=wfilters(wvlt);
N=length(ga);
twoL=2^L;
NL=(N-1)*(twoL-1)+1;
m=0;
n=0;
M(1)=m;
while m < Mr(1)
    m=m+1;
    n=n+twoL;
    if W(m) ~= 0
        nmin=min(nmin,n-NL+1);
        nmax=max(nmax,n);
        M(1)=m;
    end
end
for k=2:(L+1)
    M(k)=m;
    n=0;
    while m < sum(Mr(1:k))
        m=m+1;
        n=n+twoL;
        if W(m) ~= 0
            nmin=min(nmin,n-NL+1);
            nmax=max(nmax,n);
            M(k)=m;
        end
    end
    twoL=twoL/2;
    NL=(N-1)*(twoL-1)+1;
end
%
%   determine time associated with x(1) and length of x
%
twoL=2^L;
nstart=1;
nblocks=0;
while nmin < nstart
    nstart=nstart-twoL;
    nblocks=nblocks+1;
end
xlength=nmax-nstart+2;
%
%   structure values of Mr and W appropriately
%
Mrfix(L+2)=xlength;
for k=(L+1):-1:2
    Mrfix(k)=floor((Mrfix(k+1)+N-1)/2);
end
Mrfix(1)=Mrfix(2);
Wfix=zeros(1,nblocks);
if M(1) > 0
    Wfix=[Wfix W(1:M(1))];
end
if length(Wfix) >= Mrfix(1)
    Wfix=Wfix(1:Mrfix(1));
else
    Wfix=[Wfix zeros(1,Mrfix(1)-length(Wfix))];
end
for k=2:(L+1)
    Wtemp=[zeros(1,nblocks) W(sum(Mr(1:(k-1)))+1:M(k))];
    if length(Wtemp) >= Mrfix(k)
        Wtemp=Wtemp(1:Mrfix(k));
    else
        Wtemp=[Wtemp zeros(1,Mrfix(k)-length(Wtemp))];
    end
    nblocks=2*nblocks;
    Wfix=[Wfix Wtemp];
end
%
%   do inverse transform
%
x=waverec(Wfix,Mrfix,wvlt);
if x(xlength) ~= 0
    fprintf('\na check failed in call to waverecfix\n')
    return
end
x=x(1:(xlength-1));
Tx=1-nblocks;


