1 function [W, R, P] = miso_firwiener(N, X, y)
   2 %MISO_FIRWIENER Optimal FIR Wiener filter for multiple inputs.
   3 %   MISO_FIRWIENER(N, X, Y) computes the optimal FIR Wiener filter of 
   4 %   order N, given any number of (stationary) random input signals as 
   5 %   the columns of matrix X, and one output signal in column vector Y.
   6 
   7 %   Author: Keenan Pepper
   8 %   Last modified: 2007-12-21
   9 
  10 %   References:
  11 %     [1] Y. Huang, J. Benesty, and J. Chen, Acoustic MIMO Signal
  12 %     Processing, Springer−Verlag, 2006, page 48
  13 
  14 % Number of input channels.
  15 M = size(X, 2);
  16 
  17 % Input covariance matrix, in abbreviated block Toeplitz form.
  18 R = zeros(M*(N+1), M);
  19 for m = 1:M
  20     for i = 1:M
  21          rmi = xcorr(X(:,m), X(:,i), N);
  22          Rmi = flipud(rmi(1:N+1));
  23          top = (m-1) * (N+1) + 1;
  24          bottom = m * (N+1);
  25          R(top:bottom,i) = Rmi;
  26     end
  27 end
  28 
  29 R = reshape(R, [N+1,M,M]);      % This just permutes the indices to
  30 R = permute(R, [2,1,3]);        % change R from Toeplitz-block to
  31 R = reshape(R, [M*(N+1),M]);    % block-Toeplitz.
  32 
  33 % Cross−correlation vector.
  34 P = zeros(1, M*(N+1));
  35 for i = 1:M
  36     top = (i-1)*(N+1)+1;
  37     bottom = i * (N+1);
  38     p = xcorr(y, X(:,i), N);
  39     P(top:bottom) = p(N+1:2*N+1)';
  40 end
  41 
  42 P = reshape(P, [N+1,M]);        % More index rearrangement.
  43 P = reshape(P', [M*(N+1),1]);
  44 
  45 W = block_levinson(P, R);
  46 
  47 W = reshape(W, [M,N+1]);        % Make output compatible with earlier
  48 W = reshape(W', [1, M*(N+1)]);  % version.