function SNRvsBER_BPSKdiscrete(numBits)
  if nargin < 1
    numBits = 1000;
  end

  if isOctave()
    pkg load communications
  end

  %% SNR: Eb/N0
  snr_db = 0:0.2:10;
  snr = 10 .^ (snr_db ./ 10);

  %% Eb is assumed to be 1
  %% noise N ~ Gaussian(mean=0, variance=N0/2)
  N0 = 1 ./ snr;

  ber = zeros(1, length(snr));

  data = randomPSK(2, numBits);

  for i = 1:length(snr)
    noisyData = discreteAWGN(data, N0(i));
    decodedData = decodeNoisyBPSK(noisyData);

    bitErrors = sum(data ~= decodedData);
    ber(i) = bitErrors / numBits;
  end

  fig1 = figure(1);
  clf;

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

  %% Plot theoretical curve
  %% BPSK: bit error when noise N > sqrt(Eb)
  %%   Pr(N > sqrt(Eb))
  %% = Pr(Z > sqrt(Eb) / sqrt(N0/2))
  %% = Q(sqrt(2 * SNR))
  ber_th = qfunc(sqrt(2 * snr));
  semilogy(snr_db, ber_th, 'b', 'LineWidth', 1);

  grid on;
  xlabel('$E_b/N_0$ (dB)');
  ylabel('BER');
  legend('Simulated', '$Q\left(\sqrt{2 E_b / N_0}\right)$');

  formatFigure;
  saveas(gcf, strcat('BPSK_BER_SNR_', num2str(numBits), '.svg'));

end