
function [HIT, FA, HITFA, ACC]= bmdnnte(TefPath,TmPath,parameters,DNNPara)
%% This function implements the system described in JIANG and Wang,
% "Binaural Classification for Reverberant Speech Segregation Using Deep Neural Networks," submitted to IEEE
% TASLP, 2014.
%
% Inputs:
%    TfPath             binaural wav train file path;include the '.wav' 
%    TmPath             train files' IBM path, include the '.mat' 
%   parameters          optional parameter settings.  note that parameter settings must match those used for model training.  default parameters match those used for the provided models.
%       .fn            amount train files (default the number of '.wav' file in the train file path )  
%       .fs             sample rate (default = 16 kHz)
%       .nc             number of frequency channels (default = 64)
%       .fRange         lower and upper center frequencies (default = [80 5000] Hz)
%       .winL        	window length in s (default = 0.02)
%       .step           step size in s (default = 0.01)
%       .maxitd         maximum ITD value in ms (default = 1)
%       .maxiid         IID dimensions (default = 2)
%     DNNPara           trained DNN parameters 
% Output:
%    HIT: the hit rat of the IBM one
%    FA:  the false alarm of the IBM zero
%    HITFA: hit-fa
%    ACC: the estimation accuracy of the IBM
%
% Written by Yi Jiang (2013)
%%
addpath(genpath(pwd));
TfLIST=dir(fullfile(TefPath,'*.wav'));
TmLIST=dir(fullfile(TmPath,'*.mat'));
%Runing parameter
if nargin < 2
    parameters = struct;
end
% Parameter configuration
if ~isfield(parameters,'fn')
    parameters.fn = size(fLIST,1);      % tes files number
end
if ~isfield(parameters,'fs')
    parameters.fs = 16000;              % sample rate
end
if ~isfield(parameters,'nc')
    parameters.nc = 64;                 % number of filter channels
end
if ~isfield(parameters,'fRange')
    parameters.fRange = [50 8000];      % low- and high- center frequencies for filterbank
end
if ~isfield(parameters,'winL')
    parameters.winL = 0.020;            % window length (in s)
end
if ~isfield(parameters,'step')
    parameters.step = 0.010;            % step size (in s)
end
if ~isfield(parameters,'maxitd')
    parameters.maxitd = 1;              % maximum ITD value (in ms)
end
if ~isfield(parameters,'maxiid')
    parameters.maxiid = 2;              % maximum IID size (in groups)
end

%%set the test environments
Analyzecorpus=TefPath;
trainDNN=DNNPara;
channelnumber=parameters.nc ;
fs= parameters.fs;                      % sampling rate
fRange = parameters.fRange;             % default frequency range in Hz
FranmeT=parameters.winL;
winLength = FranmeT*fs;                 % default window length in samp
winShift =parameters.step*fs;           % 0.5*winLength;frame shift (default is half frame)

egroupsize=floor(winLength/parameters.maxiid);  %IID feature size
tdnumber=floor(parameters.fs/1000*parameters.maxitd)+1; %cross correlation feature size
tdlength= winLength-tdnumber;
unitsize= winLength/egroupsize;%units numberh
groupsize=1;


outputpreda=[];
testay=[];
outenc=0.0;
outenn=0.0;
tarenc=0.0;
sigx1=[];
sigx2=[];
outputpreda=[];
framen=0;
fsM=[];

matlabpool local;

%% test part binary mask estimated
testfn=parameters.fn
testy=[];
testya=[];
fifran=[];
parfor tn1=1:testfn
[tgsig,fs,nbits]=wavread([Analyzecorpus TfLIST(tn1).name]);%load the target sound
  gfx1= gammatone(tgsig(:,1),channelnumber, fRange,fs);   
  gfx2= gammatone(tgsig(:,2),channelnumber, fRange,fs); 
  sigLength= size(gfx1,2);
  M = floor((sigLength-winLength+winShift)/winShift);     % number of time frames
  sigLength=winLength+(M-1)*winShift;
  
  sigx1=[sigx1 gfx1(:,1:sigLength)];
  sigx2=[sigx2 gfx2(:,1:sigLength)];
  framen=framen+M;
  fifran=[fifran M];% each file frame.  
  masklab=load([Analyzecorpus TmLIST(tn1).name]);
  mclab=masklab.mask';
  testy=[testy;mclab];
