0001 function [bstats,wstats,cstats,sstats]=nt_idxx(fname,iname,blksize,channels_to_keep,nfft,chunksize)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027 nt_greetings;
0028
0029
0030
0031 assert(nargin>0, '!');
0032 if nargin<2 ; iname=[]; end
0033 if nargin<3 || isempty(blksize)
0034 blksize.wav=100;
0035 end
0036 if nargin<4; channels_to_keep=[]; end
0037 if nargin<5 || isempty(nfft); nfft=1024; end
0038 if nargin<6 || isempty(chunksize); chunksize=500000; end
0039
0040 if isnumeric(blksize); tmp=blksize; blksize=[]; blksize.wav=tmp; end
0041 if ~isempty(iname) && ~ischar(iname); error('!'); end
0042
0043
0044 try, ft_version; catch, disp('You must download FieldTrip from http://www.fieldtriptoolbox.org'); return; end
0045
0046
0047 bstats=[];
0048 wstats=[];
0049 cstats=[];
0050 sstats=[];
0051 bstats.fname=fname;
0052
0053
0054 h=ft_read_header(fname);
0055 bstats.header=h;
0056 bstats.sr=h.Fs;
0057 bstats.nsamples=h.nSamples;
0058 bstats.label=h.label;
0059 bstats.nchans=h.nChans;
0060
0061 if isempty(channels_to_keep); channels_to_keep=1:bstats.nchans; end
0062 if any(channels_to_keep>bstats.nchans); error('!'); end
0063 bstats.channels_to_keep=channels_to_keep;
0064 bstats.nchans=numel(channels_to_keep);
0065
0066
0067 nbasic=ceil(bstats.nsamples/blksize.wav);
0068 wstats.min=zeros(nbasic,bstats.nchans);
0069 wstats.max=zeros(nbasic,bstats.nchans);
0070 wstats.mean=zeros(nbasic,bstats.nchans);
0071 wstats.rms=zeros(nbasic,bstats.nchans);
0072 wstats.card=zeros(nbasic,1,'uint32');
0073
0074 chunksize=floor(chunksize/blksize.wav)*blksize.wav;
0075
0076
0077 if isfield(blksize,'cov')
0078 tmp=log2(blksize.cov/blksize.wav);
0079 assert(tmp==round(tmp), ...
0080 'blksize.cov should be power of 2 times blksize.wav');
0081 ncov=ceil(bstats.nsamples/blksize.cov);
0082 cstats.cov=zeros(ncov,bstats.nchans,bstats.nchans);
0083 cstats.card=zeros(ncov,1,'uint32');
0084 chunksize=floor(chunksize/blksize.cov)*blksize.cov;
0085 end
0086
0087
0088 if isfield(blksize,'psd')
0089 if blksize.psd < nfft; error('!'); end;
0090 tmp=log2(blksize.psd/blksize.wav);
0091 assert(tmp==round(tmp), ...
0092 'blksize.psd should be power of 2 times blksize.wav');
0093 npsd=ceil(bstats.nsamples/blksize.psd);
0094 sstats.psd=zeros(npsd,bstats.nchans,nfft/2+1);
0095 sstats.card=zeros(npsd,1,'uint32');
0096 sstats.nfft=nfft;
0097 chunksize=floor(chunksize/blksize.psd)*blksize.psd;
0098 end
0099
0100
0101 foffset=0;
0102 boffset=0;
0103 coffset=0;
0104 soffset=0;
0105
0106 while true
0107
0108
0109
0110
0111 begsample=foffset+1;
0112 endsample=min(foffset+chunksize,bstats.nsamples);
0113 x=ft_read_data(fname, 'begsample',begsample,'endsample',endsample);
0114 x=x';
0115 x=x(:,channels_to_keep);
0116
0117
0118 n=floor(size(x,1)/blksize.wav);
0119 xb=x(1:n*blksize.wav,:);
0120 xb=reshape(xb,[blksize.wav,n,bstats.nchans]);
0121 wstats.min(boffset+(1:n),:)=min(xb);
0122 wstats.max(boffset+(1:n),:)=max(xb);
0123 wstats.mean(boffset+(1:n),:)=mean(xb);
0124 wstats.rms(boffset+(1:n),:)=sqrt(mean(xb.^2));
0125 wstats.card(boffset+(1:n),:)=blksize.wav;
0126 boffset=boffset+n;
0127
0128
0129 if size(x,1)>n*blksize.wav
0130 tmp=x(n*blksize.wav+1:end,:);
0131 wstats.min(boffset+1,:)=min(tmp);
0132 wstats.max(boffset+1,:)=max(tmp);
0133 wstats.mean(boffset+1,:)=mean(tmp);
0134 wstats.rms(boffset+1,:)=sqrt(mean(tmp.^2));
0135 wstats.card(boffset+1,:)=size(tmp,1);
0136 end
0137
0138
0139 foffset=foffset+n*blksize.wav;
0140
0141 if ~isempty(cstats) && isfield(cstats, 'cov')
0142 n=floor(size(x,1)/blksize.cov);
0143 xb=x(1:n*blksize.cov,:);
0144 xb=reshape(xb,[blksize.cov, n, bstats.nchans]);
0145 for iBlock=1:n
0146 tmp=squeeze(xb(:,iBlock,:));
0147 tmp=nt_demean(tmp);
0148 cstats.cov(coffset+iBlock,:,:) = tmp'*tmp;
0149 cstats.cardcov(coffset+iBlock,:)=blksize.cov;
0150 end
0151 coffset=coffset+size(xb,2);
0152 if size(x,1)>n*blksize.cov
0153 tmp=x(n*blksize.cov+1:end,:);
0154 tmp=nt_demean(tmp);
0155 cstats.cov(coffset+1,:,:)=tmp'*tmp;
0156 cstats.cardcov(coffset+1,:)=size(tmp,1);
0157 end
0158 end
0159
0160 if ~isempty(sstats) && isfield(sstats, 'psd')
0161 n=floor(size(x,1)/blksize.psd);
0162 xb=x(1:n*blksize.psd,:);
0163 xb=reshape(xb,[blksize.psd, n, bstats.nchans]);
0164 for iBlock=1:n
0165 tmp=squeeze(xb(:,iBlock,:));
0166 tmp=nt_demean(tmp);
0167 sstats.psd(soffset+iBlock,:,:) = pwelch(tmp, nfft, 'power')';
0168 sstats.cardpsd(soffset+iBlock,:,:)=blksize.psd;
0169 end
0170 soffset=soffset+size(xb,2);
0171 if size(x,1)>n*blksize.psd
0172 tmp=x(n*blksize.psd+1:end,:);
0173 if size(tmp,1)<nfft; break; end
0174 tmp=nt_demean(tmp);
0175 sstats.psd(soffset+1,:,:) = pwelch(tmp, nfft, 'power')';
0176 sstats.cardpsd(soffset+1,:)=size(tmp,1);
0177 end
0178 end
0179
0180 nt_whoss
0181 disp([num2str(foffset), '/', num2str(h.nSamples), ' (', num2str(foffset/h.nSamples*100), '%)']);
0182 disp([boffset, coffset, soffset]);
0183
0184 if endsample>=bstats.nsamples; break; end;
0185 end
0186
0187 if ~nargout
0188 if isempty(iname)
0189 [FILEPATH,NAME,EXT]=fileparts(fname);
0190 if isempty(FILEPATH); FILEPATH=pwd; end
0191 if ~exist([FILEPATH,filesep,'idxx'], 'dir')
0192 mkdir([FILEPATH,filesep,'idxx']);
0193 end
0194 iname=[FILEPATH,filesep,'idxx',filesep,NAME,EXT,'.idxx'];
0195 end
0196 wstats.min=nt_double2int(wstats.min);
0197 wstats.max=nt_double2int(wstats.max);
0198 wstats.mean=nt_double2int(wstats.mean);
0199 wstats.rms=nt_double2int(wstats.rms);
0200 save(iname, 'bstats', 'wstats','cstats', 'sstats','-v7.3');
0201 clear bstats wstats cstats sstats;
0202 end
0203
0204
0205
0206