function [Qk,Q1,Q0,Beta]=wlms_design(gamma,k,C,D) 
%function [Qk,Q1,Q0,Beta]=wlms_design(gamma,k,C,D) 
%
% Design of filters used in the Wiener LMS (WLMS) adaptation law.
% WLMS tracks the parameter vector h(t) in a linear regression 
% model of a signal y(t), based on an ARIMA model of h(t):
%
%          C(1/z)
%   h(t)=  ------ I e(t) ;  y(t)= regressors*h(t) + noise
%          D(1/z)  
%
% This function calculates the following Wiener LMS filter polynomials:
%
%       Qk(1/z) - for arbitrary horizon k, with k being an integer
%                  (k>0,prediction; k<0, fixed-lag smoothing)
%       Q1(1/z) - for horizon k=1 (one-step ahead prediction)
%
%       Q0(1/z) - for horizon k=0 (filtering) 
%                (Note: mu = leading coefficient of Q0.)
%
%     Beta(1/z) - spectral factor, denominator of the Learning filter.   
%
% Inputs:
%
% 1. gamma - parameter drift-to-noise ratio 
%        (can be used as a tuning knob)
% 2.     k - horizon parameter
% 3.     C - numerator polynomial of the ARIMA-model
% 4.     D - denominator of the ARIMA-model
%
% Uses the functions spefac.m and QLpoly.m
%
% Reference: IEEE Trans. on Communications, vol 49, Dec. 2001. 
% L. Lindbom (2001-12-20), M. Sternad (2002-01-05)


%-------- solve spectral factorization -----------------
%
%    r*Beta(1/z)Beta(z)=gamma*C(1/z)C(z)+D(1/z)D(z)
%
%-------------------------------------------------------

nC=length(C)-1;
nD=length(D)-1;
if nC>nD
   D=[D zeros(1,nC-nD)];
elseif nC<nD
   C=[C zeros(1,nD-nC)];
end

rhs=gamma*xcorr(C)+xcorr(D);
[Beta,r]=spefac(rhs);

%----------- determine the WLMS filters ---------------

Qk=QLpoly(D,Beta,r,k);
Q1=QLpoly(D,Beta,r,1);
Q0=QLpoly(D,Beta,r,0);