% Copyright (c) 2005, International Computer Science Institute
%
% All rights reserved.
% 
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions
% are met:
% 
% Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
% 
% Redistributions in binary form must reproduce the above copyright
% notice, this list of conditions and the following disclaimer in the
% documentation and/or other materials provided with the distribution.
% 
% Neither the name of ICSI nor the names of its contributors may be
% used to endorse or promote products derived from this software
% without specific prior written permission.
% 
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
% "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
% LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
% A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
% OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
% SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
% TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

function output=pbf(alData, fs, gamma, k)
%Multiple-microphone phase-error based processing.

%June 2005, Marc Ferras, International Computer Science Institute

gamma
k

%Determine the number of input channels
nChannels=length(alData);
%Length of the input waves
lData=length(alData{1});

%Frame length for STFT processing
%frameLengthSec=30e-3;
frameLength=1024;
fftLength2=floor(frameLength/2)+1;

%Frame shift for STFT processing
frameShiftSec=10e-3;
frameShift=floor(frameShiftSec*fs);

%Window function for STFT processing
window=hanning(frameLength)./sum(abs(hanning(frameLength)));

%Initialize output vector
output=zeros(lData,1);

%Initialize number of processed frames
nFrame=0;

iniFrameSample=1;
endFrameSample=iniFrameSample+frameLength-1;
while endFrameSample<lData
    %A new frame to process
	nFrame=nFrame+1;
    
    %Get short-time magnitude and phase spectrums for each input channel
	for ch=1:nChannels
		frame=alData{ch}(iniFrameSample:endFrameSample);
		wFrame=window.*frame;
		fdFrame{ch}=fft(wFrame);
		magFrame{ch}=abs(fdFrame{ch});
		phaseFrame{ch}=angle(fdFrame{ch});
    end

	
    %Initialize product of masks
	prodHij=ones(size(frame));
    
    %Initialize enhanced frame vector 
	enhFrame=zeros(size(frame));
    
    %Process all possible pairs of channels
	for ch1=1:nChannels
		for ch2=1:nChannels
			if ch1 ~= ch2
                
                %First way of computing phase error
                %phaseError=phaseFrame{ch1}-phaseFrame{ch2};
				%phaseError=acos(cos(phaseError));
                
                
				%Second way of computing phase error
				phaseError1way=abs(phaseFrame{ch1}-phaseFrame{ch2});
				
				negAngles1 = phaseFrame{ch1}<0;
				tmp1=phaseFrame{ch1};
				tmp1(negAngles1) = tmp1(negAngles1)+2*pi;
				negAngles2 = phaseFrame{ch2}<0;
				tmp2=phaseFrame{ch2};	
				tmp2(negAngles2) = tmp2(negAngles2)+2*pi;
				phaseError2way=abs(tmp1-tmp2);

				phaseError=min(phaseError1way,phaseError2way);

                %Compute mask for ch1-ch2 pair 
                h12=1./(1+gamma.*phaseError.*phaseError);
                
                %Update product of masks
				prodHij=prodHij.*h12;
			end
        end
        %Compute geometrically averaged mask
		phii=prodHij.^(1/k);	
        
        %Apply mask onto input frame
		enhFrame=enhFrame+phii.*fdFrame{ch1};
    end

    %Get time-domain frame
	tdFrame=real(ifft(enhFrame));
    %Window it to avoid time-aliasing
	wtdFrame=window.*tdFrame(1:frameLength);

    %Overlap-add to synthesize the output waveform
	output(iniFrameSample:endFrameSample)=output(iniFrameSample:endFrameSample)+wtdFrame;

    %Update frame boundaries
	iniFrameSample=iniFrameSample+frameShift;
	endFrameSample=endFrameSample+frameShift;
end

return;