Code covered by the BSD License

# Precedence-based cross-correlograms

27 May 2010 (Updated )

Calculate cross-correlograms with a wide range of options

ch_xcorr(hc_L,hc_R,frame_length,noverlap,maxlag,tau,varargin)
```function [ccg,ic] = ch_xcorr(hc_L,hc_R,frame_length,noverlap,maxlag,tau,varargin)
% Calculate cross-correlograms with a wide range of options.
%
%   ccg = ch_xcorr(hc_L,hc_R,frame_length,noverlap,...
%       maxlag,tau)
%   ccg = ch_xcorr(...,inhib)
%   ccg = ch_xcorr(...,ic_t,norm_flag)
%   ccg = ch_xcorr(...,inhib_mode)
%   [ccg,ic] = ch_xcorr(...)
%
%   ccg = ch_xcorr(hc_L,hc_R,frame_length,noverlap,maxlag,
%   tau) calculates cross-correlograms with a range of
%   cross-correlation parameters.
%
%   The function cross-correlates the input 2-D matrices
%   hc_L and hc_R over frames of length frame_length
%   samples. It is assumed that the number of frequency
%   channels is min(size(hc_L)) and hence hc_L and hc_R can
%   be in either orientation. The cross-correlograms consist
%   of cross-correlations for every frame and frequency
%   channel. ccg has dimensions [lag,frequency,frame]. The
%   function calculates running cross-correlations for every
%   sample and integrates these cross-correlations over
%   noverlap frames. The lags are determined by maxlag and
%   the cross-correlation is computed over the range of lags
%   -maxlag:maxlag, i.e., 2*maxlag+1 lags. The number of
%   frames frame_count is calculated thus:
%
%   frame_count =...
%       floor(...
%       (max(size(hc_L))-maxlag-1)/...
%       frame_length...
%        )-noverlap+1;
%
%   The underlying cross-correlation algorithm is based on
%   that proposed by Faller & Merimaa [1]. In this
%   implmentation, the time constant of the backward
%   infinite exponential window is given by tau (in
%   samples).
%
%   ccg = ch_xcorr(...,inhib) multiplies the running
%   cross-correlations with inhib before the
%   cross-correlations are averaged. If inhib is
%   un-specified it is set to ones(size(hc_L)).
%
%   ccg = ch_xcorr(...,ic_t,norm_flag) allows an IC
%   threshold to be specified and normalisation to be turned
%   on and off. The algorithm calculates Interaural
%   Coherence (IC) according to [1]. The cross-correlogram
%   is calculated for each frame by averaging only those
%   cross-correlations within the frame for which the
%   corresponding IC exceeds the IC threshold ic_t
%   (0<=ic_t<=1). Setting this value to 0 will result in all
%   cross-correlations contributing to the average
%   (default). This implementation also permits the
%   normalisation to be turned off by specifying norm_flag;
%   a non-zero value means normalisation will be applied
%   (default). This means, for example, that an IC theshold
%   could be employed, but the un-normalised
%   cross-correlations would contribute to the
%   cross-correlograms. Note that either both or neither of
%   these values must be specified.
%
%   ccg = ch_xcorr(...,inhib_mode) uses the mode inhib_mode
%   to apply the inhibition. The options are 'multiply'
%   (default), whereby the inhibitory signal inhib will be
%   multiplied with the running cross-correlation, or
%   'subtract', whereby inhib will be subtracted from the
%   running cross-correlation.
%
%   [ccg,ic] = ch_xcorr(...) returns the calculated IC to
%   the matrix IC. Although the matrix returned is the same
%   size as hc_L, IC is only calculated for samples
%   1:frame_count*fraeme_length, other values will be set to
%   0.
%
%   Algorithm
%
%   See the enclosed documentation for more details on the
%   workings of the algorithm and an important caveat.
%
%   References
%
%   [1] C. Faller and J. Merimaa, "Source localization in
%   complex listening situations: Selection of binaural cues
%   based on interaural coherence", The Journal of the
%   Acoustical Society of America, vol. 116, pp.3075-3089,
%   Nov. 2004.
%
%
%   C. Hummersone, R. Mason, and T. Brookes, "A comparison
%   of computational precedence models for source separation
%   in reverberant environments", in 128th Audio Engineering
%   Society Convention, London, May 2010, paper 7981.

% !---
% ==========================================================
% Last changed:     \$Date: 2013-03-08 11:04:10 +0000 (Fri, 08 Mar 2013) \$
% Last committed:   \$Revision: 235 \$
% Last changed by:  \$Author: ch0022 \$
% ==========================================================
% !---

assert(all(size(hc_L)==size(hc_R)),'''hc_L'' and ''hc_R'' must be the same size')

% Calculate frame count
frame_count = floor(max(size(hc_L))/(frame_length));
frame_count = frame_count-noverlap+1;

% Calculate number of frequency channels
numchans = min(size(hc_L));
numsamples = max(size(hc_L));

% Check orientation of HC data (i.e. that frequency runs across the rows)
dims = size(hc_L);
hc_L = check_input(hc_L,2,numchans);
hc_R = check_input(hc_R,2,numchans);
% set a flag if data has been transposed in this way
if dims(1)~=size(hc_L,1)
rot = true;
else
rot = false;
end

% Find inhibition data, if any (IS2D and FINDINPUTS defined below)
inhib = find_inputs(@(x)(isnumeric(x) & all(numel(x)>size(x)) & length(size(x))<3),...
varargin,'''inhib'' must be a matrix (same size as ''hc_L''), ''ic_t'' and ''norm_flag'' must both be specified and be scalars');
if isempty(inhib)
% If no data then the matrix is set to all ones
inhib = ones(size(hc_L));
else
% Check matrix data is in correct orientation
inhib = check_input(inhib,2,numchans);
assert(all(size(inhib)==size(hc_L)),'''inhib'' must be a matrix the same size as ''hc_L'' or ''hc_R''')
end

% Append HC and inhibition data with zeros for cross-correlation
hc_L = [hc_L; zeros(maxlag+1,numchans)];
hc_R = [hc_R; zeros(maxlag+1,numchans)];
inhib = [inhib; zeros(maxlag+1,numchans)];

% Find IC threshold and normalisation flag
[ic_t,norm_flag] = find_inputs(@isscalar,varargin,'''inhib'' must be a matrix (same size as ''hc_L''), ''ic_t'' and ''norm_flag'' must both be specified and be scalars');
if isempty(ic_t)
ic_t = 0;
norm_flag = 1;
end

% Find inhibition mode
inhib_mode = find_inputs(@ischar,varargin,'''inhib_mode'' must be a character array');
if isempty(inhib_mode)
inhib_mode = 'multiply';
end
switch inhib_mode
case 'multiply'
inhib_mode_ID = 1;
case 'subtract'
inhib_mode_ID = 2;
otherwise
error('inhib_mode must be set to ''multiply'' or ''subtract''')
end

% tau must be greater than or equal to 1
tau = max(tau,1);

% Check source file is compiled
check_mex_compiled('-largeArrayDims','ch_xcorr_c.c')

% Calculate cross-correlograms
[ccg,ic] = ch_xcorr_c(hc_L,hc_R,frame_count,frame_length,noverlap,maxlag,tau,inhib,ic_t,norm_flag,inhib_mode_ID);

% Correct orientation of IC data, if data was transposed, and crop to remove appended zeros
ic = ic(1:numsamples,:);
if rot
ic = ic';
end

% end of ch_xcorr()

% ----------------------------------------------------------
% Local functions:
% ----------------------------------------------------------

% ----------------------------------------------------------
% check_input: check input is correct orientation
% ----------------------------------------------------------
function output = check_input(input,dim,target)

if size(input,dim)~=target
output = input';
assert(size(output,dim)==target,'Input invalid')
else
output = input;
end

% ----------------------------------------------------------
% find_inputs: validate and provide inputs (if any)
% ----------------------------------------------------------
function varargout = find_inputs(fhandle,input,msg)

indices = cellfun(fhandle,input);

if any(indices) % inputs are specified
if sum(indices)~=nargout
error(msg) % return error message
end
varargout(1:nargout) = input(indices);
else % unspecified returns empty matrix
varargout(1:nargout) = {[]};
end

% [EOF]
```