function [W,R,P] = miso_firwiener(N,X,y)
%MISO_FIRWIENER Optimal FIR Wiener filter for multiple inputs.
%   MISO_FIRWIENER(N,X,Y) computes the optimal FIR Wiener filter of order
%   N, given any number of (stationary) random input signals as the columns
%   of matrix X, and one output signal in column vector Y.

%   Author: Keenan Pepper
%   Last modified: 2007/08/02

%   References:
%     [1] Y. Huang, J. Benesty, and J. Chen, Acoustic MIMO Signal
%     Processing, Springer-Verlag, 2006, page 48

% Number of input channels.
M = size(X,2);

% Input covariance matrix.
R = zeros(M*(N+1),M*(N+1));
for m = 1:M
    for i = m:M
        rmi = xcorr(X(:,m),X(:,i),N);
        Rmi = toeplitz(flipud(rmi(1:N+1)),rmi(N+1:2*N+1));
        top = (m-1)*(N+1)+1;
        bottom = m*(N+1);
        left = (i-1)*(N+1)+1;
        right = i*(N+1);
        R(top:bottom,left:right) = Rmi;
        if i ~= m
            R(left:right,top:bottom) = Rmi';  % Take advantage of hermiticity.
        end
    end
end

% Cross-correlation vector.
P = zeros(1,M*(N+1));
for i = 1:M
    top = (i-1)*(N+1)+1;
    bottom = i*(N+1);
    p = xcorr(y,X(:,i),N);
    P(top:bottom) = p(N+1:2*N+1)';
end

% The following step is very inefficient because it fails to exploit the
% block Toeplitz structure of R. It's done the same way in the built-in
% function "firwiener".
W = P/R;