function theoreticalPSK(EbN0_db, M, varargin)
  %% Plot theoretical curve
  EbN0 = 10 .^ (EbN0_db ./ 10);
  if M == 2 || M == 4
    %% 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)
    ber_th = qfunc(sqrt(2 * EbN0));
    semilogy(EbN0_db, ber_th, varargin{:});
  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, varargin{:});
  end
end