rNyquistM.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% root-Nyquist (M) filter design                        %
% Source:                                               %
% https://my.ece.utah.edu/~farhang/Nyquist_M_r1.pdf     %
% Edited by Wojciech SP5WWP to work with rcosdesign()   %
%                                                       %
% parameters:                                           %
% N: filter order (filter length = N+1)                 %
% M: number of samples per symbol period                %
% alpha: rolloff 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)                                     %
% itns: Number of iterations for the least-squares      %
% optimization                                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function h=rNyquistM(N,M,alpha,gmaZ,gmaT,eta,itns)
	%%% Set up 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)];
	Gamma2=Gamma.^2; Gamma2=diag(Gamma2);
	%%% Initial filter %%%
	h=rcosdesign(alpha,N/M,M);
	if rem(N+1,2)==0 h1=h(1:(N+1)/2);else h1=h(1:N/2+1); end
	Lh1=length(h1);
	%%% Set up 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 up the matrix Phi %%%
	Phi=zeros(N+1,N+1); fo=(1/2/M)*(1+alpha);
	Phi=[1-2*fo -2*fo*sinc(2*fo*[1:N])]; Phi=toeplitz(Phi);
	Phi=Phi+1e-10*eye(size(Phi)); %to stabilize Chol fac.
	%%% Form the matrices S′ 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(Lh1,1); X(1:end-M)=eta; Phi1=Phi1+diag(X);
	%%% Iterative least-squares optimization %%%
	C=chol(Phi1); % Cholesky factorization
        h1=h1';
 	for kk=1:itns
 		B=kron(eye(N+1),h1')*S1; D=[B; C]; u=[1;zeros(N+Lh1,1)];
 		h1=(h1+inv(D'*Gamma2*D)*(D'*Gamma2*u))/2;
 	end
 	h=E*h1;
end