function [dp, ep, dq, eq] = saoudi_diag(alphap, alphaq);
% Syntax: [dp, ep, dq, eq] = saoudi_diag(alphap, alphaq);
% 
% Calculates the coefficients necessary for the the main and secondary (lower)
% diagonal of Mp/2 and Mstarp/2 matrices.
%
% Function based on article Saoudi, Boucher, Guyader:
% "A new efficient algorithm to compute the LSP parameters for speech coding" 
% Signal Processing, no. 28, 1992, pp. 201-212.
% 
% alphap is vector of reccurence coefficients of the symetric polynomial,
%   size P+1, value of alpha(1) without importance (alpha0 does not exist)
% alphaq is vector of reccurence coefficients of the antisymetric polynomial,
%   size P+1, value of alpha(1) without importance (alpha0 does not exist)
% dp, dq are the elements of the main diagonal of matrice Mp/2, resp. Mstarp/2
%   length P/2
% ep, eq are the elements of the second. diagonal of matrice Mp/2, 
%   resp. Mstarp/2, length P/2, but the 1st element does not exist
%
% - all vectors are COLUMN
% - The function works only for P even
%
% The implementation in Matlab is not "vectorized", so it may be very slow

P=length(alphap)-1;
if(P/2 ~= floor(P/2) )
 error ('P is not even');
end

%%% init %%%
md=P/2;
dp=zeros(md,1); dq=zeros(md,1);		% both matrices have size md x md
ep=zeros(md-1,1); eq=zeros(md-1,1);

%%% 1) the sym. part - matrice Mp/2. 1/2 of the function
%%%    symetrique of the article

alphap(2) = 2 * alphap(2);
dp(1) = alphap(2) + alphap(3) - 2;

for i=1:md-1,
  j=i*2;
  dp(i+1)=alphap(j+2) + alphap(j+3) -2;
  ep(i+1)=alphap(j+1) * alphap(j+2);
end

%%% 2) the antisym. part - matrice Mstarp/2. 1/2 of the function
%%%    anti_symetrique of the article 

dq(1) = alphaq(3) - 2;

for i=1:md-1,
  j=i*2;
  dq(i+1)=alphaq(j+2) + alphaq(j+3) -2;
  eq(i+1)=alphaq(j+1) * alphaq(j+2);
end