1%% Griffiths-Jim Beamformer
  2
  3% *Required Inputs*
  4% * x-Matrix of raw audio data, each column a mic array track
  5% * fs-sampling rate (Hertz)
  6% * spos-3x1 source location point (meters)
  7% * mpos-column matrix of 3D mic positions (meters)
  8% * c- speed of sound (meters per second)
  9% * p-Constant delay placed on all tracks entering the BM
 10% * q-Constant delay placed on DSB output just before final summer
 11% * xPrev-Matrix of audio data from a previous window.
 12% * iw-(optional) if set to 1, the dsb will apply an inverse distance 
 13%          weighting on the microphones, useful for near-field/immersive
 14%          applications.  Any other value or its absense will result in
 15%          the uniform weight dbs.
 16
 17% *INPUT TO BMLMS*
 18% * mu-step size of the LMS filters in [0, 2]
 19% * order-order of the LMS filters
 20% * beta-forgetting factor of the LMS filters in [0, 1]
 21% * phi and psi-Upper and lower bounds for a CCAF blocking matrix 
 22% * bPrev-Vector of DSB output from a previous audio window to the BM LMS filters 
 23% * bmWinit-initial tap weights for the BM LMS filters. Set to [] if not needed  
 24%
 25% * bmWForce and mcWForce-Cubic matrices of tap weights to force
 26%     the BM and MC LMS filters to use.  Useful for recreating a
 27%     previous GJBF run where taps have already been determined.
 28%
 29% * |snrThresh| and |snrRate| - Parameters for selective adaptation.
 30%     |snrRate| times the SNR for each track is estimated as the
 31%     10*log10 ratio of the power of the fixed beamformer to each of
 32%     the blocking matrix output tracks.  If these estimated SNR's
 33%     exceed |snrThresh| then only the BM LMS coefficients are allowed
 34%     to adapt while the taps of the MC are kept constant; otherwise
 35%     the BM taps are locked and the MC coefficients adapt.  Must
 36%     specify both |snrThresh| and |snrRate| or neither.
 37% * |snrInit| - Assume this vector of SNR values at the start
 38
 39
 40% OPTION FOR MCLMS
 41% * |K| - Apply a norm constraint to the multiple-input canceller (MC)
 42%    LMS filters.  Leave blank for no constaint.
 43% * |zPrev| - Matrix of BM output for a previous audio window.  Used
 44%      as initial input to the MC LMS filters.  Supply zeros if no
 45%      previous data is available.
 46% * |mcWinit| - Matrix of initial tap weights for the MC LMS
 47%      filters.  Set to [] if not needed.
 48% OUTPUT
 49% * |b| - Vector of DSB output
 50
 51% *BMLMS Outputs*
 52% * |z| - Matrix of output tracks from the blocking matrix.approximate the noise 
 53%         in the simulation as closely as possible.  This can be checked 
 54%         to see how much target signal leakage occurs in the beamformer.
 55%      
 56% * |bmWall| - Cubic matrix of taps saved from the blocking matrix.
 57%      Each matrix (1st and 2nd dims) of the cube will be |order| tall
 58%      and |M| wide such that each column is the taps for an adaptive
 59%      filter for an individual channel.  The 3rd dimension specifies
 60%      the iteration of the taps.  This matrix can be fed into another
 61%      run as |bmWForce|.
 62%    mcAdapt - Binary matrix instructing the corresponding filters in
 63%     the MC to adapt.
 64% * |snrAll| - Matrix of estimated SNR's at every iteration in dB
 65%      (matrix of size NxM)
 66
 67
 68% *MCLMS Outputs*
 69% * |y| - Beamformed output audio track.
 70% * |mcWall| - Cubic matrix of taps saved from the multiple-input
 71%      canceller, referenced in the same mammer as |bmWall| and
 72%      applicable in subsequent runs as |mcWForce|.
 73
 74%% Function Declaration
 75function [y, bmw, mcWall, snrAll, b, z,mcAdapt] = ...
 76    gjbf(x, fs, spos, mpos, c, p, q, mu, order, beta, phi, psi, K, ...
 77     xPrev, bp, zp, bw1, mw1, snrThresh, ...
 78     snrRate, snr1, bmF, mcWForce,iw)
 79
 80%% Argument Error Checking
 81if nargin == 23
 82    iw = 0;
 83end
 84if isempty(x) || ~isreal(x) || ~all(all(isfinite(x)))
 85    error('x must be a real matrix')
 86elseif isempty(fs) || ~isreal(fs) || ~isscalar(fs) || ~isfinite(fs) ...
 87    || fs <= 0 
 88    error('fs must be positive real scalar');
 89elseif ~isreal(spos) || ~all(isfinite(spos)) || size(spos,1) ~= 3 || ...
 90        size(spos,2) ~= 1  || length(size(spos)) ~= 2
 91    error('spos must be a real 3x1 vector');
 92elseif ~isreal(mpos) || ~all(all(isfinite(mpos))) || size(mpos,1) ~= 3
 93    error('mpos must be real with 3 rows')
 94elseif ~isreal(c) || ~isscalar(c) || ~isfinite(c) || c <= 0
 95    error('c must be positive real scalar');
 96elseif isempty(p) || ~isscalar(p) || ~isreal(p) || ...
 97        ~isfinite(p) || p < 0 || abs(mod(p,floor(p))) > 0
 98    error('p must be a positive integer');
 99elseif isempty(q) || ~isscalar(q) || ~isreal(q) || ...
