function [SGRAM,FFT]=sgram_3d(fft_root,Tlims,jrec,f_lims,tau_lims,parms, ...
	iifig,clim);
%
% Call [SGRAM,FFT]=sgram_3d(fft_root,Tlims,jrec,f_lims,tau_lims,parms, ...
%	iifig,clim);
%
% fft_root	= Root of FFT file (suffix = _fft) ([] to be prompted);
% Tlims		= FFT sequences to output ([] for all);  Column vector of sequences, or
%			  [Tmin,Tmax], where T in hod for iihr=1, or minutes into FFT file 
%			  for iihr=0, or range in km for iihr=3;
% jrec		= List of receiver numbers to output;
% f_lims	= [fmin,fmax], or column vector of frequency indices, to use in IFFT;
% tau_lims	= [tau_min,tau_max] time interval to output ([] for all) in seconds;
% parms		= [iienv,fzpad,cref], where:  
% 	dfnorm	= 1 for data FFT, 0 for simulation, where range is used on vertical axis;
% 	iihr	= vertical axis units: 0 = min into FFT, 1 = HOD, 2 = range in km;
%	iienv	= flag to envelope (0=no;1=|env|;2=dB(|env|);3=arg(env)/pi;4=unwrap(arg(env))/pi);
%	fzpad 	> 0 to zero pad by 2^fzpad (fzpad<0 ==> reduce fs),
%			< 0 to reduce sample rate and FFT size (when enveloping);
% 	iihann	= 1 to put a Hann window in the frequency domain before inverse FFTing (def=1);
%	iinorm	= 1 to normalize each time series to unity (0 dB);
%	ns_avg	= number of sequences to average before outputing (default = 1);
%	cref	= Ref sound speed to align arrivals (0 for none) (not implemented yet);
% iifig		= 1 to call view3(SGRAM);
% clim		= [cmin,cmax] for colorbar.
%
% OUTPUT is structure SGRAM.  Use view5 to produce image or line plots.
% SGRAM.dims		= [time, sequence or time, receiver or beam number];
%
%
% 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 nargin<1, fft_root=[]; end
if nargin<2, Tlims=[]; end
if nargin<3, jrec=[]; end
if nargin<4, f_lims=[]; end
if nargin<5, tau_lims=[]; end
if nargin<6 | isempty(parms), parms=[1 1 2 -3 1 0]; end
if nargin<7, iifig=1; end
if nargin<8, clim=[]; end
dfnorm=parms(1); iihr=parms(2); iienv=parms(3); fzpad=parms(4); 
if length(parms)>=5, iihann=parms(5); else, iihann=1; end
if length(parms)>=6, iinorm=parms(6); else, iinorm=0; end
if length(parms)>=7, ns_avg=parms(7); else, ns_avg=1; end
if length(parms)>=8, f_bb_shift=parms(8); else, f_bb_shift=0; end
fft_file=[fft_root,'_fft'];
if exist(fft_file)~=2,
	disp(['Note: ',fft_file,' not found locally.']);
	[fft_file,path_name]=uigetfile('*_fft','Please Choose _fft File');
	if fft_file==0, SGRAM=[]; return; end
	lfft=length(fft_file); fft_root=fft_file(1:lfft-4);
	if ~strcmp(path_name,pwd), eval(['cd ',path_name]); end
end
%
FFT=fft_read1_open(fft_root,dfnorm);
fax=(FFT.imap(:,1) - 1) * FFT.df;
[nr,nc]=size(f_lims);
if nr==1 & nc==2,
	jfreq=find(fax>=f_lims(1)-FFT.df/2 & fax<=f_lims(2)+FFT.df/2);
elseif isempty(f_lims), 
	jfreq=[1:FFT.nfreq];
end
if isempty(jfreq),
	disp('No overlap between desired frequency band and FFT band: ');
	f_lims,FFT.f1,FFT.f2,
	keyboard
end
nfreq=length(jfreq); fax_get=fax(jfreq); 
if isempty(jrec), jrec=[1:FFT.nrec]; end
%
%
[nr,nc]=size(Tlims);
if isempty(Tlims), 
	jseq=[1:FFT.nseq];
