function D=coralMRA(S,varargin)
% USAGE
% D=coralMRA(S,varargin)
% D=coralMRA(S,J,filtertype);
%
% INPUT:  
% S: A coral structure S of seismograms.  Works best when the number of
%    data samples per seismogram is less than 10,000.
% J: The 'level' of the transform.  A level 4 transform returns 4 detail
%    coefficients, 1 smooth.  The 1st coefficient contains energy in the
%    frequency band [100,50]Hz.  The 2nd coefficient contains data in the
%    [50,25]Hz band, the 3rd coefficient contains data in the [25,12.5]Hz
%    band, and so on.  The smooth gives the base-line drift.  RECOMMENDED:
%    Use a value between 4 and 6. DEFAULT: 5.
% filtertype: Which scaling filter coefficient you use.  DEFAULT: C6.
%             Other choices are available, but not yet coded in.  Just use
%             C6!
% OUTPUT:
% D:  A structure with fields: details, smooth, periodogram.  The
%     periodogram is a length(J) vector giving the square of the 2-norm of
%     each wavelet coefficient (not detail coefficient).  The sum of the
%     periodogram coefficients is the sample variance.
%% Description
% Wavelet detail coefficients give a multi-resolution analysis of a time
% series by decomposing the original time series into a dyadic sum of
% band-pass filtered time series.  Each sequence gives the contribution of
% the signal from the respective band.  For example, a level 4 MRA produces
% 5 time series. 
%%
[mx,nx]=size(S);


if(mx<nx);
    disp('A ROW vector of structures?  Are you kidding?');
    disp('Well, you are going to GET a COLUMN vector of structures back');
    S=S';
    [mx,nx]=size(S);
end;


for(i=1:mx);
    
    X=S(i).data;
    
    % I HATE ROW VECTORS
    [xrow,xcol]=size(X);
    if(xcol>xrow); X=X'; end;
    
    sint=S(i).recSampInt;
    
    if(nargin==1);
        filtertype='C6_scaling.txt';
        J=5;
    end;
    
    if(nargin==2);
        J=varargin{1};
        filtertype='C6_scaling.txt';
    end;
    
    if(nargin==3);
        J=varargin{1};
        filtertype=varargin{2};
    end;
    
    % GET WAVELET COEFFICIENTS
    [W,V] = wavmodwt(X,J,filtertype);
    N=length(W);
    
    % TIME SCALE INDICES
    j=1:J+1;
    tau=2.^(j)*sint;
    scale=log2(tau);
    
    % ENERGY PER UNIT SCALE
    P=[W;V(end,:)]*[W;V(end,:)]'; P=(1/N)*diag(P); svar=sum(P);
    D(i).periodogram=P;

    % COMPUTES EACH DETAIL COEFFICIENT PER SEISMOGRAM
    for(k=J:-1:1)
        Wv=W(k,:);
        if(k==J);
            Dj=wavmodwtdet(Wv,k,filtertype);
        else
        Djtemp=wavmodwtdet(Wv,k,filtertype);
        Dj=cat(1,Dj,Djtemp);
        end;
    end;
  
    detailSum=sum(Dj);
    
    % THE SMOOTH IS THE DATA LESS THE DETAILS
    
    SJ=X'-detailSum;

    D(i).details=flipud(Dj);
    D(i).smooth=SJ;
end;