1 function x = block_levinson(y, L)
   2 %BLOCK_LEVINSON Block Levinson recursion for efficiently solving 
   3 %symmetric block Toeplitz matrix equations.
   4 %   BLOCK_LEVINSON(Y, L) solves the matrix equation T * x = y, where T 
   5 %   is a symmetric matrix with block Toeplitz structure, and returns the 
   6 %   solution vector x. The matrix T is never stored in full (because it 
   7 %   is large and mostly redundant), so the input parameter L is actually 
   8 %   the leftmost "block column" of T (the leftmost d columns where d is 
   9 %   the block dimension).
  10 
  11 %   Author: Keenan Pepper
  12 %   Last modified: 2007-12-21
  13 
  14 %   References:
  15 %     [1] Akaike, Hirotugu (1973). "Block Toeplitz Matrix Inversion".
  16 %     SIAM J. Appl. Math. 24 (2): 234-241
  17 
  18 s = size(L);
  19 d = s(2);                 % Block dimension
  20 N = s(1) / d;             % Number of blocks
  21 
  22 B = reshape(L, [d,N,d]);  % This is just to get the bottom block row B
  23 B = permute(B, [1,3,2]);  % from the left block column L
  24 B = flipdim(B, 3);
  25 B = reshape(B, [d,N*d]);
  26 
  27 f = L(1:d,:)^-1;          % "Forward" vector
  28 b = f;                    % "Backward" vector
  29 x = f * y(1:d);           % Solution vector
  30 
  31 for n = 2:N
  32     ef = B(:,(N-n)*d+1:N*d) * [f;zeros(d)];
  33     eb = L(1:n*d,:)' * [zeros(d);b];
  34     ex = B(:,(N-n)*d+1:N*d) * [x;zeros(d,1)];
  35     A = [eye(d),eb;ef,eye(d)]^-1;
  36     fn = [[f;zeros(d)],[zeros(d);b]] * A(:,1:d);
  37     bn = [[f;zeros(d)],[zeros(d);b]] * A(:,d+1:end);
  38     f = fn;
  39     b = bn;
  40     x = [x;zeros(d,1)] + b * (y((n-1)*d+1:n*d) - ex);
  41 end