1%% CCAF Bound Calculator
 2% This function calculates the constant upper and lower matrices phi
 3% and psi for CCAF's for the Blocking Matrix (BM) 
 4%
 5% *Syntax*
 6% |[phi, psi] = ccafbounds(m, fs, c, p, order)|
 7% *Inputs*
 8% *|mpos|-microphone positions in R^3 (meters)
 9% *|fs|-sample rate (Hertz)
10% *|soundc|-Speed of sound (meters/sec)
11% *|prop|-Estimated propagation time across the array in samples
12% *|order|-Order of the adaptive filters
13%
14% *Outputs*
15% * |phi| - Matrix of upper bounds, where each column is a vector of
16%      bounds for a the adaptive filter of a single track in the BM.
17% * |psi| - Matrix of lower bounds with the same structure as psi.
18%
19% *Notes*
20%
21% This code is based on equations derived for a linear microphone
22% array where the effect of the coefficient bounds
23% would be to expand the main lobe of the beamformer depending on the
24% sine of a parameter delta-theta, the supposed error in steering
25% angle of the beamformer.  However, here we consider a beamformer of
26% arbitary geometry and thus the parameter "delta-theta" no longer
27% makes sense.  Our current adaptation of the algorithm is
28% 
29% # Consider the center of the array as the centroid of array,
30%   calculated as the arithmetic mean of the mic coordinates.
31% # Hard-code a value for sin(delta-theta), which we know must be
32%   bounded [-1 1].  At present we've selected .05
33% 
34% Note that these adaptations will still recreate the original
35% performance of the algorithm for a linear array, where our selection
36% of .05 for sin(delta-theta) corresponds to a steering angle error of
37% about +/- 3 degrees.
38
39
40%% Function Declaration
41function [phi, psi] = ccafbounds(mpos, fs, soundc, prop, order)
42
43%% Argument Error Checking
44narginchk(5, 5);
45if ~isreal(mpos) || ~all(all(isfinite(mpos))) || size(mpos,1) ~= 3
46    error('mpos must be a real matrix with 3 rows')
47elseif ~isscalar(fs) || ~isreal(fs) || ~isfinite(fs) || ...
48        fs < 0 || abs(mod(fs, floor(fs))) > 0
49    error('fs must be a positive integer')
50elseif ~isreal(soundc) || ~isscalar(soundc) || ~isfinite(soundc) || soundc <= 0
51    error('soundc must be positive real scalar');
52elseif ~isscalar(prop) || ~isreal(prop) || ~isfinite(prop) || ...
53        prop < 0 || abs(mod(prop, floor(prop))) > 0
54    error('prop must be a positive integer')
55elseif ~isscalar(order) || ~isreal(order) || ~isfinite(order) || ...
56        order < 0 || abs(mod(order, floor(order))) > 0
57    error('order must be a positive integer')
58end
59
60%% Calculate Bounds
61% In the original paper, the bound vectors for each adaptive filter are
62% calculated as
63%
64% $$\phi_{m,n} = \frac{1}{\pi\ max(.1, \ 
65% (n\ ^\_ \ P)\ ^\_ \ T_m, \ ^\_ (n\ ^\_ \ P)\ ^\_ \ T_m)} \quad
66% \psi = \ ^\_ \phi\ \forall\ m, n$$
67%
68% $$T_m = \frac{b_mf_s}{c}\sin\Delta\theta $$
69%
70% where bm is the distance of the mth microphone to the center of the
71% array.  Remember that we must fudge for sin(delta-theta) in R^3.
72
73M = size(mpos,2);  % number of microphones in the array
74phi = zeros(order, M);  % initialze upper bound matrix for iteration
75
76for mIter = 1:M   % iterate over all microphones
77    sinDt = .05;  % kludge for 3-D (see Notes section above)
78    arrayCentroid = mean(mpos,2);  % use centroid as "center" of array
79    bm = sqrt(sum((mpos(:,mIter)-arrayCentroid).^2));  % Get mic distance
80                                                    % from centroid
81    Tm = bm*fs*sinDt/soundc;  % Hoshuyama equation
82    for nIter = 1:order  % Set bound for each tap of this adaptive filter
83        phi(nIter, mIter) = 1/(pi*max([.1, (nIter-prop)-Tm, -(nIter-prop)-Tm])); 
84         % directly from Hoshuyama paper
85    end
86end
87psi = -phi;  % psi is simply the opposite of phi
88
89end