100        ~isfinite(q) || q < 0 || abs(mod(q,floor(q))) > 0
101    error('q must be a positive integer');
102elseif isempty(mu) || ~isscalar(mu) || ~isreal(mu) || ...
103        ~isfinite(mu)
104    error('mu must be a real scalar');
105elseif isempty(order) || ~isscalar(order) || ~isreal(order) || ...
106        ~isfinite(order) || order < 0 || ...
107        abs(mod(order,floor(order))) > 0
108    error('order must be a positive integer');
109elseif isempty(beta) || ~isscalar(beta) || ~isreal(beta) || ...
110        ~isfinite(beta) || beta < 0 || beta > 1
111    error('beta must be a real scalar in [0,1]');
112elseif ~isempty(phi) && (length(size(phi))~=2 || ...
113             ~isreal(phi) || ~all(all(isfinite(phi))))
114    error('phi must be a real matrix')
115elseif ~isempty(psi) && (length(size(psi))~=2 || ...
116             ~isreal(psi) || ~all(all(isfinite(psi))))
117    error('psi must be a real matrix')
118elseif xor(isempty(phi), isempty(psi))
119    error('phi and psi must both be specified or both empty')
120elseif ~isempty(K) && ...
121    (~isscalar(K) || ~isreal(K) ||  ~isfinite(K) || K < 0)
122    error('K must be a real positive scalar');
123elseif isempty(xPrev) || ~isreal(xPrev) || ~all(all(isfinite(xPrev)))
124    error('xPrev must be a real matrix')
125elseif ~isvector(bp) || ~isreal(bp) || ~all(isfinite(bp))
126    error('bp must be a real vector');
127elseif isempty(zp) || ~isreal(zp) || ~all(all(isfinite(zp)))
128    error('zp must be a real matrix')
129elseif ~isempty(bw1) && ...
130    (~isreal(bw1) || ~all(all(isfinite(bw1))))
131    error('bw1 must be a real matrix')
132elseif ~isempty(mw1) && ...
133    (~isreal(mw1) || ~all(all(isfinite(mw1))))
134    error('mw1 must be a real matrix')
135elseif ~isempty(snrThresh) && ...
136    (~isreal(snrThresh) || ...
137     ~isscalar(snrThresh) || ~isfinite(snrThresh))
138    error('snrThresh must be a positive real scalar')
139elseif ~isempty(snrRate) && ...
140    (~isreal(snrRate) || ...
141     ~isscalar(snrRate) || ~isfinite(snrRate) || snrRate <= 0)
142    error('snrRate must be a positive integer')
143elseif xor(isempty(snrThresh), isempty(snrRate))
144    error(['snrThresh and snrRate must both be specified or' ...
145       ' both empty'])
146end
147
148%% Setup
149[N, M] = size(x);  % number of samples in track, mics in array
150% SNR parameters
151if isempty(snrRate)  % no SNR thresholding, so assign dummy vals
152    snrAll = [];
153else  % apply SNR thresholds
154    snrAll = zeros(N, M);  % all estimated SNRs mat
155end
156%% Fixed Beamformer (FBF)
157% Find the distances of each mic to the source, calculate the
158% time-difference of arrival, and then add shifted waveforms. 
159[xShift,wghts] = delaytracking(x, xPrev, fs, spos, mpos, c, iw);
160b = sum(xShift, 2);  % add data across columns
161
162%% Blocking Matrix (BM)
163if isempty(phi)&& isempty(psi)  % traditional GJBF BM
164    z = diff(xShift, [], 2); % column differences
165    mcAdapt = []; % mc always adapts
166    bmw=[];  % Matlab complains if an output arg isn't set..
167else  % call function for adaptive BM
168    
169    % Apply p delays to each xShift column  
170    for k=1:size(x,2)  
171       xShift(:,k) = wghts(k)*[xPrev(end-p+1:end, k); ...
172               x(1:end-p, k)];
173    end
174    % Call external function to handle BM calculations
175    [z,bmw, mcAdapt, snrAll] =  ...
176    bmlms(xShift, b, mu, order, beta, phi, psi, bp, ...
177          bw1, bmF, snrThresh, snrRate, snr1);
178end
179
180
181%% Multiple Input Canceller (MC)
182% Apply q delay to the DSB vector before running the MC
183b = [bp(end-q+1:end); b(1:end-q)];
184% Call external function for calculations
185[y, mcWall] = mclms(z, b, mu, order, beta, zp, mw1, ...
186                    mcWForce, mcAdapt, K);
187
188end  % function gjbf