1%% Multiple-Input Canceller LMS Filter
 2
 3%% Function Declaration
 4function [y, mcWall] = mclms(z, b, mu, order, beta, zPrev, mcWinit, ...
 5    mcWForce, mcAdapt, K)   
 6
 7%% Argument Error Checking
 8narginchk(10, 10);
 9if ~isreal(z) || length(size(z)) ~= 2 || ~all(all(isfinite(z)))
10    error('z must be a real matrix');
11elseif ~isvector(b) || ~isreal(b) || ~all(isfinite(b))
12    error('b must be a real vector');
13elseif isempty(mu) || ~isscalar(mu) || ~isreal(mu) || ~isfinite(mu)
14    error('mu must be a real scalar');
15elseif isempty(order) || ~isscalar(order) || ~isreal(order) || ...
16        ~isfinite(order) || order < 0 || ...
17    abs(mod(order,floor(order))) > 0
18    error('order must be a positive integer'); 
19elseif isempty(beta) || ~isscalar(beta) || ~isreal(beta) || ...
20        ~isfinite(beta) || beta < 0 || beta > 1
21    error('beta must be a real scalar in [0,1]');
22elseif length(size(zPrev)) ~= 2 || ~isreal(zPrev) || ...
23        ~all(all(isfinite(zPrev)))
24    error('zPrev must be a real matrix');
25elseif ~isempty(mcWinit) && ...
26    (length(size(mcWinit)) ~= 2 || ~isreal(mcWinit) || ...
27        ~all(all(isfinite(mcWinit))))
28    error('mcWinit must be a real matrix');
29elseif ~isempty(mcWForce) && ...
30    (~isreal(mcWForce) || length(size(mcWForce)) ~= 3 || ...
31     size(mcWForce,1) ~= order || ...
32     size(mcWForce,2) ~= size(z,2) || ...
33     ~all(all(all(isfinite(mcWForce)))))
34    error('mcWForce must be a real cubic matrix');
35elseif ~isempty(mcAdapt) && ...
36    (~isreal(mcAdapt) || ~all(all(isfinite(mcAdapt))))
37    error('mcAdapt must be a real matrix');
38elseif ~isempty(K) && ...
39    (~isscalar(K) || ~isreal(K) ||  ~isfinite(K) || K < 0)
40    error('K must be a real positive scalar');
41end
42
43%% Adaptive Filtering
44% Apply the LMS algorithm with possibly some tap controls.  The LMS
45% adaptation algorithm applied to each audio track is
46%
47% $${\bf w}_i(n+1) = \beta {\bf w}_i(n) +
48% \mu\frac{y(n)}{{\bf z}_i^T(n){\bf z}_i(n)}{\bf z}_i(n),
49% \quad \beta \in [0, 1] , \quad \mu \in [0, 2], \quad i \in [1, M] $$
50%
51% where the algorithm is "leaky" if beta < 1.  Note that, as in all
52% GSC algorithms, the output signal y is used directly as the error
53% signal in these filters, thus minimizing the overall output power of
54% the beamformer.  In addition, if K is specified, we have a
55% norm-constrained adaptive filter (NCAF), which restrains the total
56% power that can pass through the filter if the tap norm exceeds K, as
57% described by
58%
59% $$\Omega_i = {\bf w}_i(n+1){\bf w}_i^T(n+1)$$
60%
61% $${\bf w}_i(n+1) = \sqrt{\frac{K}{\Omega_i}}{\bf w}_i(n+1), 
62% \quad \Omega_i > K$$
63%
64% The output of the beamformer is calculated as the Fixed-Beamformer
65% Output minus the sum of the filter outputs.
66%
67% $$ y(n) = b(n) \ ^\_ \ \sum_{i=1}^{M} {\bf w}_i'(n){\bf z}_i(n) $$
68%
69% Each column of w holds the taps for an adaptive filter on a column
70% of the input microphone tracks x.  mcWall saves a 3-D matrix of all
71% values of w.  Use mcWinit for initial value of w if supplied; zeros
72% otherwise.
73%% inital values of MCLMS 
74
75if isempty(zPrev),  zPrev = zeros(order, size(z,2));  end
76mcWall = zeros(order, size(z,2), length(b));
77
78% initial targmc
79y = zeros(size(b));
80if ~isempty(mcWinit)
81w = mcWinit; else w = zeros(order, size(z,2)); 
82end
83
84
85% ord--one for each track of BM
86for n = 1:length(b)   
87% calculate targmc output and ready to opt
88[zWin,y,w1]=targmc(zPrev,z,b,n,order,w,y);
89% Tap update. If mcAdapt specified, adapt only requested columns
90[w]=optwmc(zWin,n,w1,mcAdapt,mu,order,beta,z,mcWForce);
91% NCAF (Norm-Constrained Adaptive Filter)
92[w]=NCAF(K,w,order,n,mcAdapt);
93% Save the taps used on this iteration
94mcWall(:,:,n) = w;
95end  % function mclms