function [LGRAM,FFT]=cts2fft(cts_root,fft_root,T_int,T_step,zpfac,t_lims,f_lims, ...
	ph_lims,iihann,iilg,cal_root);
%
% Call [LGRAM,FFT]=cts2fft(cts_root,fft_root,T_int,T_step,zpfac,t_lims,f_lims, ...
%	ph_lims,iihann,iilg);
% T_int		= FFT Integration time (T_int<0 ==> set T_int>=|T_int| such power of 2 NFFT);
% T_step	= Time between starts of FFTs (T_step=999 ==> 50% ovl, Hann=1, zpfac=0,
%			  -100<T_step<=0 ==> Use |T_step| as percent overlap)
%			  <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 to output ([] = all frequencies);
% ph_lims	= Phone numbers to output ([] = all phones);
% iihann	= 1 to shade with a Hann window in time domain;
% iilg		= 1 to call lgram_3d, 2 to view LG also;
% cal_root  = [cal_root '_CAL.mat'] = calibration file for VLAD
%
% See:  cts_details, cts_open, cts_read, cts_write, cts_write_open, cts_combine, cts2fft, cts2spec, 
%       cts_from_real, cts_shot, cts_wf_make, fft2cts, mcs_gui.
%
if nargin<2, fft_root=[]; end
if nargin<3, T_int=-50; end
if nargin<4, T_step=999; end
if nargin<5, zpfac=1; end
if nargin<6 | isempty(t_lims), t_lims=[-24,48]; end
if nargin<7, f_lims=[]; end
if nargin<8, ph_lims=[]; end
if nargin<9, iihann=1; end
if nargin<10, iilg=2; end
if nargin<11, cal_root=[]; end
%
% Check for existence of input data file:
[data_file,cts_root,iif,path_back]=adi_check_file(cts_root,'.cts', ...
	'DATA FILE');
if ~iif, disp('Exiting cts2fft ...'); return; end
%
if zpfac<0 | round(zpfac)~=zpfac | zpfac>10,
	disp(['zpfac = ',num2str(zpfac)]),
	error('USER ERROR: Bad zpfac: FFT size to be increased by 2^zpfac (0<=zpfac<=10)'); 
end
%
t0 = clock;
CTS=cts_open(cts_root);
data_type=CTS.header(7);  %tms:  data_type = 7 is ADAR, 8 = VLAD
%
nrec=CTS.nrec;
fs=CTS.fs;
%%if calfac<=0, error('USER ERROR: calibration factor must be > 0'); end
%
dt=1/fs; 
% For T_int<0, set T_int to get nT_int a power of 2:
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>=-100,
	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;
