
function data= bmdnntr(TfPath,TmPath,parameters)
%% 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)
%       .maxitd         step size in ms (default = 1)
%       .maxiid         IID dimensions (default = 2)
%       
% Output:
%   data                The trained DNN paprameters                         
%       
% Written by Yi Jiang (2013)

addpath(genpath(pwd));
fLIST=dir(fullfile(TfPath,'*.wav'));
mLIST=dir(fullfile(TmPath,'*.mat'));
%%Runing parameter
if nargin < 2
    parameters = struct;
end
% Parameter configuration
if ~isfield(parameters,'fn')
    parameters.fn = size(fLIST,1);      % sample rate
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 train environments
Traincorpus=TfPath;
channelnumber=parameters.nc;
fs=parameters.fs;                       % wav file sampling rate
fRange=parameters.fRange;               % default frequency range in Hz
FranmeT=parameters.winL;                %frame is 20ms
winLength = FranmeT*fs;                 % default window length in samp
winShift =parameters.step*fs;       

% frame shift (default is half frame)floor_value = 1e-5;

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=[];
optlrtheta=[];
trainay=[];
sigx1=[];
sigx2=[];

matlabpool local;
%%%%%%%% prepare train materials
trainfn=parameters.fn ;         % Train files amount
trainty=[];  
fsM=[];
trainpara=[];
gfccf=[];
k=0;                            %the train number 
framen=0;
for gtn=1:trainfn

   [tgsig,fs,nbits]=wavread([TfPath fLIST(gtn).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;
    fsM=[fsM M];                                                % record file frame number
  
    masklab=load([TmPath  mLIST(gtn,1).name]);
    mclab=masklab.mask';
    trainty=[trainty;mclab];
    gtn;
end;

%% train by channels
parfor chnumber=1:parameters.nc;
  traintx=[];  
  ra=sigx1(chnumber,:); 
  ra2=sigx2(chnumber,:);    
  %deal with each train file
  filestartp=0
  for tf =1:trainfn 
       startpoint =1;
       for m=1:fsM(tf)
            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;             %calculaty energy of speech
            %%You can add the T-F  feature  there
            %gfccf=get_gfccfeature(r,16);%GFCC FEATURE   
     
            %% extreact IID feature
            IID=zeros(1,unitsize);
            ER1=zeros(1,unitsize);
            ER2=zeros(1,unitsize);
            for gi=1:unitsize;         
                 ER1(gi)=sum(energyr((gi-1)*egroupsize+1:gi*egroupsize));
                 ER2(gi)=sum(energyr2((gi-1)*egroupsize+1:gi*egroupsize));
                 IID(gi)=sum(energyr((gi-1)*egroupsize+1:gi*egroupsize))/(sum(energyr2((gi-1)*egroupsize+1:gi*egroupsize))+eps);
                 IID(gi)=10*log10(IID(gi)+eps);         
            end;   
            %% extrance ITD featuere              
            nITD=zeros(1,tdnumber);
            ITD=zeros(1,tdnumber);
            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;        
                nITD(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;        
                ITD(gi)=sl2*sr2'/(energynr*energynr2+eps);
            end
            ITDA=[fliplr(nITD(1,2:size(nITD,2))) ITD(1:size(nITD,2)-1)];% CCFR feature
            %% get the features
            traintf= [IID ITDA];       
            %traintf= [gfccf.feat(:,m)' IID ITDA];
            traintx=[traintx;traintf]; 
       end
       filestartp=startpoint+winLength;%  the file end add one shift.
  end
 
%set up DNN train set
train_x=traintx;                %training data sets
train_y=trainty(:,chnumber);    %Training goal

%% DNN train on each channel
[pred acc optnet optlrtheta lastact para outputhit outputfa outputacc  outputpred]=DNNTraina(train_x,train_y,train_x,train_y);%train and test are same

%% result output
outputpreda=[outputpreda outputpred];       %estimated binary mask
trainay=[trainay train_y];
trainDNN(chnumber).optlrtheta=optlrtheta;
trainDNN(chnumber).optnet=optnet;
trainDNN(chnumber).para=para;
chnumber
end;
matlabpool close;
%% output  train result, if need.
[acc hit fa]=getHITFA(trainay(:),outputpreda(:))
data=trainDNN;
end
