baseband.m

   1 function baseband(rolloff, M, numSymbs)
   2   %% Set defaults for inputs
   3   if nargin < 3
   4     numSymbs = 1000;
   5   end
   6   if nargin < 2
   7     M = 2;
   8   end
   9   if nargin < 1
  10     rolloff = 0.5;
  11   end
  12
  13
  14   %% https://www.mathworks.com/help/comm/examples/passband-modulation-with-adjacent-channel-interference.html
  15   Rsym = 1e6; % symbol rate (sym/sec)
  16
  17   span = 6; % filter span
  18   sps = 4; % samples per symbol
  19
  20   txFilter = comm.RaisedCosineTransmitFilter...
  21                  ('Shape', 'Square root', ...
  22                   'RolloffFactor', rolloff, ...
  23                   'FilterSpanInSymbols', span, ...
  24                   'OutputSamplesPerSymbol', sps);
  25   rxFilter = comm.RaisedCosineReceiveFilter...
  26                  ('Shape', 'Square root', ...
  27                   'RolloffFactor', rolloff, ...
  28                   'FilterSpanInSymbols', span, ...
  29                   'InputSamplesPerSymbol', sps, ...
  30                   'DecimationFactor', 1);
  31
  32   fs = Rsym * sps; % sampling freq (Hz)
  33
  34   t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs)';
  35
  36
  37   EbN0_db = 0:0.2:10;
  38   EbN0 = 10 .^ (EbN0_db ./ 10);
  39   Es = 1;
  40   Eb = Es / log2(M);
  41   N0 = Eb ./ EbN0;
  42
  43   EsN0 = EbN0 .* log2(M);
  44   EsN0_db = 10 .* log10(EsN0);
  45
  46   plotlen = length(EbN0);
  47   ber = zeros(1, plotlen);
  48
  49   data = randi([0 M - 1], numSymbs, 1);
  50   modData = pskmod(data, M, 0, 'gray');
  51
  52   xBaseband = txFilter([modData; zeros(span, 1)]);
  53
  54
  55
  56   for i = 1:plotlen
  57     snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps); % why sps?
  58     noiseEnergy = 10 ^ (-snr / 10);
  59
  60     yBaseband = awgn(xBaseband, snr, 'measured');
  61
  62     rBaseband = rxFilter([yBaseband; zeros(span, 1)]);
  63     %% truncate filter transients
  64     rBaseband = rBaseband(span * sps / 2 + 1 : end);
  65     %% normalize to unit energy
  66     rBasebandEnergy = sum(abs(rBaseband) .^ 2) / numSymbs;
  67     rBaseband = rBaseband .* sqrt((1 + noiseEnergy) / rBasebandEnergy);
  68
  69     rSampled = rBaseband(sps*span/2+1:sps:(numSymbs+span/2)*sps);
  70
  71     demodData = pskdemod(rSampled, M, 0, 'gray');
  72
  73     [bitErrors, ber(i)] = biterr(data, demodData);
  74   end
  75
  76   fig1 = figure(1);
  77   clf;
  78
  79   %% Plot simulated results
  80   semilogy(EbN0_db, ber, 'r', 'LineWidth', 2);
  81   hold on;
  82
  83   %% Plot theoretical curve
  84   %% BPSK: bit error when noise Nr > sqrt(Eb)
  85   %%   Pr(Nr > sqrt(Eb))
  86   %% = Pr(Z > sqrt(Eb) / sqrt(N0/2))
  87   %%
  88   %% QPSK = 2 BPSKs, one real and one imaginary, each with one bit
  89   %% so BER is the same as BPSK (assuming Gray code)
  90   if M == 2 || M == 4
  91     ber_th = qfunc(sqrt(2 * EbN0));
  92     semilogy(EbN0_db, ber_th, 'b', 'LineWidth', 1);
  93     legend('Simulated', 'Discrete');
  94   else
  95     %% Approximation: J.G. Proakis and M. Salehi, 2000, Contemporary
  96     %%                Communication Systems using MATLAB (Equations
  97     %%                7.3.18 and 7.3.19), Brooks/Cole.
  98     ber_ap = 2 * qfunc(sqrt(EbN0 * log2(M) * 2) * sin(pi / M)) / log2(M);
  99     semilogy(EbN0_db, ber_ap, 'b', 'LineWidth', 1);
 100     legend('Simulated', 'Discrete');
 101   end
 102
 103   title(strcat(num2str(M), '-PSK with Gray code'));
 104   grid on;
 105   xlabel('$E_b/N_0$ (dB)');
 106   ylabel('BER');
 107
 108   formatFigure;
 109   %saveas(gcf, strcat('BER_SNR_', num2str(M), 'PSK_', num2str(numSymbs), ...
 110   %                   '.svg'));
 111
 112   %scatterplot(rxFilt);
 113   %eyediagram(rxFilt, sps);
 114
 115 end