discretePSK_BER_SNR.m

   1 function discretePSK_BER_SNR(M, numBits)
   2   %% Set defaults for inputs
   3   if nargin < 2
   4     numBits = 1000;
   5   end
   6   if nargin < 1
   7     M = 2;
   8   end
   9
  10   if isOctave()
  11     pkg load communications
  12   end
  13
  14   EbN0_db = 0:0.2:10;
  15   EbN0 = 10 .^ (EbN0_db ./ 10);
  16
  17   Es = 1;
  18   Eb = Es / log2(M);
  19   N0 = Eb ./ EbN0;
  20
  21   EsN0 = EbN0 .* log2(M);
  22   EsN0_db = 10 .* log10(EsN0);
  23
  24   plotlen = length(EbN0);
  25
  26   ber = zeros(1, plotlen);
  27
  28   data = randi([0, M - 1], 1, numBits);
  29   txsig = pskmod(data, M, 0, 'gray');
  30
  31   for i = 1:plotlen
  32     rxsig = awgn(txsig, EsN0_db(i));
  33     demodData = pskdemod(rxsig, M, 0, 'gray');
  34
  35     [bitErrors, ber(i)] = biterr(data, demodData);
  36   end
  37
  38   fig1 = figure(1);
  39   clf;
  40
  41   %% Plot simulated results
  42   semilogy(EbN0_db, ber, 'r', 'LineWidth', 2);
  43   hold on;
  44
  45   %% Plot theoretical curve
  46   %% BPSK: bit error when noise Nr > sqrt(Eb)
  47   %%   Pr(Nr > sqrt(Eb))
  48   %% = Pr(Z > sqrt(Eb) / sqrt(N0/2))
  49   %%
  50   %% QPSK = 2 BPSKs, one real and one imaginary, each with one bit
  51   %% so BER is the same as BPSK (assuming Gray code)
  52   if M == 2 || M == 4
  53     ber_th = qfunc(sqrt(2 * EbN0));
  54     semilogy(EbN0_db, ber_th, 'b', 'LineWidth', 1);
  55     legend('Simulated', 'Theoretical');
  56   else
  57     %% Upper bound: R. Venkataramanan, Lent 2018, 3F4 Examples Paper 2
  58     %%              (Question 5), CUED.
  59     %% Approximation: J.G. Proakis and M. Salehi, 2000, Contemporary
  60     %%                Communication Systems using MATLAB (Equations
  61     %%                7.3.18 and 7.3.19), Brooks/Cole.
  62     ber_ub = 2 * qfunc(sqrt(EbN0 * log2(M)) * sin(pi / M));
  63     ber_ap = 2 * qfunc(sqrt(EbN0 * log2(M) * 2) * sin(pi / M)) / log2(M);
  64     semilogy(EbN0_db, ber_ub, 'b', 'LineWidth', 1);
  65     semilogy(EbN0_db, ber_ap, 'g', 'LineWidth', 1);
  66     legend('Simulated', 'Upper bound', 'Approximation');
  67   end
  68
  69   title(strcat(num2str(M), '-PSK with Gray code'));
  70   grid on;
  71   xlabel('$E_b/N_0$ (dB)');
  72   ylabel('BER');
  73
  74   formatFigure;
  75   saveas(gcf, strcat('BER_SNR_', num2str(M), 'PSK_', num2str(numBits), ...
  76                      '.svg'));
  77 end