function phasenoise_AWGN(rolloff, M, numSymbs)
  %% Set defaults for inputs
  if nargin < 3
    numSymbs = 1000;
  end
  if nargin < 2
    M = 2;
  end
  if nargin < 1
    rolloff = 0.25;
  end

  plotted = 0;

  Rsym = 2.5e10; % symbol rate (sym/sec)

  span = 6; % filter span
  sps = 4; % samples per symbol

  fs = Rsym * sps; % sampling freq (Hz)
  Tsamp = 1 / fs;

  t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs)';


  EbN0_db = 0:0.2:14;
  EbN0 = 10 .^ (EbN0_db ./ 10);

  Es = 1;
  Eb = Es / log2(M);
  N0 = Eb ./ EbN0;

  EsN0 = EbN0 .* log2(M);
  EsN0_db = 10 .* log10(EsN0);

  plotlen = length(EbN0);

  ber = zeros(1, plotlen);


  data = randi([0 M - 1], numSymbs, 1);
  modData = pskmod(data, M, 0, 'gray');

  x = txFilter(modData, rolloff, span, sps);

  linewidthTx = 0;%1e5; % Hz
  linewidthLO = 1e6; % Hz
  %%linewidthTx = Rsym * 1e-4; % Hz
  %%linewidthLO = Rsym * 1e-3; % Hz

  totalIterations = 1;

  for iter = 1:totalIterations
    [xPN, pTxLO] = phaseNoise(x, linewidthTx, linewidthLO, Tsamp);

    for i = 1:plotlen
      snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps);
      noiseEnergy = 10 ^ (-snr / 10);

      y = awgn(xPN, snr, 'measured');

      r = rxFilter(y, rolloff, span, sps);
      %% normalize energy
      %r = normalizeEnergy(r, numSymbs, 1 + noiseEnergy);

      [rPhaseEq, phiests] = phaseNoiseCorr(r, M, 40 * sps);

      rSampled = rPhaseEq(sps*span/2+1:sps:(numSymbs + span/2) * sps);
      demodData = pskdemod(rSampled, M, 0, 'gray')';

      %%[bitErrors, ber(i)] = biterr(data, demodData);
      [zzz, thisBER] = biterr(data, demodData);
      ber(i) = ber(i) + thisBER / totalIterations;


      if plotted == 0 && EbN0_db(i) >6 && ber(i) > 1e-1
        plotted=1
        figure(1234);
        plot(-phiests);
        hold on;
        plot(pTxLO);
        legend('estimate', 'actual');
        hold off;

        figure(100);
        %plot(t(1:length(x)), real(x));
        %%plot(t, real(x(1:length(t))));
        length(t)
        length(x)
        hold on
        %%plot(t(1:length(xPhaseNoise)), real(xPhaseNoise));

        plot(t, real(r), 'g')

        sampledTimes = t(sps*span/2+1:sps:(numSymbs+span/2)*sps);

        plot(sampledTimes, real(rSampled), 'x')
        hold off

      end

    end
  end


  figure(1);
  clf;

  %% Plot simulated results
  semilogy(EbN0_db, ber, 'r', 'LineWidth', 2);
  hold on;

  theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1);
  legend({'Simulated phase noise', 'Without phase noise'}, 'Location', 'southwest');

  title(strcat(num2str(M), '-PSK with phase noise and correction'));
  grid on;
  xlabel('$E_b/N_0$ (dB)');
  ylabel('BER');

  formatFigure;
end