function y = comp_lp(r,M,L);

% 	function y = comp_lp(r,M,L);
%
%		r : [r(0) ..  r(N)]  autocorrelation sequence
%		M : number of channels
%		L : number of freq. samples
%		y : output product filter, [y(0) y(1) ... y(N)]' 
%		    whose FT is nonnegative at L points	
%
%   	modifications:
%  8/4/96 	   the number of frequencies L has been 
%		   modified to be the same as in other programs
%                  tolerance to nonnegativity, delta has been removed			   

%  8/9/96	   only L/2 frequencies are considered because of the symmetry

%  8/12/96	   output is made one-sided: [y(0) y(1) ... y(N)]'


%    v is the nonzero components of g: [g1..gM-1 gM+1..g(2M-1) .. g(kM-1)] 
%    aim:   maximize   r'*v
%	    subject to [cosw cos2w .. com(M-1)w cos(M+1)w ..  cos(kM-1)w]*v >= -.5

%    equivalent problem:
%                     minimize r'*v
%		subject to [cosw cos2w .. com(M-1)w cos(M+1)w ..  cos(kM-1)w]*v <= .5


% if frequencies are evaluated at L points, then
%          constraint becomes W*v <= .5[1 1  ... 1]'
% where
%          W is Lxk(M-1) matrix for frequencies w=0, 2pi/L, ..., 2pi(L-1)/L
             

N = length(r)-1;		  % order of h

W1 = cos(2*pi*[0:L-1]'*[0:N]/L);    % Fourier matrix

W = W1;
W(:,1:M:N+1) = '';		  	% annihilate every Mth column

b = [.5*ones(L,1)];	  	

r(1:M:N+1)='';				%autocorrelation seq. to be used in opt.

v = lp(r,W,b);

gr = zeros(N,1);                    % right hand side of g
k=1:N;
k(M:M:N)='';
gr(k)=-v;			  % inserts zeros to g (actually inserts v to zeros :))	
y = [1;gr];



