phasenoise_AWGN.m

   1 numSymbs = 5e3;
   2 M = 4;
   3 rolloff = 0.5;
   4
   5 Rsym = 2.5e10; % symbol rate (sym/sec)
   6
   7 span = 6; % filter span
   8 sps = 4; % samples per symbol
   9
  10 fs = Rsym * sps; % sampling freq (Hz)
  11 Tsamp = 1 / fs;
  12
  13 t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).';
  14
  15
  16 EbN0_db = 0:0.2:14;
  17 EbN0 = 10 .^ (EbN0_db ./ 10);
  18
  19 Es = 1;
  20 Eb = Es / log2(M);
  21 N0 = Eb ./ EbN0;
  22
  23 EsN0 = EbN0 .* log2(M);
  24 EsN0_db = 10 .* log10(EsN0);
  25
  26 plotlen = length(EbN0);
  27
  28 berPSK = zeros(1, plotlen);
  29 berDEPSK = zeros(1, plotlen);
  30 berDPSK = zeros(1, plotlen);
  31
  32 data = randi([0 M - 1], numSymbs, 1);
  33
  34 pskSym = pskmod(data, M, pi / M, 'gray');
  35 %% DEPSK: Part VII, M.G. Taylor (2009)
  36 depskSym = pskmod(data, M, 0, 'gray');
  37 for i = 2:numSymbs
  38   depskSym(i) = depskSym(i) * depskSym(i-1);
  39 end
  40
  41 dpskSym = dpskmod(data, M, pi / M, 'gray');
  42
  43 xPSK = txFilter(pskSym, rolloff, span, sps);
  44 xDEPSK = txFilter(depskSym, rolloff, span, sps);
  45 xDPSK = txFilter(dpskSym, rolloff, span, sps);
  46
  47 linewidthTx = 0; % Hz
  48 linewidthLO = 5e6; % Hz
  49 %linewidthLO = Rsym * 1e-3;
  50
  51 iterations = 25;
  52 avgSa = 40;
  53
  54 for it = 1 : iterations
  55   [xPSKpn, pTxLoPSK] = phaseNoise(xPSK, linewidthTx, linewidthLO, Tsamp);
  56   [xDEPSKpn, pTxLoDEPSK] = phaseNoise(xDEPSK, linewidthTx, linewidthLO, Tsamp);
  57   [xDPSKpn, pTxLoDPSK] = phaseNoise(xDPSK, linewidthTx, linewidthLO, Tsamp);
  58
  59   for i = 1:plotlen
  60     snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps);
  61     noiseEnergy = 10 ^ (-snr / 10);
  62
  63     yPSK = awgn(xPSKpn, snr, 'measured');
  64     yDEPSK = awgn(xDEPSKpn, snr, 'measured');
  65     yDPSK = awgn(xDPSKpn, snr, 'measured');
  66
  67     rPSK = rxFilter(yPSK, rolloff, span, sps);
  68     rDEPSK = rxFilter(yDEPSK, rolloff, span, sps);
  69     rDPSK = rxFilter(yDPSK, rolloff, span, sps);
  70
  71     rPSKSamp = rPSK(sps*span/2+1:sps:(numSymbs+span/2)*sps);
  72     rDEPSKSamp = rDEPSK(sps*span/2+1:sps:(numSymbs+span/2)*sps);
  73     rDPSKSamp = rDPSK(sps*span/2+1:sps:(numSymbs+span/2)*sps);
  74
  75     [rPSKSampEq, phiestsPSK] = phaseNoiseCorr(rPSKSamp, M, pi/M, avgSa);
  76     [rDEPSKSampEq, phiestsDEPSK] = phaseNoiseCorr(rDEPSKSamp, M, 0, avgSa);
  77
  78     demodPSK = pskdemod(rPSKSampEq, M, pi/M, 'gray').';
  79     %% The decoding method described in Taylor (2009)
  80     %% works on the complex symbols, i.e. after taking
  81     %% the nearest symbol in the constellation, but before
  82     %% converting them back to integers/bits.
  83     %% MATLAB's pskdemod() does not provide this intermediate
  84     %% result, so to be lazy, a pskmod() call is performed
  85     %% to obtain the complex symbols.
  86     demodDEPSK = pskdemod(rDEPSKSampEq, M, 0, 'gray').';
  87     remodDEPSK = pskmod(demodDEPSK, M, 0, 'gray');
  88     delayed = [1; remodDEPSK(1:end-1)];
  89     demodDEPSK = pskdemod(remodDEPSK .* conj(delayed), M, 0, 'gray');
  90
  91     demodDPSK = dpskdemod(rDPSKSamp, M, pi/M, 'gray');
  92
  93     [~, ber] = biterr(data, demodPSK);
  94     berPSK(i) = berPSK(i) + ber / iterations;
  95     [~, ber] = biterr(data, demodDEPSK);
  96     berDEPSK(i) = berDEPSK(i) + ber / iterations;
  97     [~, ber] = biterr(data, demodDPSK);
  98     berDPSK(i) = berDPSK(i) + ber / iterations;
  99
 100     if EbN0_db(i) == 8 && it == 1
 101       figure(1234);
 102       plot(repelem(-phiestsPSK, sps));
 103       hold on;
 104       plot(pTxLoPSK);
 105       legend('estimate', 'actual');
 106       hold off;
 107
 108       figure(1);
 109       scatterplot(rPSKSampEq);
 110       title('rPSKSampEq');
 111     end
 112
 113   end
 114 end
 115
 116
 117 figure(1);
 118 clf;
 119
 120 %% Plot simulated results
 121 semilogy(EbN0_db, berPSK, 'r', 'LineWidth', 1.5);
 122 hold on;
 123 semilogy(EbN0_db, berDEPSK, 'c', 'LineWidth', 2);
 124 semilogy(EbN0_db, berDPSK, 'Color', [0, 0.6, 0], 'LineWidth', 2.5);
 125
 126 theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1);
 127 DEPSKTheoretical = berawgn(EbN0_db, 'psk', M, 'diff');
 128 semilogy(EbN0_db, DEPSKTheoretical, 'Color', [1, 0.6, 0], 'LineWidth', 1);
 129 DPSKTheoretical = berawgn(EbN0_db, 'dpsk', M);
 130 semilogy(EbN0_db, DPSKTheoretical, 'm', 'LineWidth', 1);
 131
 132 legend({'PSK with Viterbi-Viterbi', ...
 133         'DEPSK with Viterbi-Viterbi', ...
 134         'DPSK', ...
 135         'Theoretical PSK over AWGN', ...
 136         'Theoretical DEPSK over AWGN', ...
 137         'Theoretical DPSK over AWGN'}, ...
 138        'Location', 'southwest');
 139
 140 title({'QPSK with phase nosie and correction', ...
 141        strcat('$10^{', num2str(log10(numSymbs * log2(M))), ...
 142               '}$~bits, LO~', ...
 143               num2str(linewidthLO / 1e6), '~MHz, blocksize~', ...
 144               num2str(avgSa), '~Sa')});
 145 grid on;
 146 xlabel('$E_b/N_0$ (dB)');
 147 ylabel('BER');
 148
 149 formatFigure;