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