function [rPhaseEq, phiests] = phaseNoiseCorr(r, M, phoffset, blocksize)
  %% phase noise correction
  phiests = zeros(1, length(r));
  rPhaseEq = zeros(1, length(r));
  for l = 1:blocksize:length(r)
    block = r(l : min(l + blocksize - 1, length(r)));

    sum_M = sum(block .^ M);
    phi_est = angle(sum_M .* exp(j * M * phoffset)) / M;

    if l > 1
      %% phase unwrapping
      phi_prev = phiests(l - 1);
      m = floor(0.5 + (phi_prev - phi_est) * M / (2 * pi));
      phi_est = phi_est + m * 2 * pi / M;
    end

    phiests(l:min(l+blocksize-1, length(r))) = phi_est * ones(1, min(blocksize, length(r) - l+1));

    block = block .* exp(j * -phi_est);
    rPhaseEq(l : min(l + blocksize-1, length(r))) = block;
  end
end