elseif nr==1 & nc==2,
% convert Tlims in min after start of FFt file to hod;
	if iihr==0, Tlims=mod(FFT.moy0+Tlims,1440)/60; end
	if iihr==2,
		jseq=find(FFT.range>=Tlims(1) & FFT.range<=Tlims(2));
	else	
		jseq=find(FFT.hodax>=Tlims(1) & FFT.hodax<=Tlims(2));
	end
else,
	jseq=Tlims(:)';
end
if isempty(jseq),
	Tlims,FFT.hodax(1),FFT.hodax(length(FFT.hodax)),
	disp('USER ERROR: No FFT sequences to be processed.'); keyboard; return;
end
%
SGRAM.title=fft_root;
%
if ns_avg<1 | round(ns_avg)~=ns_avg,
	disp(['USER ERROR: ns_avg must be integer>=1: ',int2str(ns_avg)]);
	return
end
%
nrec=length(jrec); nseq=length(jseq); jseq_out=jseq(ns_avg:ns_avg:nseq);
nseq_out=floor(nseq/ns_avg);
%
if iienv==0,
    fac=2^fzpad;
else
    if fzpad<0, disp('fzpad should now be >=0 for enveloping also'); fzpad=0; end
    fband=diff(fax_get([1,end]));
    fac=2^(-(nextpow2(FFT.fs/fband)-1-fzpad)); fsx=FFT.fs*fac; %%disp(['fsx = ',num2str(fsx)]);
end
nfftx=round(FFT.nfft*fac); fsx=FFT.fs*fac; nT_intx=round(fac*FFT.nT_int);
if isempty(tau_lims), tau_lims=[0,FFT.Tfft-1/fsx]; end
tau_lims(1)=max(tau_lims(1),0); tau_lims(2)=min(tau_lims(2),FFT.Tfft-1/fsx);
jt=round(tau_lims*fsx)+1; iit=jt(1):jt(2); ntau=length(iit);
t=[iit-1]'*(1/fsx);
%
SGRAM.ndim=3; SGRAM.dims=[ntau,nseq_out,nrec];
SGRAM.labels={'Relative Time - s' [] 'Receiver Number'};
SGRAM.scales={t [] jrec};
SGRAM.data=zeros(ntau,nseq_out,nrec);
if iihr==2,
	SGRAM.scales{2}=FFT.range(jseq);
	SGRAM.labels{2}='Range - km'
elseif iihr==1,
	SGRAM.scales{2}=FFT.hodax(jseq_out);
	SGRAM.labels{2}=['Time - hour of day ',int2str(FFT.jday0)];
else
	SGRAM.scales{2}=FFT.jday0*1440 + FFT.hodax(jseq_out)*60 - FFT.moy0;
	SGRAM.labels{2}=['Time - min'];
end
SGRAM.dlab='Pressure - \muPa';
%
if iihann==1,
	delf=fax_get(nfreq) - fax_get(1);
	hann=.5*(1 - cos((2*pi/delf)*(fax_get-fax_get(1)))); 
% Hanning window widens the mainlobe, reduces it by 1.76 dB, and 
% lowers the sidelobes.
	hann_norm=sqrt(sum(hann.^2)/length(hann));
% hann_arr is nfreq by nrec:
	hann_arr=repmat((hann/hann_norm),[1,nrec]);
elseif iihann>0,
	f0=round(iihann); w=mod(iihann*10000,10000); 
	disp(['Gaussian f0,w = ',int2str([f0,w])]);
	hann=exp(-((fax_get-f0)/w).^2);
	hann_norm=sqrt(sum(hann.^2)/length(hann));
% hann_arr is nfreq by nrec:
	hann_arr=repmat((hann/hann_norm),[1,nrec]);
