function Hy = WienerFilter( u, y, Hn )
% Hy = WienerFilter( u, y, Hn )
% identify or simulate a MIMO system modeled by the discrete time convolution:
%
%     y(j) = sum_{k=1}^N H(k) u(j-k+1)
%
%
%     INPUTS            DESCRIPTION
%   ==========        ===============
%       u             sequence of input data                       (r x N)
%       y             sequence of outputs                          (m x N)
%      Hn             if Hn is a scalar it is interpeted as the order 
%                     of Markov parameters, n, and the fitted Markov
%                     parameters, H, of dimension (m x rn),  are returned.
%                     otherwise,  Hn is interpreted as a set of Markov 
%                     parameters, H, of dimension (m x rn), and the
%                     model output y, of dimension (m x N), is returned
%
%                     note, the Wiener filter coefficients are 
%                     returned in reverse order [ H(n), H(n-1), ..., H(1) ] 
%
%    OUTPUTS            DESCRIPTION
%   ==========        ===============
%      Hy             Hn is a scalar, Hy contains the Markov parameters (m x rn)
%                     otherwise, Hy contains the model output (m x N)
%
% Wiener, Norbert, The extrapolation, Interpolation, and Smoothing of 
% Stationary Time Series, John Wiley, New York, 1949.


 [r,Nu] = size(u);	% number of  inputs, duration of  inputs
 [m,Ny] = size(y);	% number of outputs, duration of outputs

 N = min(Nu,Ny);

 if max(size(Hn)) == 1	% fit the model to input/output data sequences
    n = Hn;
    if ( r*n > N-n+1 )
      error(' WienerFilter: not enough data to fit');
    end
 else			% simulate the model response with input data sequences
    H = Hn;
    n = size(H,2)/r;
 end

 U = zeros(r*n,N-n+1);
 for k=1:n            % assemble the Hankel matrix of the time sequence u
    U ( r*(k-1)+1:r*k , : ) = u ( : , k:N-n+k );
 end	

 if max(size(Hn)) == 1     % --- fit the model to data

   Ruu = (1/(N-n)) * U*U';	  %  auto-correlation of input data
   Ryu = (1/(N-n)) * y(:,n:N)*U'; % cross-correlation of input/output data

   Hy  = Ryu / Ruu;               % Wiener-Hopf equations 

 else			   % --- evaluate the model output 

   Hy = [ zeros(m,n-1) H*U ];

 end

 clear U

%{
 Su = svd(Ruu);
 figure(2)
  clf
  plot(Su/Su(1,1),'o')
   ylabel('\sigma_i/\sigma_1,   singular values of Ruu')
   xlabel('i')
%  axis([0 m*n+1 -1 1.1*Su(1,1)])
%}
% --------------------------------------- WienerFilter   HP Gavin 
% System Identification, Duke University, Fall 2013, 
% updated:  2013-09-04, 2017-10-06
