function [Qk,Lk_star]=QLpoly(D,Beta,r,k) 
%function [Qk,Lk_star]=QLpoly(D,Beta,r,k) 
%
% This function determines the polynomials Qk and Lk_star, 
% associated with the design of WLMS algorithms, by using a
% closed-form solution to a related Diophantine equation.
%
% The design of the WLMS is based on ARIMA-modelling of
% time-varying FIR parameters.
% 
% 
% Inputs: 
%
%          D - denominator of the ARIMA model
%       Beta - stable spectral factor, see "spefac.m"
%          r - scalar associated with "spefac.m".
%          k - prediction horizon (k>0), fixed-lag smoothing (k<0)
%              filtering (k=0)
% 
% 
%       [Qk,Lk_star]=QLpoly(D,Beta,r,k) 
% 
%        
% L. Lindbom (2001-12)
 
nB=length(Beta)-1; 
nD=length(D)-1; 
nC=nB-nD;

if nC>0, D=[D zeros(1,nC)];end

%*********** initialize with k=0 (filtering) *************

Qk=(r*Beta(1:nB)-D(1:nB))/r;
Lk_star=conj(Beta(2:nB+1)-D(2:nB+1)); 
  

if k > 0                %Prediction
   
   for i=1:k 
       
      Qk_0=Qk(1); 
           tmpQ=[Qk 0]-D*Qk_0; 
      Qk=tmpQ(2:nB+1); 
       
      nLk=max(nC+i,nB)-1; 
      tmpL=[0 Lk_star]+[conj(r*Beta)*Qk_0 zeros(1,nLk-nB)]; 
      Lk_star=tmpL(1:nLk+1);
   end 
    
elseif k < 0            %Fixed-lag smoothing
    
   for i=1:abs(k) 
       
      Lk_0=Lk_star(1); 
      tmpL=[Lk_star 0]-conj(Beta)*Lk_0; 
      Lk_star=tmpL(2:nB+1); 
       
      nQk=max(nC+i,nD-1); 
      tmp=Lk_0/r; 
       
      tmpQ=[0 Qk]+[D*tmp zeros(1,nQk-nD)]; 
      Qk=tmpQ(1:nQk+1); 
       
   end 
          
end