end
%
% In case either of tau_lims < 0:
iineg=find(iit<1); iit(iineg)=nfftx + iit(iineg);
fft_full=zeros(nfftx,nrec) + i*zeros(nfftx,nrec);
if iienv,
% Shift band so it is centered on f=0:
%%	iih=ceil(nfreq/2); iix=nfftx-nfreq+iih + 1;
%%	pser=fft_full; x=zeros(ntau,nrec) + i*ones(ntau,nrec);
%
    if f_bb_shift==0, f_bb_shift=fax_get(ceil(nfreq/2)); end
    jf_bb=round((f_bb_shift - fax_get(1))/FFT.df) + 1;
    f_bb_shift=fax_get(jf_bb);
	if jf_bb<=0 | jf_bb>length(fax_get),
		error(['Desired baseband frequency not in FFT file band.  Aborting.']);
	end
	if abs(f_bb_shift-fax_get(jf_bb)) > .00001,
		disp(['NOTE: basebanding bin not equal to desired f_bb: ', ...
		num2str([f_bb_shift,fax_get(jf_bb),FFT.df,f_bb_shift-fax_get(jf_bb)])])
	end
	pser=fft_full; x=zeros(ntau,nrec) + i*ones(ntau,nrec);
	iipos=[jf_bb:nfreq]; npos=length(iipos); iipos_bb=[1:npos];
	iineg=[1:jf_bb-1]; nneg=length(iineg); iineg_bb=[nfftx-nneg+1:nfftx];
else
	imap1=FFT.imap(jfreq,1); imap2=nfftx + 2 - imap1;
	pser=zeros(nfftx,nrec); x=zeros(ntau,nrec); f_bb_shift=0;