nT_step=round(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;
%
nfft=2 ^ ( zpfac + ceil(log2(T_int/dt)) );
nT_zero=nfft - nT_int;
T_fft=nfft/fs; df=fs/nfft;
%
ii_to=1:nT_int-nT_step;
ii_from=nT_step+1:nT_int;
ii_end=nT_int-nT_step+1:nT_int;
%
if isempty(ph_lims),
   iirec=[1:nrec];
else
   [nr,nc]=size(ph_lims);
   if nr==1 & nc==2,
      iirec=[ph_lims(1):ph_lims(2)];
   elseif nc==1,
      iirec=ph_lims';
   else
      disp('USER ERROR: ph_lims should be [jrmin,jrmax] or [jr1;jr2;...jrn]');
      return;
   end
end
nrec_out=length(iirec);     %disp(iirec),
if any(iirec<1) | any(iirec>nrec),
	disp(ph_lims); error('Illegal receiver specification in cts2fft');
	return
end
rec_indx_re=(iirec-1)*2 + 1; rec_indx_im=rec_indx_re + 1;
if iihann,
	hann=.5*(1 - cos((2*pi/T_int)*t_out)); 
% 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 window widens the mainlobe, reduces it by 1.76 dB, and lowers the sidelobes.
%%	hann_norm=sqrt(sum(hann.^2)/length(hann)),
	hann_norm=sqrt(3/8);
	hann_arr=(hann/hann_norm) * ones(1,nrec_out);
end
% Account for fact that we may want to output a frequency band of width greater than fs.
% f_lims_ok are those frequency limits inside f_bb_shift +- fs/2;
% f_lims can go beyond, and iif_ok are indices in cspec to place OK frequencies.
%f_lims_ok(1)=max(f_lims_all(1),f_lims(1));
%f_lims_ok(2)=min(f_lims_all(2),f_lims(2));
%nf1_ok=round(f_lims_ok(1)/df) + 1; nf2_ok=round(f_lims_ok(2)/df) + 1; 
%nf1=round(f_lims(1)/df) + 1; nf2=round(f_lims(2)/df) + 1;
%nfreq_ok=nf2_ok - nf1_ok + 1;
%iif_ok=nf1_ok - nf1 + [1:nfreq_ok];
%
%fmax=f_bb_shift + max(f_lims(2)-f_bb_shift,f_bb_shift-f_lims(1));   % Maximum frequency in CTS
%fs_new=2*fmax;
% Factor such that true fs covers basebanded band
%fac_bb=2^ceil(log2(fs_new/df))/nfft,
%
% Factor such that true fs covers basebanded band
f_bb_shift=CTS.f_bb_shift;
if isempty(f_lims),	f_lims=[f_bb_shift - fs/2 + df,f_bb_shift + fs/2]; end
fac_bb=2^ceil(log2(2*f_lims(2)/fs));
% Adjust basebanded nfft,fs such that true band is given:
SPEC.nfft=nfft*fac_bb; SPEC.fs=fs*fac_bb; SPEC.nT_int=nT_int*fac_bb;
SPEC.T_step=T_step; SPEC.iihann=iihann;
%
%
% Find frequency indices to output (should be very nearly an integer value):
f_bb_des=round(f_bb_shift/df)*df;
%%%f_bb_des=round(f_bb_shift/df)*df + 10*df;  % This was done to prove that sign below is correct 8-23-01
delf=f_bb_shift - f_bb_des; % OLD - NEW
dw=0; T0=0;
disp(sprintf('f_bb_shift = %1.5f, f_bb_des = %1.5f, delf = %1.5f',[f_bb_shift,f_bb_des,delf]));
if abs(delf) > 5e-5,
	disp('NOTE: f_bb_shift not a multiple of df!!!  Will multiply by exp(i*dw*t)');
% This agrees with bbdop_wf_prop_DB, fft_dtau_xx_init and has been proved correct using a sin wave:
	dw=2*pi*delf;
% This doesn't seem to agree with fft_dtau_bank_init!!  Try minus sign: DOESN'T WORK!!
%%    dw=2*pi*(f_bb_des - f_bb_shift);
%%   dw=-dw; disp('NEGATING dw to TEST!!');
else
    disp('NOTE: f_bb_shift very close to a multiple of df.  No exp(i*dw*dt) necessary');
end	
%
% Redo iif_out calculation a better way:
iif_out_all=[nfft/2+2:nfft,1:nfft/2+1];
fax_out_all_bb=[(-nfft/2)*df + [1:nfft]*df];
fax_out_all=fax_out_all_bb + f_bb_des;
ii=find(fax_out_all >= f_lims(1)-df/2 & fax_out_all <= f_lims(2)+df/2);
iif_out=iif_out_all(ii); 
fax_out=fax_out_all(ii); 
%% disp(iif_out([1:5,end/2-4:end/2+4,end-4:end])); disp(fax_out([1:5,end/2-4:end/2+4,end-4:end])); f_bb_des,

 %
nT_fac=2;	% Read in real and imaginary parts of complex time series
%
SPEC.Tw=T_fft;
% Compute imap for inverse FFT calculations:
%(fax_out(1)/df), (fax_out(end)/df),
nf1=round(fax_out(1)/df)+1; nf2=round(fax_out(end)/df)+1;
%
%%disp('TEMP nf1 adjust: '); nf1=nf1-10, nf2=nf2-10,
SPEC.nfreq=nf2-nf1+1; 
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_out; REC.zr=zeros(1,nrec_out);

%tms
if data_type==7  % ADAR
	f_ax=[nf1-1:nf2-1]'*df;
    f0_cal = (f_ax(1)+f_ax(end))/2;  % mid band should be 325
    dB_octave=-6;
    calfac = 10^(82.31/20);  % ?? Ron said -82.31;  he must be talking dB
    freq_cal=zeros(SPEC.nfreq,nrec_out);
    for jrec=1:nrec_out,
        dB_cal=log2(max(.1,f_ax/f0_cal)) * dB_octave;
        freq_cal(:,jrec)=10.^(dB_cal/20);       %figure; plot(f_ax,freq_cal);
    end
    freq_cal = freq_cal*calfac;
elseif data_type==8 % VLAD
%     iidepth=CTS.header(7+data_type);
	f_ax=[nf1-1:nf2-1]'*df;
    log_fax=log2(f_ax);
	if ~isempty(cal_root), 
        eval(['load ',cal_root,'.mat']);
        if CAL.dims(2) ~= nrec_out,
            error(['CAL File data should have nrec calibration numbers: ',int2str([CAL.dims(2),nrec_out])]);
        end
	end
    freq_cal=ones(SPEC.nfreq,nrec_out);
    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_out,
        dB_cal=interp1(log2(CAL.scales{1}),cal_file_fac(:,jrec),log_fax,'linear','extrap');
        freq_cal(:,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(:,jrec)=freq_cal(:,jrec) .* exp(i*phase_cal*pi/180);
        end
    end
end
%
% Open output FFT file:
if isempty(fft_root),
	fft_root=[cts_root,'_',int2str(T_int),'_',int2str(T_step)]; 
end
FFT=fft_write1_open(fft_root,SPEC,REC);
%%	fft_norm=.5*sqrt(2/(nT_int*nfft)); IS THIS RIGHT???
% This is standard norm used in processing real time series:
fft_norm=sqrt(2/(nT_int*nfft));		% THIS IS HOW IT WAS ORIGINALLY
% For real time series we use sqrt(2) factor to account for negative frequencies
% that are thrown away.  Factor not needed here since all energy is represented in
% complex time series and whole spectrum is shifted to positive frequencies: 
fft_norm=fft_norm/sqrt(2);			
%%fft_norm=sqrt(2/(FFT.nT_int*FFT.nfft));
% Position file to starting point:
nT_t1=round(3600*(t_lims(1)-CTS.Jhour1)/dt);
fmt='real*4'; nbyte_sample=8;
nbyte=nrec*nbyte_sample*nT_t1;		% nbyte_sample bytes per sample
if nbyte>0,
	fseek(CTS.fid_cts,nbyte,0);		% Seek from current position since cts file has header
else
	t_lims(1)=CTS.Jhour1;
end
%%'t_lims,CTS.T_lims: ', disp([t_lims,CTS.T_lims]);
% Adjust last requested time so that projected CPU times are valid:
t_lims(2)=min(t_lims(2),CTS.T_lims(2));
t_hr=t_lims(1);
% Read in first FFT's worth for all receivers:
npt_read=2*nrec*nT_int;
[pser,count]=fread(CTS.fid_cts,[2*nrec,nT_int],fmt);
if count<npt_read,
%    size(pser),size(zeros(2*nrec,(npt_read-count)/(2*nrec))),
%    keyboard
    pser=[pser, zeros(2*nrec,(npt_read-count)/(2*nrec))];
    count=npt_read;
    disp('Zero filled to make first FFT');
end
pser=pser(rec_indx_re,:)' + i*pser(rec_indx_im,:)';		% pser is nT_int by nrec_out complex;
if dw~=0,
% Change basebanding frequency from f_bb_shift to f_bb_des:
	pser=pser.*repmat(exp(i*dw*(T0 + [0:nT_int-1]'*dt)),[1,nrec_out]);
	T0=T0 + nT_int*dt;
end
nT_read=count/(2*nrec);
%
ptot=zeros(nT_int,nrec_out); pshade=zeros(nfft,nrec_out);
cspec=zeros(SPEC.nfreq,nrec_out);
% 
done=0;
if nT_read<nT_int,
	disp('Warning: CTS file ran out of data on first FFT: nT_int,nT_read,nrec = ');
	disp([nT_int,nT_read,nrec]);
	ptot(iit(1:nT_read),:)=pser;
%%	done=1;
	LGRAM=[]; return
else
	ptot(iit,:)=pser;
end
%
%ptot(iit,:)=calfac*pser;
%
nseq=0;
% Phase factor associated with basebanding:
exp_fac=1;
% This works now (10-30-00):
%%exp_const=exp(i*2*pi*f_bb_shift*T_step), 
% FIX(??) 2-8-02:
%%disp('Setting exp_const according for f_bb_des rather than CTS.f_bb_shift');
% Keeping the old way raised sidelobes on correlation considerably (using 12-s FFTs).
% This is correct:
exp_const=exp_const_calc(f_bb_des,T_step,df,dt,1); 
%%%%angle(exp(i*2*pi*f_bb_des*T_step))*180/pi,
% 2-12-02: in ABX, it was found that conjugating exp_hann_const made simulated time series smooth,
% but data time series jump.  Problem must be here, where wrong sign being applied:
%%disp('NOTE: sign of exp_const NOT switched');
%%exp_const=exp(-i*2*pi*f_bb_des*T_step); disp('NOTE: sign of exp_const switched');
disp(sprintf('exp_const angle = %1.5f',angle(exp_const)*180/pi)); %% disp(T_step), disp(f_bb_des),
while 1==1,
	nseq=nseq + 1;
	if iihann,
		pshade(iit,:)=ptot(iit,:).*hann_arr;
	else
		pshade(iit,:)=ptot(iit,:);
	end
	spec=fft(pshade,nfft,1); 
	cspec(:,:)=(exp_fac*fft_norm)*spec(iif_out,:);
    if any(data_type == [7 8])
        cspec=cspec.*freq_cal; 
    end
	exp_fac=exp_fac*exp_const;
%
	t_start=t_hr*3600;		% t_start in seconds of Julian day Jday;
	stat=fft_write1(FFT,cspec,nseq,[],CTS.Jday,t_start);
%
%%	if done,
%%		disp(['Stopping because end of file reached: ', ...
%%			sprintf('%3.3f,%3.3f',[t_hr,t_lims(2)])]);
%%		break;
%%	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
% Read in next chunk of data:
	[pser,count]=fread(CTS.fid_cts,[2*nrec,nT_step],fmt);
	nT_read=count/(2*nrec);
	if nT_read<nT_step,
		disp(['Stopping because end of file reached: ', ...
			sprintf('%3.3f,%3.3f',[t_hr,t_lims(2)])]);
		break;
	end
	pser=pser(rec_indx_re,:)' + i*pser(rec_indx_im,:)';		% pser is nT_step by nrec_out complex
	if dw~=0,
% Change basebanding frequency from f_bb_shift to f_bb_des:
		pser=pser.*repmat(exp(i*dw*(T0 + [0:nT_step-1]'*dt)),[1,nrec_out]);
		T0=T0 + nT_step*dt;
	end
	ptot(ii_to,:)=ptot(ii_from,:);
	ptot(ii_end,:)=pser;
	t_hr=t_hr + T_step/3600;
	if rem(nseq,10)==1,
		ctime=etime(clock,t0)/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
%
fclose(CTS.fid_cts);
fclose(FFT.fft_fid);
disp(['# Sequences output: ',int2str(nseq)]);
ctime=etime(clock,t0)/60;
disp(['CPU Time = ',sprintf('%2.2f',ctime), ...
	' min. Output file = ',fft_root,'_fft']);
%
if iilg
	if iilg>1, iifig=1; else, iifig=0; end
%	[LGRAM,FFT]=lgram_3d(fft_root,t_lims,[],f_lims,1,1,iifig);
	[LGRAM,FFT]=lgram_3d(fft_root,t_lims,[],[],[1 1 0 1],iifig,[]);
else
	LGRAM=[];
end
%
%%%
%
function exp_const=exp_const_calc(f_bb,T_step,df,dt,sg);
%
% Computes exp(i*2*pi*f_bb*tscan) using facts that f_bb is an integral
% multiple of df and tscan is an integral multiple of dt, and that
% NFFT=1/(df*dt).
%
M_df=round(f_bb/df); N_dt=round(T_step/dt); NFFT=round(1/(df*dt));
disp([M_df,N_dt,NFFT]);
while mod(M_df,2)==0 & mod(NFFT,2)==0,
    M_df=M_df/2; NFFT=NFFT/2;
end
disp([M_df,N_dt,NFFT]);
while mod(N_dt,2)==0 & mod(NFFT,2)==0,
    N_dt=N_dt/2; NFFT=NFFT/2;
end
disp([M_df,N_dt,NFFT]);
dfac=mod(M_df/NFFT*N_dt,1);
phase=sg*2*pi*dfac; phase*180/pi;
exp_const=cos(phase) + i*sin(phase);
%
% 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.
