function h=rNyquistMG(N,M,alpha,gmaZ,gmaT,eta);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Square-root Nyquist (M) filter design
%
% h=srNyquistM(N,M,alpha,gmaZ,gmaT,eta);
%
% parameters:
%   N: filter order (filter length = N+1)
%   M: number of samples per symbol period 
%   alpha: roll-off factor (range 0 to 1)
%   gmaZ: Weight factor for middle tap and zero crossings
%   gmaT: Weight factor for the tails of g=h*h (used for designs with
%   robust behavor against timing jitter)
%   eta:  Weight factor for tails of h (used for designs with reduced PAR)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Setup the weight matrix Gamma
Gamma=zeros(1,1+N); Gamma(M+2:end)=gmaT; Gamma(1:M:end)=gmaZ;  
Gamma=[Gamma ones(1,1+N/2)];Gamma=diag(Gamma);
%
%   Initial filter
h=sr_cos_p(N,M,alpha);
if rem(N+1,2)==0  h1=h(1:(N+1)/2); else h1=h(1:N/2+1); end
Lh1=length(h1);
%
% Set constraint matrices, Sn.
S=zeros(N+1,N+1,N+1);
temp=ones(N+1,1);
for n=1:N+1
    S(:,:,n)=spdiags(temp,-(n-1),N+1,N+1);
end
%
% Set the matrix Phi
Phi=zeros(N+1,N+1);
f2=(1/2/M)*(1+alpha);
Phi=[1-2*f2 -2*f2*sinc(2*f2*[1:N])]; Phi=toeplitz(Phi);
Phi=Phi+1e-10*eye(size(Phi)); %to stabilize the Cholosky factorization
%
% Form the matrices Sn' and Phi'.
I=eye(Lh1);
J=hankel([zeros(Lh1-1,1); 1]);
if rem(N+1,2)==1 J=J(2:end,:); end
E=[I; J]; Phi1=E'*Phi*E;
S1=[];
for n=1:N+1
    S1=[S1; E'*S(:,:,n)*E];
end

%
%   Add tail constraint to reduce PAR
%X=zeros(size(Phi1)); X=diag(X); 
X=zeros(Lh1,1); X(1:end-M)=eta; Phi1=Phi1+diag(X);
% 
%   Iterative lease-squares optimization
C=chol(Phi1);   % Choloskey factorization
for kk=1:100
    B=kron(eye(N+1),h1')*S1;
    D=Gamma*[B; C];
    u=zeros(N+1+Lh1,1);
    u(1)=1;
    u=Gamma*u;
    h1=(h1+(D'*D)\(D'*u))/2;
end
h=E*h1;

%
% SQUARE-ROOT RAISED-COSINE PULSE: h=sr_cos_p(N,L,alpha)
% This function generates a square-root raised-cosine pulse of length N+1. 
% There are L samples per symbol period.
% alpha is the roll-off factor.
%
function h=sr_cos_p(N,L,alpha)
t=[-N/2:1:N/2]/L;
h=zeros(size(t));
for k=1:length(t)
    if t(k)==0
        h(k)=1-alpha+4*alpha/pi;
    elseif (t(k)==(-1/(4*alpha)))|(t(k)==(1/(4*alpha)))
        h(k)=(alpha/sqrt(2))*((1+2/pi)*sin(pi/4/alpha)+(1-2/pi)*cos(pi/4/alpha));
    else
        h(k)=(sin(pi*(1-alpha)*t(k))+4*alpha*t(k)*cos(pi*(1+alpha)*t(k)))/(pi*t(k)*(1-(4*alpha*t(k))^2));
    end
end
h=h'/sqrt(L);