function t_hr=fft_make(data_roots,data_suf,fft_dir,fft_root, ...
	T_int,T_step,zpfac,t_lims,f_lims,ph_lims,iihann,ii_app, ...
    spike_parms,TAPE,Jday_start)
%
% Call t_hr=fft_make(data_roots,data_suf,fft_dir,fft_root, ...
% 	T_int,T_step,zpfac,t_lims,f_lims,ph_lims,iihann,ii_app, ...
%   spike_parms,TAPE,Jday_start)
%
% data_roots = Root of data file (cell array of file roots for appending);
% T_int		= FFT Integration time (T_int<0 ==> set T_int such that no zeros inserted: 
%			  zpfac set to zero, T_step set to T_int);
% T_step	= Time between start of FFTs (T_step<0 ==> overlap by T_int/2 (50%), set
%			  iihann=1);
% zpfac		= Zero padding factor: increase FFT by factor of 2^zpfac;
% t_lims	= Time interval to output [t1,t2] in hours of day ([] = all time);
% f_lims	= Frequency band [fmin,fmax] to output ([] = all frequencies, 
%             multiple rows for multiple FFT files with different bands);
% ph_lims	= Phone numbers to output ([] = all phones, [phmin,phmax,-dph], or [ph1;ph2;...phn]);
% iihann	= 1 to shade with a Hann window in time domain;
% ii_app	= 1 to append to output FFT file given by root fft_root;
% spike_parms   = [ii_spike_rem,spike_dB];
% TAPE      = TAPE structure for tapes.  Contains array type, sample rate, calibration factor, ...
% Jday_start= Day of year to start
% Outputs:
% t_hr		= Hour of day last processed;
%
% COPYRIGHT (c) 2008
% 
% Applied Research Laboratories, The University of Texas at Austin
% (ARL:UT)
% All rights reserved.
% 
% This material may be reproduced by or for the U.S. Government
% only. This notice must appear in all copies of this file and its
% derivatives. Further distribution outside the U.S. Government is
% prohibited without the explicit written authorization from the
% Director of the Environmental Sciences Laboratory, ARL:UT.
% 
% THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND,
% EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO THE
% IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
% PURPOSE, OR NONINFRINGEMENT.
% 
% IN NO EVENT WILL ARL:UT BE LIABLE FOR ANY DAMAGES WHATSOEVER ARISING
% OUT OF THE USE OR INABILITY TO USE THIS SOFTWARE, EVEN IF ARL:UT HAS
% BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.

if zpfac<0 | round(zpfac)~=zpfac | zpfac>4,
	zpfac,
	error('USER ERROR: Bad zpfac: FFT size to be increased by 2^zpfac (0<=zpfac<=4)'); 
end
%
if iscell(data_roots),
    nfile_app=length(data_roots); jfile_app=1;
    data_root=data_roots{jfile_app};
else
    nfile_app=0; jfile_app=1;
    data_root=data_roots;
end
data_file=[data_root,data_suf];
data_file; nfile_app;

TAPE.data_file=data_file;
[TAPE] = ts_open(TAPE);
if TAPE.status <= 0
    return;
end

t_hr=0;
nchan=TAPE.nchan;
fs=TAPE.fs;
calfac=TAPE.calfac;
f0_cal=100; dB_octave=0; ii_dB_octave=0;    % For f-dependent calibrations;

iiph=adi_lim_eval(ph_lims);
if isempty(iiph), iiph=[1:nchan]'; end
if isfield(TAPE,'phone2chan'),
    if max(iiph)>length(TAPE.phone2chan), 
        disp('Error: Phone #s asked for exceed MCS.TAPE.phone2chan!!'); 
    end
    iichan=TAPE.phone2chan(iiph);
else
    disp('NOTE: No channel to phone mapping available for this data type ...');
    iichan=iiph;
end
nrec=length(iichan);

if any(calfac<=0), error('USER ERROR: calibration factor must be > 0'); end


% For T_int<0, set T_int such that no zeros used to pad:
dt=1/fs;
if T_int<0,
   T_int=dt * 2^ceil(log2(-T_int/dt));
