function [bstats,cstats,sstats]=nt_idxx(fname,iname,bsize,channels_to_keep,nfft,chunksize)
%nt_idxx(fname,iname,bsize,chunksize,nfft) - create an index file to summarize large data file
%
%  fname: name of data file to index
%  iname: name of index file to create [default fname with EXT = .idx]
%  bsize: size of blocks over which to calculate stats [default: 100]
%  channels_to_keep: ignore other channels
%  nfft: fft size for psd [default: 1024]
%  chunksize: size of chunks to read from disk [default: 500000]
%
% If bsize is a struct, the following fields are expected:
%   bsize.basic:  blocksize to calculate basic statistics
%   bsize.cov:  blocksize to calculate covariance [default: none]
%   bsize.psd: blocksize to calculate psd [default: none]
% If bsize is numeric, it refers to basic statistics.
%
% Usage:
%   nt_idx(fname,...): calculate index structs, store in index file
%
%   [bstats,cstats,sstats]=nt_idx(fname,...): return index structs:
%     bstats: basic stats (min, max, etc.)
%     cstats: covariance
%     sstats: psd
%
% NoiseTools
nt_greetings;



assert(nargin>0, '!');
if nargin<2 ; iname=[]; end
if nargin<3 || isempty(bsize)
    bsize.basic=100;
end
if nargin<4; channels_to_keep=[]; end
if nargin<5 || isempty(nfft); nfft=1024; end
if nargin<6 || isempty(chunksize); chunksize=500000; end

if isnumeric(bsize); tmp=bsize; bsize=[]; bsize.basic=tmp; end

% check for FieldTrip
try, ft_version; catch, disp('You must download FieldTrip from http://www.fieldtriptoolbox.org'); return; end

% use separate structs to make it easy to read just one kind of stats from file
bstats=[]; % index structure for basic stats
cstats=[]; % index structure for covariance
sstats=[]; % index structure for spectrogram
bstats.fname=fname; 

% read header 
h=ft_read_header(fname);
bstats.header=h;
bstats.sr=h.Fs;
bstats.nsamples=h.nSamples;
bstats.label=h.label;
bstats.nchans=h.nChans;

if isempty(channels_to_keep); channels_to_keep=1:bstats.nchans; end
if any(channels_to_keep>bstats.nchans); error('!'); end
bstats.channels_to_keep=channels_to_keep;
bstats.nchans=numel(channels_to_keep);

% allocate basic stats arrays:
nbasic=ceil(bstats.nsamples/bsize.basic); % total number of blocs for basic stats
bstats.min=zeros(nbasic,bstats.nchans);
bstats.max=zeros(nbasic,bstats.nchans); 
bstats.mean=zeros(nbasic,bstats.nchans);
bstats.rms=zeros(nbasic,bstats.nchans);
bstats.card=zeros(nbasic,1,'uint32');

chunksize=floor(chunksize/bsize.basic)*bsize.basic; 

% allocate covariance array
if isfield(bsize,'cov')
    tmp=log2(bsize.cov/bsize.basic);
    assert(tmp==round(tmp), ...
        'bsize.cov should be power of 2 times bsize.basic');
    ncov=ceil(bstats.nsamples/bsize.cov);
    cstats.cov=zeros(ncov,bstats.nchans,bstats.nchans);
    cstats.card=zeros(ncov,1,'uint32');
    chunksize=floor(chunksize/bsize.cov)*bsize.cov;
end

% allocate psd array
if isfield(bsize,'psd'); 
    tmp=log2(bsize.psd/bsize.basic);
    assert(tmp==round(tmp), ...
        'bsize.psd should be power of 2 times bsize.basic');
    npsd=ceil(bstats.nsamples/bsize.psd);
    sstats.psd=zeros(npsd,bstats.nchans,nfft/2+1);
    sstats.card=zeros(npsd,1,'uint32');
    sstats.nfft=nfft;
    chunksize=floor(chunksize/bsize.psd)*bsize.psd;