end
emaska=[];
outputpreda=[];
outacc=[];
outfa=[];
outhit=[];
%% DNN classifuication by channel
parfor chnumber=1:64

optnet=trainDNN(chnumber).optnet;
optlrtheta=trainDNN(chnumber).optlrtheta;
para=trainDNN(chnumber).para;
ra=sigx1(chnumber,:); 
ra2=sigx2(chnumber,:); 
  
testx=[];
filestartp=0;
for testfni=1:testfn
    startpoint =1;
    for m=1:fifran(testfni)
    startpoint =filestartp+(m-1)*winShift; 

    r=ra(startpoint+1:startpoint+winLength);
    r2=ra2(startpoint+1:startpoint+winLength);
    energyr=r.^2;                   %calculaty energy of speech
    energyr2=r2.^2;                 
    %gfccf=get_gfccfeature(r,16);   %You can use the GFCC add this   
    tER1=zeros(1,unitsize); 
    tIID=zeros(1,unitsize);    
    tnITD=zeros(1,tdnumber);
    tITD=zeros(1,tdnumber); 
	%extract IID feature  
    for gi=1:unitsize
	        tER1(gi)=sum(energyr((gi-1)*egroupsize+1:gi*egroupsize));
        tIID(gi)=sum(energyr((gi-1)*egroupsize+1:gi*egroupsize))/(sum(energyr2((gi-1)*egroupsize+1:gi*egroupsize))+eps);
        tIID(gi)=10*log10(tIID(gi)+eps);        
     end;
     %extract ITD feature   
    for gi=1:tdnumber  
       sl1=r((gi-1)*groupsize+1:(gi-1)*groupsize+tdlength)-mean(r((gi-1)*groupsize+1:(gi-1)*groupsize+tdlength));
       sr1=r2(1:tdlength)-mean(r2(1:tdlength));
        energynr=(sl1*sl1')^0.5;
        energynr2=(sr1*sr1')^0.5;        
        tnITD(gi)=sl1*sr1'/(energynr*energynr2+eps);
        
        sl2=r(1:tdlength)-mean(r(1:tdlength));
        sr2=r2((gi-1)*groupsize+1:(gi-1)*groupsize+tdlength)-mean(r2((gi-1)*groupsize+1:(gi-1)*groupsize+tdlength));
        energynr=(sl2*sl2')^0.5;
        energynr2=(sr2*sr2')^0.5;        
        tITD(gi)=sl2*sr2'/(energynr*energynr2+eps);
        end
        
        tITDA=[fliplr(tnITD(1,2:size(tnITD,2))) tITD(1:size(tITD,2)-1)];% CCFR feature
        %you can add other monaural features and binaural features there
        %energynra=(r*r')^0.5;
        %energynra2=(r2*r2')^0.5;        
        %tIIDA=20*log10(energynr/energynr2);
    
        testf= [tIID tITDA];    
        testx=[testx;testf];  
    end
    filestartp=startpoint+winLength;%  the file end add one shift.
    
end

test_x=testx;
test_y=testy(:,chnumber);

pred=DNNTESTa(para,optnet,optlrtheta,test_x,test_y);

outputpreda=[outputpreda pred];
[acc1 hit1 fa1]=getHITFA(test_y(:),pred(:));
outacc=[outacc acc1];
outhit=[outhit hit1];
outfa=[outfa fa1];
testya=[testya test_y];
emaska=[emaska pred];
chnumber
end;
matlabpool close;
%% result analyze
[acc hit fa]=getHITFA(testya(:),outputpreda(:));
HIT=hit*100;
FA=fa*100;
ACC=acc*100;
HITFA=(hit-fa)*100;
end