Code covered by the BSD License

# Precedence-based cross-correlograms

### Christopher Hummersone (view profile)

27 May 2010 (Updated )

Calculate cross-correlograms with a wide range of options

ch_xcorr(hc_L,hc_R,fs,varargin)
```function [ccg,ic] = ch_xcorr(hc_L,hc_R,fs,varargin)
% Calculate cross-correlograms with a wide range of options.
%
%   ccg = ch_xcorr(hc_L,hc_R,fs)
%   ccg = ch_xcorr(hc_L,hc_R,fs,'parameter',value)
%   [ccg,ic] = ch_xcorr(...)
%
%   ccg = ch_xcorr(hc_L,hc_R,fs) cross-correlates the input
%   2-D matrices hc_L and hc_R over 10ms frame with a
%   maximum lag of 1ms. 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 each frame. The number of
%   frames frame_count is calculated thus:
%
%   frame_count = ...
%       floor((max(size(hc_L))-maxlag-1)/frame_length);
%
%   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(hc_L,hc_R,fs,'parameter',value) allows a
%   number of options to be specified. The options are:
%
%   ({} indicates the default value)
%
%   'frame_length'   : {round(0.01*fs)} | scalar
%       The length of frames used to calculate for
%       integrating cross-correlations.
%   'noverlap'       : {1} | scalar
%       The number of frames over which to integrate the
%       cross-correlations. Note that the frame count is
%       reduced accordingly.
%   'maxlag'         : {round(0.001*fs)} | scalar
%       The maximum lag of the cross-correlation.
%   'tau'            : {round(0.01*fs)} | scalar
%       The time constant of the exponential window used to
%       calculate running cross-correlations.
%   'inhib'          : {[]} | array
%       Specificies an array with which to multiply the
%       cross-correlations before they are integrated. The
%       value defaults to an empty array, meaning that no
%       inhibition will be applied.
%   'ic_t'           : {0} | scalar
%       Specifies the interaural coherence (IC) threshold.
%       Only samples for which the IC exceeds this threshold
%       will be used to integrate cross-correlations. The
%       algorithm calculates Interaural Coherence (IC)
%       according to [1]. The value should be in the range
%       [0,1];
%   'norm_flag'      : {0} | scalar
%       Specifies whether the cross-correlograms are
%       calculated using normalised cross-correlations. A
%       non-zero value indicates that normalised
%       cross-correlations are used.
%   'inhib_mode'     : {'subtract'} | 'multiply'
%       Specify how the inhibition is applied. The default
%       'subtract' will subtract inhib from the running
%       cross-correlations; 'multiply' will multiply inhib
%       with the running cross-correlations.
%
%   [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*frame_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: 2015-02-24 12:53:09 +0000 (Tue, 24 Feb 2015) \$
% Last committed:   \$Revision: 320 \$
% Last changed by:  \$Author: ch0022 \$
% ==========================================================
% !---

assert(nargin>=3,'Number of input arguments must be greater than or equal to three.')

if isparameter(varargin,'inhib_mode') && ~isparameter(varargin,'inhib')
warning('''inhib_mode'' specified, but no inhibition array ''inhib''.');
end

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

options = struct(...
'frame_length',round(0.01*fs),...
'noverlap',1,...
'maxlag',round(0.001*fs),...
'tau',round(0.01*fs),...
'inhib',[],...
'ic_t',0,...
'norm_flag',0,...
'inhib_mode','subtract');

if nargin > 3 % if parameters are specified
optionNames = fieldnames(options);
% count arguments
nArgs = length(varargin);
if round(nArgs/2)~=nArgs/2
error('CH_XCORR needs propertyName/propertyValue pairs')
end
% overwrite defults
for pair = reshape(varargin,2,[]) % pair is {propName;propValue}
IX = strcmpi(pair{1},optionNames); % find match parameter names
if any(IX)
% do the overwrite
options.(optionNames{IX}) = pair{2};
else
error('%s is not a recognized parameter name',pair{1})
end
end
end

% assign options to variables
frame_length = options.frame_length;
noverlap = options.noverlap;
maxlag = options.maxlag;
tau = options.tau;
inhib_mode = options.inhib_mode;
norm_flag = options.norm_flag;
ic_t = options.ic_t;
inhib = options.inhib;

% check inputs
assert(all(size(hc_L)==size(hc_R)),'''hc_L'' and ''hc_R'' must be the same size')
assert(round(frame_length)==frame_length && isscalar(frame_length) && frame_length>0,'''frame_length'' must be an integer greater than zero')
assert(round(noverlap)==noverlap && isscalar(noverlap) && noverlap>0,'''noverlap'' must be an integer greater than zero')
assert(round(maxlag)==maxlag && isscalar(maxlag) && maxlag>0,'''maxlag'' must be an integer greater than zero')
assert(isscalar(tau) && tau>=1,'''tau'' must be a scalar greater than or equal to one')
assert(isscalar(norm_flag),'''norm_flag'' must be a scalar')
assert(isscalar(ic_t) && ic_t>=0 && ic_t<=1,'''ic_t'' must be a scalar in the range [0,1]')
assert(ischar(inhib_mode),'''inhib_mode'' must be a char array (string)')

% 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 and inhib 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

% set inhibition mode ID
switch inhib_mode
case 'multiply'
inhib_mode_ID = 1;
if isempty(inhib)
inhib = ones(size(hc_L));
end
case 'subtract'
inhib_mode_ID = 2;
if isempty(inhib)
inhib = zeros(size(hc_L));
end
otherwise
error('''inhib_mode'' must be set to ''multiply'' or ''subtract''')
end

inhib = check_input(inhib,2,numchans);

% 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)];

assert(all(size(inhib)==size(hc_L)),'''inhib'' must be a matrix the same size as ''hc_L'' or ''hc_R''')

% 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

% ----------------------------------------------------------
% isparameter: check for input parameter
% ----------------------------------------------------------
function set = isparameter(input,parameter)

set = any(strcmpi(input(cellfun(@ischar,input)),parameter));

% [EOF]
```