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

  if isOctave()
    pkg load communications
  end

  %% https://www.mathworks.com/help/signal/ref/rcosdesign.html
  %% https://www.mathworks.com/help/comm/ug/pulse-shaping-using-a-raised-cosine-filter.html
  span = 6; % filter span
  sps = 4;

  rrcFilter = rcosdesign(rolloff, span, sps, 'sqrt');

  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], numSymbs, 1);
  modData = pskmod(data, M, 0, 'gray');

  txSig = upfirdn(modData, rrcFilter, sps);

  for i = 1:plotlen
    snr = EbN0_db(i) + 10 * log10(log2(M));% - 10 * log10(sps); % why sps?
    rxSig = awgn(txSig, snr);

    rxFilt = upfirdn(rxSig, rrcFilter, 1, sps);
    rxFilt = rxFilt(span + 1 : end - span); % remove filter delay

    demodData = pskdemod(rxFilt, M, 0, 'gray');

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

  fig1 = figure(1);
  clf;

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

  %% Plot theoretical curve
  %% BPSK: bit error when noise Nr > sqrt(Eb)
  %%   Pr(Nr > sqrt(Eb))
  %% = Pr(Z > sqrt(Eb) / sqrt(N0/2))
  %%
  %% QPSK = 2 BPSKs, one real and one imaginary, each with one bit
  %% so BER is the same as BPSK (assuming Gray code)
  if M == 2 || M == 4
    ber_th = qfunc(sqrt(2 * EbN0));
    semilogy(EbN0_db, ber_th, 'b', 'LineWidth', 1);
    legend('Simulated RRC', 'Discrete');
  else
    %% Approximation: J.G. Proakis and M. Salehi, 2000, Contemporary
    %%                Communication Systems using MATLAB (Equations
    %%                7.3.18 and 7.3.19), Brooks/Cole.
    ber_ap = 2 * qfunc(sqrt(EbN0 * log2(M) * 2) * sin(pi / M)) / log2(M);
    semilogy(EbN0_db, ber_ap, 'b', 'LineWidth', 1);
    legend('Simulated RRC', 'Discrete');
  end

  title(strcat(num2str(M), '-PSK RRC 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'));

  %scatterplot(rxFilt);
  %eyediagram(rxFilt, sps);

end