end
%
% First term is to reverse normalization applied to forward FFTs,
% second term numerator reverses norm of this ifft, and second term
% denominator applies (1/N) correct inverse FFT that would be applied
% if no change in FFT order were done:
fft_norm=sqrt(FFT.nfft*FFT.nT_int/2)*(nfftx/FFT.nfft);
% In case of enveloping, account for negative frequencies that are filtered out:
if iienv, fft_norm=sqrt(2)*fft_norm; end
%%if dfnorm, fft_norm=fft_norm/sqrt(FFT.df); end
%
% The integral of spectral level over frequency of the FFT (LGRAM) 
% must be the same as the average squared pressure level over time:
% For continuous (non-transient) data signals:
% zrms(SGRAM.data,SGRAM.iidB) = zint(LGRAM.data,LGRAM.iidB,diff(LGRAM.scales{1}(1:2)))
%
% For TL & Transfer functions, want sum(SG(t)) * dt = sum(LG(t) / nfreq
% so that, for one arrival, max(SG(t)) = avg(TL(f)).
% Call lgram_3d with dfnorm=0, rgram_3d with dfnorm=1 (doesn't matter), and mult SG.data by the
% factor dt/nfft.
%
% For iienv=0,1,-1, SGRAM.data are real or complex time series in linear units, 
% in which case,
% 10*log10(mean(zmagsq(SGRAM.data),1))=
%	10*log10(sum(10.^(LGRAM.data/10))*diff(LGRAM.y(1:2)))
% For iienv=2, SGRAM.data are in dB, in which case,
% 10*log10(mean(10.^(SGRAM.data/10),1))=
%	10*log10(sum(10.^(LGRAM.data/10))*diff(LGRAM.y(1:2)))
%
ns_out=0; 
for js=1:nseq,
	if js==1,
		nrec_out=nrec;
		if ns_avg>1, cspec_avg=zeros(nfreq,nrec); end
		cspec=zeros(nfreq,nrec);
	end
	js_out=jseq(js);
	[XH,cspec(:,:)]=fft_read1(FFT,js_out,jrec,jfreq);
% Must hann right away for the case iinorm;
% Otherwise, correlations won't be normalized:
	if iihann, cspec=cspec.*hann_arr; end
% Average over sequence number (meant for averaging cross_spectra):
	if ns_avg>1,
		cspec_avg(:,:)=cspec_avg(:,:) + cspec;
		if mod(js,ns_avg)==0, cspec(:,:)=cspec_avg/ns_avg; cspec_avg(:,:)=0; end
	end
	if mod(js,ns_avg)==0,
		ns_out=ns_out + 1;
% Do inverse FFTs:
		if iienv==0,
        	fft_full(imap1,:)=fft_norm*cspec;
			fft_full(imap2,:)=conj(fft_full(imap1,:));
			pser(:,:)=ifft(fft_full,nfftx,1);
			x(:,:)=real(pser(iit,:));
		else
%			fft_full(1:nfreq-iih,:)=fft_norm*cspec(iih+1:nfreq,:);
%			fft_full(nfftx-iih+1:nfftx,:)=fft_norm*cspec(1:iih,:);
%			pser(:,:)=ifft(fft_full,nfftx,1);
%			x(:,:)=pser(iit,:);            
            fft_full(iipos_bb,:)=fft_norm*cspec(iipos,:);
			fft_full(iineg_bb,:)=fft_norm*cspec(iineg,:);
			pser(:,:)=ifft(fft_full,nfftx,1);	% nfftx by nrec;
			x(:,:)=pser(iit,:);					% ntau by nrec;
		end
		if iienv==0,
%%			x(:,:)=real(x);		% Already done above;
		elseif iienv==1,
			x(:,:)=abs(x);
		elseif iienv==2,
			x(:,:)=10*log10(max(1e-30,zmagsq(x)));
		elseif iienv==3,
			x(:,:)=angle(x)/pi;
		elseif iienv==4,
			x(:,:)=unwrap(angle(x))/pi;
		end
		SGRAM.data(:,ns_out,:)=x(:,:);
%%		if ~isdata, SGRAM.t_ref(ns_out)=XH.t0; end
		if iihr==1,
			SGRAM.scales{2}(ns_out)=XH.hod;
		else iihr==0,
			SGRAM.scales{2}(ns_out)=XH.moy+XH.t0/60 - FFT.moy0;
		end
	end
%
%%	[t,pser]=myffti(cspec,FFT,tau_lims,0,fzpad,iienv);
	if mod(js,10)==1, disp(['Finished sequence ',int2str(js),' of ',int2str(nseq)]); pause(.001); end
end
%
if iinorm,
	if iienv==2,
		pmax=max(SGRAM.data,[],1);
		SGRAM.data=SGRAM.data - repmat(pmax,[ntau,1,1]);
	elseif iienv<2,
		pmax=max(abs(SGRAM.data),[],1);
		SGRAM.data=SGRAM.data./repmat(pmax,[ntau,1,1]);
	end
end
%
if isempty(clim),
	cmax=max(SGRAM.data(:));
	if iienv==0
		cmin=min(SGRAM.data(:));
	elseif iienv==1,
		cmin=0;
	elseif iienv==2,
		cmax=ceil(cmax/5)*5;
		cmin=max(floor(min(SGRAM.data(:))/5)*5,cmax-60);
		SGRAM.dlab='Pressure - dB//\muPa^2';
	elseif iienv==3,
		cmin=-1; cmax=1;
		SGRAM.dlab='arg(env)/pi';
	elseif iienv==4,
		cmin=min(SGRAM.data(:));
		SGRAM.dlab='unwrap[arg(env)]/pi';
	elseif iienv==-1,
		cmax=max(dB(SGRAM.data(:)));
		cmax=ceil(cmax/5)*5;
		cmin=max(floor(min(dB(SGRAM.data(:)))/5)*5,cmax-60);
		SGRAM.dlab='Complex pressure - \muPa';
	end
	clim=[cmin,cmax];
elseif iinorm,
	if iienv==2,
		clim=[clim(1)-clim(2),0];
	end
end
SGRAM.clim=clim;
if iienv==2, SGRAM.iidB=1; elseif iienv==-1, SGRAM.iidB=-1; else, SGRAM.iidB=0; end
%
SGRAM.fzpad=fzpad; SGRAM.iienv=iienv;
%
fclose(FFT.fft_fid);
%
if iifig & ~isempty(SGRAM.data),
	view5(SGRAM); 
end
