function [mcs_bf,cbf_bf]=mcs_bf_calc(jf_bf,xspec)
%
% xspec should be nfbin by nph
% xspec should be nf_bf by nph
%
% 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.
mcs_global
%
%%jfbin=BEAMS.iifbin(jf_bf);
jf_rep=BEAMS.jf_bin2rep(jf_bf);
if jf_rep<1 | jf_rep > BEAMS.nf_rep,
	disp('ERROR: jf_rep bad in mcs_bf_calc: '); jf_bf,jfbin,jf_rep, keyboard;
end
nbeam=REP(jf_rep).nbeam;
% Find pruning frequency index (into PRUNE):
jfp=REP(jf_rep).jf_prune; nph=PRUNE(jfp).nph;
%%int2str([jf_bf,jf_rep,jfp,nph,nbeam]),
% Extract spectra of phones used at this frequency:
%%size(xspec),jfbin,jfp,PRUNE(jfp).iiph,
X=xspec(jf_bf,PRUNE(jfp).iiph).';	% Transpose makes this nph by 1;
% Compute beamformed weights and sqared magnitudes element-by-element:
bf_elements=REP(jf_rep).wtsvec .* (X*ones(1,nbeam));	% nph by nbeam;
%%bf_elements=conj(REP(jf_rep).wtsvec) .* (X*ones(1,nbeam));	% nph by nbeam;
bf_norm=real(bf_elements).^2 + imag(bf_elements).^2;
% Form cumulative sums:
bf_cum=[zeros(1,nbeam); cumsum( bf_elements,1 )];	% nph+1 by nbeam;
norm_cum=[zeros(1,nbeam); cumsum( bf_norm,1 )];		% nph+1 by nbeam;
kph1=PRUNE(jfp).kph1; kph2=PRUNE(jfp).kph2;
nsub1=PRUNE(jfp).nsub1*ones(1,nbeam); nsub2=PRUNE(jfp).nsub2*ones(1,nbeam); 
% Evaluate subaperture beamformer outputs and norms from cum sums:
% Weight normalizations by number of phones in other subarray:
bf_xc1=bf_cum(kph1(2,:)+1,:) - bf_cum(kph1(1,:),:);	% n_xc by nbeam;	
bf_xc2=bf_cum(kph2(2,:)+1,:) - bf_cum(kph2(1,:),:);
xc_norm1=nsub2 .* (norm_cum(kph1(2,:)+1,:) - norm_cum(kph1(1,:),:));	
xc_norm2=nsub1 .* (norm_cum(kph2(2,:)+1,:) - norm_cum(kph2(1,:),:));	
% Compute correlations:
% (We could also sum the normalizations before dividing by them)
xc_num=real(bf_xc1).*real(bf_xc2) + imag(bf_xc1).*imag(bf_xc2);	% Re(bf1 * bf2')
%%keyboard,
if MCS.iinm==0,
	xc_all=xc_num ./ sqrt(xc_norm1 .* xc_norm2);
	xc_all_wt=xc_all .* ( PRUNE(jfp).nph_prod' * ones(1,nbeam) );
% Sum crosscorrelations and weight according to # phones:
	mcs_bf0=sum(xc_all_wt,1) / PRUNE(jfp).nph_prod_tot;		% 1 by nbeam;
else
	xc_norm=sqrt(sum(xc_norm1,1) .* sum(xc_norm2,1));
	xc_all_wt=xc_num .* ( PRUNE(jfp).nph_prod' * ones(1,nbeam) );
% Sum crosscorrelations and weight according to # phones:
	mcs_bf0=sum(xc_all_wt,1) ./ (xc_norm * PRUNE(jfp).nph_prod_tot);		% 1 by nbeam;
end
%
% CBF:
cbf=bf_cum(nph+1,:);
cbf_bf0=real(cbf).^2 + imag(cbf).^2;
%
% If variable number of beams vs. frequency, apply quadratic-fit weights:
if BEAMS.iivar,
	bf_wts=mcs_bf0(REP(jf_rep).idx_st) .* REP(jf_rep).qfac_st;
	mcs_bf=sum( bf_wts , 1);
	bf_wts=cbf_bf0(REP(jf_rep).idx_st) .* REP(jf_rep).qfac_st;
	cbf_bf=sum( bf_wts , 1);
%%	figure; plot(BEAMS.phi_out,mcs_bf,REP(jf_rep).phi_bf,mcs_bf0); keyboard,
else
	mcs_bf=mcs_bf0; cbf_bf=cbf_bf0;
end
