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

  EbN0_db = 0:0.2:10;
  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], 1, numSymbs);
  txsig = pskmod(data, M, 0, 'gray');

  for i = 1:plotlen
    rxsig = awgn(txsig, EsN0_db(i));
    demodData = pskdemod(rxsig, M, 0, 'gray');

    [bitErrors, ber(i)] = biterr(data, demodData);
  end

  figure(1);
  clf;

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

  theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1);

  if M == 2 || M == 4
    legend('Simulated', 'Theoretical');
  else
    legend('Simulated', 'Approximation');
  end

  title(strcat(num2str(M), '-PSK with Gray code'));
  grid on;
  xlabel('$E_b/N_0$ (dB)');
  ylabel('BER');

  formatFigure;
  %saveas(gcf, strcat('BER_SNR_', num2str(M), 'PSK_', num2str(numSymbs), ...
  %                   '.svg'));
end