function [backgrounds, IndexNB] = JCSD_BackgroundExtractor(data,wavel)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% General description:
%
% This function extract background information from a set of data.
%
% Idea:
% When background is the majority of the observation data, there will be a
% majority of elements in autocorrelation matrix  close to 1 if the
% noise level is not too high. This function find a maximum clique in the
% correlation matrix, and the background can be approximated by the average
% of background components.
%
% input:
%       data - observation data
%       wavel - wavenumber 
%
% output: 
%       backgrounds - average background component
%       IndexNB - Indices of non-background data
%
% Written by: Wei Wang
% Machine Learning for Signal Processing Lab
% University of Maryland, Baltimore County
% 9/22/2007
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%       threshold - correlation value greater than 'threshold' is
%                   regarded as same component. Set to =0.9 in this
%                   function, or can be determined by a ratio of all corrs.

threshold1=0.9;  % threshold for being regarded as same spectrum
bkg_density=0.7; % a estimate ratio of background in observations, used to 
                 % determine threshold for "similarity"

[len,num]=size(data);

blocksize=500; %
blocksize=num;
blocksize=min(blocksize,num);
m=ceil(num/blocksize);
bkgs=zeros(len,m);


% determine a threshold from give background density
x1=data(:,1:blocksize);
cors=corrcoef(x1);
cor_uptri=triu(cors,1);
[ii,jj,cor1s]=find(cor_uptri);
[vcor,ids]=sort(cor1s,'descend');
pos=round(bkg_density*length(ids)); 
threshold2=vcor(pos);
threshold=max(threshold1, threshold2);

% Find the first background
G=cors - eye(blocksize)> threshold;
[MaxClique] = MaxClique_maxDegree(G);
bkg1st=mean(x1(:,MaxClique),2);

rrs=corr(bkg1st,data);
idxBkg=find(rrs>=threshold);
avgBkg1st=mean(x1(:,idxBkg),2);



% Extract non-background data to find the 2nd background
idxNB=[1:blocksize];
idxNB(idxBkg)=[];
dataNB=data(:,idxNB);
nn=size(dataNB,2);
cors2 = corrcoef(dataNB);

G=cors2 - eye(nn)> threshold;
[MaxClique2] = MaxClique_maxDegree(G);
bkg2nd=mean(dataNB(:,MaxClique2),2);
rrs2=corr(bkg2nd,dataNB);
idxBkg2nd=find(rrs2>=threshold);
avgBkg2nd=mean(dataNB(:,idxBkg2nd),2); % the second background

%
% Compare the 1st and 2nd backgound, if their peaks are not matching, then
% only keep the first background as the true background, and treat the 
% second background as non-background data
%
[posBkg1, ampsBkg1, yBkg1, y_smoothed]=extracting_peaks(avgBkg1st',wavel);   
[posBkg2, ampsBkg2, yBkg2, y_smoothed]=extracting_peaks(avgBkg2nd',wavel);
locsBkg1=wavel(posBkg1)';
locsBkg2=wavel(posBkg2)';
dt=20;
[locsData,idA,idB, ampsBkg1_all, ampsBkg2_all] ...
    = union_peaks(locsBkg1, locsBkg2, ampsBkg1, ampsBkg2, dt);
diffPeaks=setdiff(idB,idA);
if length(diffPeaks)> 0
    ampsDiff=ampsBkg2_all(diffPeaks);
    weightDiff=norm(ampsDiff)^2/norm(ampsBkg2_all)^2;
    if weightDiff > 0.2    
        backgrounds=avgBkg1st;
        IndexNB=idxNB;
        return;
    end
end

% output
backgrounds=[avgBkg1st avgBkg2nd];

IndexNB=idxNB;
IndexNB(idxBkg2nd)=[];



