function [estSFO estCFO numRoots] = OptPerSCEst(ZMat, DMat, HVec, RMat, N, Ng, symbArr, A, ...
noiseVar, maxEps, maxDelta, CSIFlag)
TAY_APPRX_SIZE = 3;
Ns = N + Ng;
idxA = A+(N/2)+1;
%idxA = A;
betaVec = zeros(size(A.'));
if (CSIFlag == 1)
alphaVec = (N-1 + 2.*(symbArr.*Ns + Ng))./N;
thetaMat = conj(ZMat(idxA,:)).*(DMat(idxA,:).*repmat(HVec(idxA),1,length(symbArr)));
Hest = HVec(idxA);
else
alphaVec = 2.*symbArr(1:end-1).*Ns./N;
thetaMat = conj(RMat);
%%% estimate the channel
YMat = ZMat(idxA,:).*conj(DMat(idxA,:));
expMat = zeros(length(A),length(symbArr));
for (kk = 1 : length(A))
for (ll = symbArr)
expMat(kk,ll) = exp(-1i*pi*(N-1+2*ll*Ns+2*Ng)*betaVec(kk)/N);
end
end
Hest = (1/length(symbArr))*(sum(YMat.*expMat,2));
end
maxVal = (maxEps + A(end)*maxDelta);
numRoots = [];
for (kk = 1:length(A))
thetaVec = thetaMat(kk,:); % vector per SC.
%%% Calculating the polynom
polyArr = zeros(1,2*(TAY_APPRX_SIZE + 1));
% Calc even powers
for (n = 0 : TAY_APPRX_SIZE)
nFactCoeff = ((-1)^n)/factorial(2*n);
polyArr(2*n+1) = nFactCoeff*imag(thetaVec)*((alphaVec*pi).^(2*n+1)).';
end
% Calc odd powers
for (n = 0 : TAY_APPRX_SIZE)
nFactCoeff = ((-1)^n)/factorial(2*n+1);
polyArr(2*n+2) = nFactCoeff*real(thetaVec)*((alphaVec*pi).^(2*n+2)).';
end
R = roots(fliplr(polyArr));
% Extract real roots.
candidates = [];
for (k = 1 : length(R))
if (isreal(R(k)) == 1)
if ((R(k) >= -maxVal) && (R(k) <= maxVal))
candidates = [candidates R(k)];
end
end
end
numRoots = [numRoots length(candidates)];
ResArr = zeros(size(candidates));
for (nn = 1 : length(candidates))
expVec = exp(1i.*pi.*alphaVec.*candidates(nn));
ResArr(nn) = sum(real(expVec.*thetaVec));
end
[~, idx] = max(ResArr);
if (isempty(idx))
est = 0;
else
est = candidates(idx);
end
NRarr = zeros(1,3);
X0arr = [-maxEps 0 maxEps];
% NRarr = zeros(1,2);
% X0arr = 0.5*[-maxEps maxEps];
% NRarr = zeros(1);
% X0arr = 0;
for (mm =1 : length(X0arr))
NRarr(mm) = NewtonRaphsonPoly(polyArr, X0arr(mm), 5);
end
% x = NewtonRaphsonPoly(polyArr, 0, 20);
% y = LaguerrePoly(polyArr, 0, 10);
% tmp = abs(est-y);
% if ((y <= -maxVal) || (y >= maxVal) || (isreal(y) == 0))
% y = 0;
% end
NRarr = unique(NRarr);
NRarr = NRarr(abs(NRarr) < maxVal);
tmpArr = zeros(size(NRarr));
for (nn = 1 : length(NRarr))
expVec = exp(1i.*pi.*alphaVec.*NRarr(nn));
tmpArr(nn) = sum(real(expVec.*thetaVec));
end
[~, idx] = max(tmpArr);
tmp = NRarr(idx);
if (isempty(idx))
tmp = 0;
end
betaVec(kk) = tmp;
%betaVec(kk) = est;
end
[estSFO estCFO] = BlueBasedOnBeta(Hest, N, Ng, symbArr, A, noiseVar, betaVec, 0);
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
