phasenoise_AWGN.m

   1 numSymbs = 5e4;
   2 M = 4;
   3 rolloff = 0.5;
   4
   5 Rsym = 2.5e10; % symbol rate (sym/sec)
   6
   7 span = 6; % filter span
   8 sps = 8; % 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 = 1;
  52 avgSa = 40;
  53
  54 TsampOrig = Tsamp;
  55
  56 for it = 1 : iterations
  57   [xPSKpn, pTxLoPSK] = phaseNoise(xPSK, linewidthTx, linewidthLO, Tsamp);
  58   [xDEPSKpn, pTxLoDEPSK] = phaseNoise(xDEPSK, linewidthTx, linewidthLO, Tsamp);
  59   [xDPSKpn, pTxLoDPSK] = phaseNoise(xDPSK, linewidthTx, linewidthLO, Tsamp);
  60
  61   for i = 1:plotlen
  62     Tsamp = TsampOrig;
  63     sps = 8;
  64
  65     snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps);
  66     noiseEnergy = 10 ^ (-snr / 10);
  67
  68     yPSK = awgn(xPSKpn, snr, 'measured');
  69     yDEPSK = awgn(xDEPSKpn, snr, 'measured');
  70     yDPSK = awgn(xDPSKpn, snr, 'measured');
  71
  72     rPSK = rxFilter(yPSK, rolloff, span, sps);
  73     rDEPSK = rxFilter(yDEPSK, rolloff, span, sps);
  74     rDPSK = rxFilter(yDPSK, rolloff, span, sps);
  75
  76     sps = 2;
  77     Tsamp = TsampOrig * 4;
  78
  79     rPSKSamp = rPSK(1:2:end);
  80     rDEPSKSamp = rDEPSK(1:2:end);
  81     rDPSKSamp = rDPSK(1:2:end);
  82
  83     [rPSKSampEq, phiestsPSK] = phaseNoiseCorr(rPSKSamp, M, pi/M, avgSa);
  84     [rDEPSKSampEq, phiestsDEPSK] = phaseNoiseCorr(rDEPSKSamp, M, 0, avgSa);
  85
  86     demodPSK = pskdemod(rPSKSampEq, M, pi/M, 'gray').';
  87     %% The decoding method described in Taylor (2009)
  88     %% works on the complex symbols, i.e. after taking
  89     %% the nearest symbol in the constellation, but before
  90     %% converting them back to integers/bits.
  91     %% MATLAB's pskdemod() does not provide this intermediate
  92     %% result, so to be lazy, a pskmod() call is performed
  93     %% to obtain the complex symbols.
  94     demodDEPSK = pskdemod(rDEPSKSampEq, M, 0, 'gray').';
  95     remodDEPSK = pskmod(demodDEPSK, M, 0, 'gray');
  96     delayed = [1; remodDEPSK(1:end-1)];
  97     demodDEPSK = pskdemod(remodDEPSK .* conj(delayed), M, 0, 'gray');
  98
  99     demodDPSK = dpskdemod(rDPSKSamp, M, pi/M, 'gray');
 100
 101     [~, ber] = biterr(data, demodPSK);
 102     berPSK(i) = berPSK(i) + ber / iterations;
 103     [~, ber] = biterr(data, demodDEPSK);
 104     berDEPSK(i) = berDEPSK(i) + ber / iterations;
 105     [~, ber] = biterr(data, demodDPSK);
 106     berDPSK(i) = berDPSK(i) + ber / iterations;
 107
 108     if EbN0_db(i) == 8 && it == 1
 109       figure(1234);
 110       plot(repelem(-phiestsPSK, 8));
 111       hold on;
 112       plot(pTxLoPSK);
 113       legend('estimate', 'actual');
 114       hold off;
 115
 116       figure(1);
 117       scatterplot(rPSKSampEq);
 118       title('rPSKSampEq');
 119     end
 120
 121   end
 122 end
 123
 124
 125 figure(1);
 126 clf;
 127
 128 %% Plot simulated results
 129 semilogy(EbN0_db, berPSK, 'r', 'LineWidth', 1.5);
 130 hold on;
 131 semilogy(EbN0_db, berDEPSK, 'c', 'LineWidth', 2);
 132 semilogy(EbN0_db, berDPSK, 'Color', [0, 0.6, 0], 'LineWidth', 2.5);
 133
 134 theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1);
 135 DEPSKTheoretical = berawgn(EbN0_db, 'psk', M, 'diff');
 136 semilogy(EbN0_db, DEPSKTheoretical, 'Color', [1, 0.6, 0], 'LineWidth', 1);
 137 DPSKTheoretical = berawgn(EbN0_db, 'dpsk', M);
 138 semilogy(EbN0_db, DPSKTheoretical, 'm', 'LineWidth', 1);
 139
 140 legend({'PSK with Viterbi-Viterbi', ...
 141         'DEPSK with Viterbi-Viterbi', ...
 142         'DPSK', ...
 143         'Theoretical PSK over AWGN', ...
 144         'Theoretical DEPSK over AWGN', ...
 145         'Theoretical DPSK over AWGN'}, ...
 146        'Location', 'southwest');
 147
 148 title({'QPSK with phase nosie and correction', ...
 149        strcat('$10^{', num2str(log10(numSymbs * log2(M))), ...
 150               '}$~bits, LO~', ...
 151               num2str(linewidthLO / 1e6), '~MHz, blocksize~', ...
 152               num2str(avgSa), '~Sa')});
 153 grid on;
 154 xlabel('$E_b/N_0$ (dB)');
 155 ylabel('BER');
 156
 157 formatFigure;