end
% For abs(T_step)=999, set T_step=T_int/2, iihann=1; iizp=0;
if abs(T_step)==999,
	T_step=T_int/2; iihann=1; zpfac=0;	
elseif T_step<=0 & T_step>=-200,
	T_step=T_int*(abs(T_step)/100);
elseif T_step<0,
	error(['Illegal T_step: ',num2str(T_step)]);
end

nT_int=round(T_int/dt); T_int=nT_int*dt;
% Use floor instead of round so that more data overlapped:
nT_step=floor(T_step/dt); T_step=nT_step*dt;
% Shading in time domain:
%%iit=[1:nT_int]'; t_out=(iit-1)*dt;
iit=[1:nT_int]'; t_out=(iit-.5)*dt;
if iihann,
	hann=.5*(1 - cos((2*pi/nT_int)*[.5:1:nT_int]'));
% This normalization preserves peak power in CW
%%%hann_norm=sum(hann)/length(hann);
% This is correct normalization that preserves total energy in the signal
% (difference of 20 * log10 ( sqrt(3/8) / (1/2) ) = 1.76 dB
% Hanning windows widen the mainlobe, reduce it by 1.76 dB, and lower the sidelobes.
	hann_norm=sqrt(sum(hann.^2)/length(hann));
	hann_vec=(hann/hann_norm);
end
%
nfft=2^(zpfac + ceil(log2(T_int/dt)) );
nT_zero=nfft - nT_int;
T_fft=nfft/fs; df=fs/nfft;
% disp(['T_int = ',num2str(T_int),'; T_step = ',num2str(T_step),'; T_zero = ',num2str(nT_zero*dt),';']);

ii_copy=nT_step<nT_int; nT_skip=max(0,nT_step-nT_int);
nT_read0=min(nT_step,nT_int);       % Number of samples to read in midst of sequences;
nT_read0_nchan = nT_read0*nchan;

if ii_copy,
    ii_to=[1:nT_int-nT_step];
    ii_from=[nT_step+1:nT_int];
end
ii_end=[max(nT_int-nT_step+1,1):nT_int];

SPEC.nfft=nfft; SPEC.fs=fs;
SPEC.Tw=T_fft; SPEC.nT_int=nT_int;
SPEC.T_step=T_step; SPEC.iihann=iihann;
if isempty(fft_root),
	fft_root=[data_root,'_',int2str(T_int),'_',int2str(T_step)]; 
end
if ischar(fft_root), fft_root={fft_root}; end
%
fft_norm=sqrt(2/(nT_int*nfft));

if isempty(f_lims), f_lims=[0,fs/2]; end
[nfft_file,ncol]=size(f_lims);
if ncol~=2, error('USER ERROR in fft_make: f_lims must be nfft_file by 2 '); end

if isfield(TAPE,'cal_table')
    cal_table=TAPE.cal_table;
    phase_table=TAPE.phase_table;
else
    cal_table=[];
end
ii_cal_table=~isempty(cal_table);
if length(dB_octave)==1 & nrec>1, dB_octave=dB_octave*ones(1,nrec); end
if any(dB_octave~=0), ii_dB_octave=1; end
%
ii_cal_file=TAPE.ii_cal & ~isempty(TAPE.cal_root);
if ii_cal_file, 
%     eval(['load ',TAPE.cal_root,'.mat']);
    load([TAPE.cal_root,'.mat']);
%     if any(CAL.dims([1,2:4])) > 1, 
%     if any(CAL.dims([3:4]) > 1), 
%         error('CAL File data not properly dimensioned!!'); 
%     end
    nph_data=length(TAPE.phone2chan);
    if CAL.dims(2)~=nph_data,
        error(['CAL File data should have nrec calibration numbers: ',int2str([CAL.dims(2),nph_data])]);
    end
end
%
nband=size(f_lims,1);
if length(fft_root)~=nband,
    disp(f_lims), disp(fft_root),
    error('fft_root must be cell array having nband entries, where nband=# rows of f_lims');
end
freq_cal=[]; ii_cal=ii_cal_table | ii_dB_octave | ii_cal_file;
for jfft=1:nfft_file,
	nf1=round(f_lims(jfft,1)/df) + 1; nf2=round(f_lims(jfft,2)/df) + 1;
	f_ax=[nf1-1:nf2-1]'*df; nfreq=nf2-nf1+1;
    log_fax=log2(f_ax);
    freq_cal{jfft}=ones(nfreq,nrec);
    if ii_cal_table,
        for jrec=1:nrec,
            dB_cal=interp1(log2(cal_table{jrec}(:,1)),cal_table{jrec}(:,2),log_fax,'linear','extrap');
            phase_cal=interp1(phase_table{jrec}(:,1),phase_table{jrec}(:,2),f_ax,'linear','extrap');
            freq_cal{jfft}(:,jrec)=10.^(dB_cal/20) .* exp(i*phase_cal*pi/180);       %figure; plot(f_ax,freq_cal);
        end
    elseif ii_dB_octave,
        for jrec=1:nrec,
            dB_cal=log2(max(.1,f_ax/f0_cal)) * dB_octave(jrec);
            freq_cal{jfft}(:,jrec)=10.^(dB_cal/20);       %figure; plot(f_ax,freq_cal);
        end
    end
    if ii_cal_file,
        freq_cal_file{jfft}=ones(nfreq,nrec);
        cal_file_fac=zeros(CAL.dims(1),CAL.dims(2));   % Fill cal_file_fac with zeros the size of all phones in array;
        cal_file_fac(:,CAL.iiph)=squeeze(CAL.data(:,:,1,1));  % Set cals for phones in CAL.data;
        if ~isempty(CAL.phase_cal)
            cal_file_phase(:,CAL.iiph)=squeeze(CAL.phase_cal(:,:,1,1));  % Set cals for phones in CAL.data;
        end
        for jrec=1:nrec,
            dB_cal=interp1(log2(CAL.scales{1}),cal_file_fac(:,iichan(jrec)),log_fax,'linear','extrap');
            freq_cal_file{jfft}(:,jrec)=10.^(dB_cal/20);
            if ~isempty(CAL.phase_cal)
                phase_cal=interp1(log2(CAL.scales{1}),cal_file_phase(:,jrec),log_fax,'linear','extrap');
                freq_cal_file{jfft}(:,jrec)=freq_cal_file{jfft}(:,jrec) .* exp(i*phase_cal*pi/180);
            end
        end
        freq_cal{jfft} = freq_cal{jfft} .* freq_cal_file{jfft};
    end
%     if ii_cal_file,
%         cal_file_fac=zeros(1,nph_data);             % Fill cal_file_fac with zeros the size of all phones in array;
%         cal_file_fac(CAL.iiph)=squeeze(CAL.data(1,:,1,1));   % Set cals for phones in CAL.data;
%         for jrec=1:nrec,
% % Fill freq_call with cals for requested phones:
%             freq_cal{jfft}(:,jrec)=freq_cal{jfft}(:,jrec) .* 10.^(cal_file_fac(iiph(jrec))/20);
%         end
% 	end

    % Compute imap for inverse FFT calculations:
	SPEC.nfreq=nfreq;
	SPEC.imap=zeros(SPEC.nfreq,3);
	SPEC.imap(:,1)=[nf1:nf2]';
	SPEC.imap(:,2)=20 + [1:2:2*SPEC.nfreq]';
	SPEC.imap(:,3)=SPEC.imap(:,2) + 1;
	REC.nrec=nrec; REC.zr=zeros(1,nrec);
	iif_out{jfft}=[nf1:nf2];
%     if nfft_file>1,
% % tms:  simple frequency labeling for reproduction in other programs
%         flo=round(f_lims(jfft,1)); fhi=round(f_lims(jfft,2));
% %         fft_rootx=[fft_root,'_',int2str(flo),'_',int2str(fhi)];
% % EKW: Even simpler for shorter names:
%         fft_rootx=[fft_dir,fft_root,'_',int2str(.5*(flo+fhi)),'Hz'];
%         disp(['Output FFT file root = ',fft_rootx]);
%     else
    fft_rootx=[fft_dir,fft_root{jfft}];
%     end
	% Open output FFT file:
	FFT{jfft}=fft_write1_open(fft_rootx,SPEC,REC,ii_app);
	% Write to output file:
	out_file=[fft_rootx,'.txt'];
	if ii_app==0,
		fid_out=fopen(out_file,'wt');
	else
		fid_out=fopen(out_file,'at');		% Open in append mode;
	end

	if fid_out < 0,
		disp(['ERROR: OUTPUT FILE ',out_file,' COULD NOT BE OPENED FOR WRITING.']);
		disp(['       CHECK THAT IT IS NOT BEING EDITED']);
		return
	end
    cl=clock;
	fprintf(fid_out,['FFT_MAKE Run Date = ',date,'; Time = %i:%i:%2.2f','\n'],cl(4:6));
    if isfield(TAPE,'rec_name'), rec_name=TAPE.rec_name{1}; else, rec_name='N/A'; end
	fprintf(fid_out,'REC = %s;\n',rec_name);
    if isfield(TAPE,'src_name') & length(TAPE.src_name)>=jfft, src_name=TAPE.src_name{jfft}; else, src_name='N/A'; end
    fprintf(fid_out,'SRC = %s;\n',src_name);
	fmtx='Array type = %s;\nFile name = %s;\nOutput FFT file = %s\n'; 
	fprintf(fid_out,fmtx,TAPE.array_type,data_file,fft_rootx);
	fmtx='fs = %1.0f; calfac = %5.1f; T_int = %5.1f; T_step = %5.1f; zpfac = %2i; T_lims = %5.2f,%5.2f\n';
	fprintf(fid_out,fmtx,fs,calfac(1),T_int,T_step,zpfac,t_lims);
	fmtx='f_lims = %2.1f,%2.1f\n'; 
	fprintf(fid_out,fmtx,f_lims');
	fmtx=['iihann = %1i; Jday_start = %3.3i; Jhour_start = %6.3f\n'];
	fprintf(fid_out,fmtx,[iihann,Jday_start,t_lims(1)]);
	fprintf(fid_out,'Phone Numbers   = ');
	fprintf(fid_out,'%3i ',iiph); fprintf(fid_out,'\n');
	fprintf(fid_out,'Channel Numbers = ');
	fprintf(fid_out,'%3i ',iichan); fprintf(fid_out,'\n');
	if ~isfield(TAPE,'blocksize'), TAPE.blocksize=0; end
	if TAPE.isatape,
        fprintf(fid_out,'blocksize = %5i; nsamp = %5i',[TAPE.blocksize,TAPE.nsamp]);
    end
	fclose(fid_out);
end

if TAPE.isatape == 0
	ii_done=0; 
	while ~ii_done
		if t_lims(1) > TAPE.Dhour_end
            jfile_app=jfile_app + 1;
            if jfile_app<=nfile_app,
                data_root=data_roots{jfile_app}; data_file=[data_root,data_suf];                
                [TAPE]=data_file_switch(TAPE,data_file);
                if TAPE.status <= 0
                    disp(['fft_make:  data_file_switch returned <= 0 for ' TAPE.datafile]);
                    keyboard
                    return
                end
            else
                TAPE=ts_close(TAPE);
                disp([t_lims(1),t_lims]);
                error(['Start time not found in specified file(s)']);
            end
		else
            ii_done=1;
		end
	end
end

% Position to start time
[TAPE]=ts_tpos(TAPE,Jday_start,t_lims(1));
if TAPE.status <= 0
    disp(['fft_make:  ts_tpos returned <= 0 for ' TAPE.data_file]);
    return
end

Jday=TAPE.day;
Jhour=TAPE.Dhour;
t_hr=Jhour;

% Read in first FFT's worth:
[pser,TAPE,count] = ts_read(TAPE,nT_int);
if count < nT_int*nchan,
    jfile_app=jfile_app + 1;
    if jfile_app<=nfile_app,  % will this also catch ii_multi and isatape?
        [pser,TAPE] = inc_read_switch_files(TAPE,pser,jfile_app,data_roots,data_suf,count, ...
            nT_int,Jday,t_hr,T_int,T_step,t_lims);
        if TAPE.status <= 0, error('Error calling inc_read_switch_files'); end
    else
		error(['Stopping because end of file reached on 1st FFT: times = ', ...
                sprintf('%3.3f,%3.3f',[t_hr,t_lims(2)])]);
    end
end
%
pser=pser(iichan,:).';
% if TAPE.status <= 0
%     disp(['fft_make:  initial read of nT_int samples failed for ' TAPE.data_file])
%     TAPE=ts_close(TAPE);
%     return;
% end
% nT_read=count/nchan;
% if nT_read<nT_int,
% 	disp('Data file does not contain enough data for one FFT: ')
% 	nT_int,nT_read,
%     keyboard
% 	return
% end

ptot=zeros(nfft,nrec); 
if iihann, pshade=zeros(nfft,1); end

if length(calfac)>1,
	for jrec=1:nrec,
		ptot(iit,jrec)=calfac(jrec)*pser(:,jrec);
	end
else
	ptot(iit,:)=calfac*pser;
end

ii_spike_rem=spike_parms(1); spike_dB=spike_parms(2);
if ii_spike_rem, 
    spike_lin=10^(spike_dB/20), 
end

tic

t_cpu0=clock; nseq=FFT{1}.nseq;
while 1==1,
	nseq=nseq + 1;
	t0=t_hr*3600;		% t0 in seconds of day;
    if ii_spike_rem,
        ptot=ptot - mean(ptot(:));
        p_samples=ptot(1:10:prod(size(ptot)));
        spike_lin=15*sqrt(median(p_samples.^2));
        iisp=find(abs(ptot)>spike_lin);
        if ~isempty(iisp),
            disp(['Spikes found ',num2str(100*length(iisp)/prod(size(ptot))),' percent of the time.']);
            ptot(iisp)=0;
        end
    end
	if iihann,
		for jrec=1:nrec,
			pshade(iit)=ptot(iit,jrec).*hann_vec;
			spec=fft(pshade,nfft);
            for jfft=1:nfft_file,
    			cspec=fft_norm*spec(iif_out{jfft});
                if ii_cal, cspec=cspec.*freq_cal{jfft}(:,jrec); end
			    stat=fft_write1(FFT{jfft},cspec,nseq,jrec,Jday,t0,nseq,jrec);
            end
		end
	else
		for jrec=1:nrec,
			spec=fft(ptot(:,jrec));
            for jfft=1:nfft_file,
    			cspec=fft_norm*spec(iif_out{jfft});
                if ii_cal, cspec=cspec.*freq_cal{jfft}(:,jrec); end
			    stat=fft_write1(FFT{jfft},cspec,nseq,jrec,Jday,t0,nseq,jrec);
            end
		end
	end

	if t_hr>t_lims(2), 
		disp(['Stopping because end time reached: ', ...
			sprintf('%3.3f,%3.3f',[t_hr,t_lims(2)])]); 
		break; 
	end
    
	t_hr=t_hr + T_step/3600;
    
	% Skip data if T_step>T_int:
    if nT_skip>0,
        [pser,TAPE,count] = ts_read(TAPE,nT_skip);
        % incomplete read
        if count < nT_skip*nchan
            jfile_app=jfile_app + 1;
            if jfile_app<=nfile_app
                [pser,TAPE] = inc_read_switch_files(TAPE,pser,jfile_app,data_roots,data_suf,count, ...
                    nT_skip,Jday,t_hr,T_int,T_step,t_lims);
                if TAPE.status <= 0
                    break;
                end
            else
        		disp(['Stopping because end of file reached: times = ', ...
                        sprintf('%3.3f,%3.3f',[t_hr,t_lims(2)])]);
                break;
            end
        end
    end
    
	% Read in next chunk of data:
    [pser,TAPE,count] = ts_read(TAPE,nT_read0);
    % incomplete read
    if count < nT_read0_nchan,
        jfile_app=jfile_app + 1;
        if jfile_app<=nfile_app,  % will this also catch ii_multi and isatape?
            [pser,TAPE] = inc_read_switch_files(TAPE,pser,jfile_app,data_roots,data_suf,count, ...
                nT_read0,Jday,t_hr,T_int,T_step,t_lims);
            if TAPE.status <= 0
                break;
            end
        else
    		disp(['Stopping because end of file reached: times = ', ...
                    sprintf('%3.3f,%3.3f',[t_hr,t_lims(2)])]);
            break;
        end
    end

    pser=reshape(pser,[nchan,nT_read0]);
    pser=pser(iichan,:).';

    if ii_copy, ptot(ii_to,:)=ptot(ii_from,:); end

	if length(calfac)>1,
		for jrec=1:nrec,
			ptot(ii_end,jrec)=calfac(jrec)*pser(:,jrec);
		end
	else
		ptot(ii_end,:)=calfac*pser;
	end
	if rem(nseq,5)==1,
		ctime=etime(clock,t_cpu0)/60;
		ptime=ctime*(t_lims(2)-t_lims(1))/(t_hr-t_lims(1));
		disp(['Finished seq ',sprintf('%4i',nseq),' at hour ', ...
		sprintf('%3.3f',t_hr),'; Cur/Proj CPU times = ', ...
        sprintf('%2.2f/%2.2f',[ctime,ptime]),' min.']);
	end
    pause(.001)
end

toc


TAPE = ts_close(TAPE);

for jfft=1:nfft_file, 
    fclose(FFT{jfft}.fft_fid); 
end
disp(['# Sequences output: ',int2str(nseq)]);
ctime=etime(clock,t_cpu0)/60;
disp(['CPU Time = ',sprintf('%2.2f',ctime),' minutes']);



function [TAPE]=data_file_switch(TAPE,data_file);
%
% Closes current input file and opens next one;
%

TAPE=ts_close(TAPE);
if TAPE.status <= 0
    disp(['data_file_switch:  ts_close returned <= 0 for ' TAPE.data_file]);
    disp('continuing');
end
TAPE.data_file=data_file;
TAPE=ts_open(TAPE);


function [pser,TAPE] = inc_read_switch_files(TAPE,pser,jfile_app,data_roots,data_suf,count, ...
    nT_read0,Jday,t_hr,T_int,T_step,t_lims);
% Call:  [pser,TAPE] = inc_read_switch_files(TAPE,pser,jfile_app,nfile_app,data_roots,data_suf,count, ...
%    nT_read0,Jday,t_hr,T_int,T_step,t_lims);
% Handle partial read and switching files for fft_make

dt=1/TAPE.fs;
disp(['Opening next data file to process continuous data: T1,T2 = ', ...
        sprintf('%3.3f,%3.3f',[t_hr,t_lims(2)])]);
data_root=data_roots{jfile_app}; data_file=[data_root,data_suf];     
[TAPE]=data_file_switch(TAPE,data_file);

count=max(count,0);
nt_ok=count/TAPE.nchan;
nt_left=nT_read0 - nt_ok;
Jhour=TAPE.Dhour; Jhour_end=TAPE.Dhour_end;
disp(['Start Julian day = ',int2str(Jday),'; T1,T2 = ',sprintf('%3.3f,%3.3f',[Jhour,Jhour_end])]);
t_cont=t_hr + (T_int-T_step)/3600 + nt_ok*dt/3600;
nT_t1=round(3600*(t_cont-Jhour)/dt);	% Number of samples to go into data file;
if nT_t1>0,
    %tms:  this skipping has caused problems; just warn
    disp(['Warning:  nT_tl >0, indicating file overlap, but no action taken'])
elseif nT_t1<0,
    disp(['Found no overlap between files.  Data may be discontinuous. Gap = ', ...
        num2str(nT_t1*dt),' s.']);
else
    disp(['Files appear to be continuous'])
end

if nt_left > 0
    [pser_left,TAPE,count] = ts_read(TAPE,nt_left);
	if count < nt_left*TAPE.nchan
        disp(['fft_make: trying to read nt_left samples failed for ' TAPE.data_file])
        TAPE=ts_close(TAPE);
        TAPE.status=-1;
        return;
	end
    if ~isempty(pser)
        pser=[pser, pser_left];
    else
        pser=pser_left;
    end
end

TAPE.status=1;