end


foffset=0;
boffset=0;
coffset=0;
soffset=0;

while true
    
    %if file_offset>=i.nsamples; break; end
    
    % read chunk from disk
    begsample=foffset+1;
    endsample=min(foffset+chunksize,bstats.nsamples);
    x=ft_read_data(fname, 'begsample',begsample,'endsample',endsample);
    x=x'; % --> time X channels
    x=x(:,channels_to_keep);
    
    % fold chunk into blocks
    n=floor(size(x,1)/bsize.basic); % number of blocks in this chunk
    xb=x(1:n*bsize.basic,:);
    xb=reshape(xb,[bsize.basic,n,bstats.nchans]);
    bstats.min(boffset+(1:n),:)=min(xb);
    bstats.max(boffset+(1:n),:)=max(xb);
    bstats.mean(boffset+(1:n),:)=mean(xb);
    bstats.rms(boffset+(1:n),:)=sqrt(mean(xb.^2));
    bstats.card(boffset+(1:n),:)=bsize.basic;
    boffset=boffset+n; 

    % extra bit at end of file?
    if size(x,1)>n*bsize.basic
        tmp=x(n*bsize.basic+1,:);
        bstats.min(boffset+1,:)=min(tmp);
        bstats.max(boffset+1,:)=max(tmp);
        bstats.mean(boffset+1,:)=mean(tmp);
        bstats.rms(boffset+1,:)=sqrt(mean(tmp.^2));
        bstats.card(boffset+1,:)=size(tmp,1);
    end
    
    foffset=foffset+n*bsize.basic;

    if isfield(bstats, 'cov')
        n=floor(size(x,1)/bsize.cov); % number of blocks
        xb=x(1:n*bsize.cov,:);        
        xb=reshape(xb,[bsize.cov, n, bstats.nchans]);
        for iBlock=1:n
            tmp=squeeze(xb(:,iBlock,:));
            cstats.cov(coffset+iBlock,:,:) = tmp'*tmp;
            cstats.cardcov(coffset+iBlock,:)=bsize.cov;
        end
        coffset=coffset+size(xb,2);
        if size(x,1)>n*bsize.cov
            tmp=x(n*bsize.cov+1:end,:);
            cstats.cov(coffset+1,:,:)=tmp'*tmp;
            cstats.cardcov(coffset+1,:)=size(tmp,1);
        end              
    end
       
    if isfield(bstats, 'psd')
        n=floor(size(x,1)/bsize.psd); % number of blocks
        xb=x(1:n*bsize.psd,:);        
        xb=reshape(xb,[bsize.psd, n, bstats.nchans]);
        for iBlock=1:n
            tmp=squeeze(xb(:,iBlock,:));
            sstats.psd(soffset+iBlock,:,:) = pwelch(tmp, nfft, 'power')';
            sstats.cardpsd(soffset+iBlock,:,:)=bsize.psd;
        end
        soffset=soffset+size(xb,2);
        if size(x,1)>n*bsize.psd
            tmp=x(n*bsize.psd+1:end,:);
            sstats.psd(soffset+1,:,:) = pwelch(tmp, nfft, 'power')';
            sstats.cardpsd(soffset+1,:)=size(tmp,1);
        end              
    end
    
    nt_whoss
    disp([foffset boffset, coffset]);
    
    if endsample>=bstats.nsamples; break; end;
end
   
if ~nargout
    if isempty(iname)
        [FILEPATH,NAME,EXT]=fileparts(fname);
    else
        [FILEPATH,NAME,EXT]=fileparts(iname);
    end
    iname=[FILEPATH,filesep,NAME,'.idxx'];
    bstats.min=nt_double2int(bstats.min);
    bstats.max=nt_double2int(bstats.max);
    bstats.mean=nt_double2int(bstats.mean);
    bstats.rms=nt_double2int(bstats.rms);
    save(iname, 'bstats', 'cstats', 'sstats');
    clear bstats cstats sstats;
end


    

