function [ estimator ] = squaringtiming( s, N )
% Squaring timing phase recovery
% Input
%   s = digital signal
%   N = number of samples per symbol
% Output
%   estimator =  phase estimate for each symbol (between 0 and N),
%           non-integer values lie between two samples
x = 0; 
for i=1:length(s)
    x = x + s(i)*s(i)*exp(-j*2*pi*i/N);
end
estimator = -myarg( x  ) / (2*pi);
end

