phasenoise_AWGN.m

   1 numSymbs = 5e5;
   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 ber = zeros(1, plotlen);
  29
  30 data = randi([0 M - 1], numSymbs, 1);
  31 %%modData = dpskmod(data, M, pi / M, 'gray');
  32 modData = dpskmod(data, M, 0, 'gray');
  33
  34 x = txFilter(modData, rolloff, span, sps);
  35
  36 linewidthTx = 0; % Hz
  37 %%linewidthLO = 1e6; % Hz
  38 linewidthLO = Rsym * 1e-3;
  39
  40
  41 avgSa = 40;
  42
  43
  44 [xPN, pTxLO] = phaseNoise(x, linewidthTx, linewidthLO, Tsamp);
  45
  46 for i = 1:plotlen
  47   snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps);
  48   noiseEnergy = 10 ^ (-snr / 10);
  49
  50   y = awgn(xPN, snr, 'measured');
  51
  52   r = rxFilter(y, rolloff, span, sps);
  53   %% normalize energy
  54   r = normalizeEnergy(r, numSymbs, 1 + noiseEnergy);
  55
  56   rSampled = r(sps*span/2+1:sps:(numSymbs+span/2)*sps);
  57   [rSaPhEq, phiests] = phaseNoiseCorr(rSampled, M, avgSa);
  58
  59   adaptiveFilterOut = adaptiveCMA(rSaPhEq.');
  60
  61   demodData = dpskdemod(rSaPhEq, M, 0, 'gray').';
  62   demodAdapt = dpskdemod(adaptiveFilterOut, M, 0, 'gray');
  63
  64   [~, ber(i)] = biterr(data, demodData);
  65
  66   if EbN0_db(i) == 8
  67     figure(1234);
  68     plot(repelem(-phiests, sps));
  69     hold on;
  70     plot(pTxLO);
  71     legend('estimate', 'actual');
  72     hold off;
  73
  74     figure(1);
  75     scatterplot(rSaPhEq);
  76     title('rSaPhEq');
  77   end
  78 end
  79
  80
  81 figure(1);
  82 clf;
  83
  84 %% Plot simulated results
  85 semilogy(EbN0_db, ber, 'r', 'LineWidth', 2);
  86 hold on;
  87
  88 theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1);
  89 legend({'Simulated phase noise + correction', 'Theoretical AWGN'}, ...
  90        'Location', 'southwest');
  91
  92
  93 title({'QPSK with phase nosie and correction', ...
  94        strcat(num2str(numSymbs * log2(M) / 1e3), '~kbit, LO~', ...
  95               num2str(linewidthLO / 1e6), '~MHz, Av~', num2str(avgSa), ...
  96               '~Sa')});
  97 grid on;
  98 xlabel('$E_b/N_0$ (dB)');
  99 ylabel('BER');
 100
 101 formatFigure;