passband.m

   1 function passband(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
  40   Es = 1;
  41   Eb = Es / log2(M);
  42   N0 = Eb ./ EbN0;
  43
  44   EsN0 = EbN0 .* log2(M);
  45   EsN0_db = 10 .* log10(EsN0);
  46
  47   plotlen = length(EbN0);
  48
  49   ber = zeros(1, plotlen);
  50
  51
  52
  53
  54   data = randi([0 M - 1], numSymbs, 1);
  55   modData = pskmod(data, M, 0, 'gray');
  56
  57   xBaseband = txFilter([modData; zeros(span, 1)]);
  58
  59
  60   %fc = 2.5e6; % Carrier freq (Hz)
  61   %carrier = sqrt(2) * exp(j * 2 * pi * fc * t);
  62
  63   %xPassbandIdeal = normalizeEnergy...
  64   %                   (real(xBaseband .* carrier(1:length(xBaseband))), numSymbs, 1);
  65
  66   %txLOFreq = [2.49e6, 2.5e6, 2.51e6];
  67   %%txLOEnergy = [0.05, 0.9, 0.05];
  68   %txLOEnergy = [0 1 0];
  69   txLOFreq = [1e8];
  70   txLOEnergy = [1];
  71
  72   carrier = zeros(length(t), 1);
  73   for i = 1 : length(txLOFreq)
  74     carrier = carrier + ...
  75               sqrt(2 * txLOEnergy(i)) * exp(j * 2 * pi * txLOFreq(i) * t);
  76   end
  77
  78   xPassband = normalizeEnergy...
  79                 (xBaseband .* carrier(1:length(xBaseband)), numSymbs, 1);
  80
  81   for i = 1:plotlen
  82     snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps); % why sps?
  83     noiseEnergy = 10 ^ (-snr / 10);
  84
  85     yPassband = awgn(xPassband, snr, 'measured');
  86
  87     rBaseband = rxFilter([yPassband .* conj(carrier(1:length(yPassband))); zeros(span * sps, 1)]);
  88     %% truncate filter transients
  89     rBaseband = rBaseband(span * sps / 2 + 1 : end);
  90     %% normalize energy
  91     rBaseband = normalizeEnergy(rBaseband, numSymbs, 1 + noiseEnergy);
  92
  93     rSampled = rBaseband(sps*span/2+1:sps:(numSymbs + span/2) * sps);
  94
  95     demodData = pskdemod(rSampled, M, 0, 'gray');
  96     [bitErrors, ber(i)] = biterr(data, demodData);
  97
  98   end
  99
 100   fig1 = figure(1);
 101   clf;
 102
 103   %% Plot simulated results
 104   semilogy(EbN0_db, ber, 'r', 'LineWidth', 2);
 105   hold on;
 106
 107   %% Plot theoretical curve
 108   %% BPSK: bit error when noise Nr > sqrt(Eb)
 109   %%   Pr(Nr > sqrt(Eb))
 110   %% = Pr(Z > sqrt(Eb) / sqrt(N0/2))
 111   %%
 112   %% QPSK = 2 BPSKs, one real and one imaginary, each with one bit
 113   %% so BER is the same as BPSK (assuming Gray code)
 114   if M == 2 || M == 4
 115     ber_th = qfunc(sqrt(2 * EbN0));
 116     semilogy(EbN0_db, ber_th, 'b', 'LineWidth', 1);
 117     legend('Simulated RRC', 'Discrete');
 118   else
 119     %% Approximation: J.G. Proakis and M. Salehi, 2000, Contemporary
 120     %%                Communication Systems using MATLAB (Equations
 121     %%                7.3.18 and 7.3.19), Brooks/Cole.
 122     ber_ap = 2 * qfunc(sqrt(EbN0 * log2(M) * 2) * sin(pi / M)) / log2(M);
 123     semilogy(EbN0_db, ber_ap, 'b', 'LineWidth', 1);
 124     legend('Simulated RRC', 'Discrete');
 125   end
 126
 127   title(strcat(num2str(M), '-PSK RRC with Gray code'));
 128   grid on;
 129   xlabel('$E_b/N_0$ (dB)');
 130   ylabel('BER');
 131
 132   formatFigure;
 133   %saveas(gcf, strcat('BER_SNR_', num2str(M), 'PSK_', num2str(numSymbs), ...
 134   %                   '.svg'));
 135
 136   %scatterplot(rxFilt);
 137   %eyediagram(rxFilt, sps);
 138
 139 end
 140
 141
 142 function y = normalizeEnergy(x, numSymbs, e)
 143   energy = sum(abs(x) .^ 2) / numSymbs;
 144   y = x * sqrt(e / energy);
 145 end