function [threshold, mu]=JCSD_threshold_background_DCCP(dataBkg, libs, clusters, blocksize, confidence)
% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% General description:
% This function determines threshold from background data for DCCP detector
%
% input:
%       dataBkg - background data
%       libs - spectrum library
%       clusters - spectrum library partitions
%       blocksize - size of data block used in detection. Usually set as
%                   the number of components. e.g., blocksize = 2 if there
%                   are 1 background and 1 possible target 
%       confidence - the possibility that samples are less than the
%                    threshold with Gaussian fit. e.g., confidence = 0.9
%
% output: 
%       threshold - see confidence above
%
% Written by: Wei Wang
% Machine Learning for Signal Processing Lab
% University of Maryland, Baltimore County
% 8/25/2007
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[n,m]=size(dataBkg);
mm=ceil(m/blocksize);
ccs=zeros(1,mm);
numclusters=size(clusters,1);

for j=1:numclusters
    icount(j)=length(find(clusters(j,:)>0));
end

for i=1:mm    
    istart=(i-1)*blocksize+1;
    iend=min(i*blocksize,m);
    x=dataBkg(:,istart:iend);

    % Remove sigularity
    while rank(x)~=size(x,2)
        x=x(:,1:end-1);
    end
    
    % Both are column vectors
    for j=1:numclusters
        cp=libs(:,clusters(j,1:icount(j)));     
        [A,B,r1] = canoncorr(x,cp);
        rr(j)=r1(1);
    end
    ccs(i)=max(rr);   
end
% 
% [mu,sigma]=normfit(ccs);
% threshold=norminv(confidence, mu, sigma);


mu=mean(ccs);
sbkg=sort(ccs,'ascend');
threshold=sbkg(round(length(sbkg)*confidence));




