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 order
   4 %   N, given any number of (stationary) random input signals as the columns
   5 %   of matrix X, and one output signal in column vector Y.
   6 
   7 %   Author: Keenan Pepper
   8 %   Last modified: 2007/08/02
   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.
  18 R = zeros(M*(N+1),M*(N+1));
  19 for m = 1:M
  20     for i = m:M
  21         rmi = xcorr(X(:,m),X(:,i),N);
  22         Rmi = toeplitz(flipud(rmi(1:N+1)),rmi(N+1:2*N+1));
  23         top = (m-1)*(N+1)+1;
  24         bottom = m*(N+1);
  25         left = (i-1)*(N+1)+1;
  26         right = i*(N+1);
  27         R(top:bottom,left:right) = Rmi;
  28         if i ~= m
  29             R(left:right,top:bottom) = Rmi';  % Take advantage of hermiticity.
  30         end
  31     end
  32 end
  33 
  34 % Cross-correlation vector.
  35 P = zeros(1,M*(N+1));
  36 for i = 1:M
  37     top = (i-1)*(N+1)+1;
  38     bottom = i*(N+1);
  39     p = xcorr(y,X(:,i),N);
  40     P(top:bottom) = p(N+1:2*N+1)';
  41 end
  42 
  43 % The following step is very inefficient because it fails to exploit the
  44 % block Toeplitz structure of R. It's done the same way in the built-in
  45 % function "firwiener".
  46 W = P/R;