phasenoise_AWGN.m

   1 function phasenoise_AWGN(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.25;
  11   end
  12
  13   plotted = 0;
  14
  15   Rsym = 2.5e10; % symbol rate (sym/sec)
  16
  17   span = 6; % filter span
  18   sps = 4; % samples per symbol
  19
  20   fs = Rsym * sps; % sampling freq (Hz)
  21   Tsamp = 1 / fs;
  22
  23   t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs)';
  24
  25
  26   EbN0_db = 0:0.2:14;
  27   EbN0 = 10 .^ (EbN0_db ./ 10);
  28
  29   Es = 1;
  30   Eb = Es / log2(M);
  31   N0 = Eb ./ EbN0;
  32
  33   EsN0 = EbN0 .* log2(M);
  34   EsN0_db = 10 .* log10(EsN0);
  35
  36   plotlen = length(EbN0);
  37
  38   ber = zeros(1, plotlen);
  39
  40
  41   data = randi([0 M - 1], numSymbs, 1);
  42   modData = pskmod(data, M, 0, 'gray');
  43
  44   x = txFilter(modData, rolloff, span, sps);
  45
  46   linewidthTx = 0;%1e5; % Hz
  47   linewidthLO = 1e6; % Hz
  48   %%linewidthTx = Rsym * 1e-4; % Hz
  49   %%linewidthLO = Rsym * 1e-3; % Hz
  50
  51   totalIterations = 1;
  52
  53   for iter = 1:totalIterations
  54     [xPN, pTxLO] = phaseNoise(x, linewidthTx, linewidthLO, Tsamp);
  55
  56     for i = 1:plotlen
  57       snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps);
  58       noiseEnergy = 10 ^ (-snr / 10);
  59
  60       y = awgn(xPN, snr, 'measured');
  61
  62       r = rxFilter(y, rolloff, span, sps);
  63       %% normalize energy
  64       %r = normalizeEnergy(r, numSymbs, 1 + noiseEnergy);
  65
  66       [rPhaseEq, phiests] = phaseNoiseCorr(r, M, 40 * sps);
  67
  68       rSampled = rPhaseEq(sps*span/2+1:sps:(numSymbs + span/2) * sps);
  69       demodData = pskdemod(rSampled, M, 0, 'gray')';
  70
  71       %%[bitErrors, ber(i)] = biterr(data, demodData);
  72       [zzz, thisBER] = biterr(data, demodData);
  73       ber(i) = ber(i) + thisBER / totalIterations;
  74
  75
  76       if plotted == 0 && EbN0_db(i) >6 && ber(i) > 1e-1
  77         plotted=1
  78         figure(1234);
  79         plot(-phiests);
  80         hold on;
  81         plot(pTxLO);
  82         legend('estimate', 'actual');
  83         hold off;
  84
  85         figure(100);
  86         %plot(t(1:length(x)), real(x));
  87         %%plot(t, real(x(1:length(t))));
  88         length(t)
  89         length(x)
  90         hold on
  91         %%plot(t(1:length(xPhaseNoise)), real(xPhaseNoise));
  92
  93         plot(t, real(r), 'g')
  94
  95         sampledTimes = t(sps*span/2+1:sps:(numSymbs+span/2)*sps);
  96
  97         plot(sampledTimes, real(rSampled), 'x')
  98         hold off
  99
 100       end
 101
 102     end
 103   end
 104
 105
 106   figure(1);
 107   clf;
 108
 109   %% Plot simulated results
 110   semilogy(EbN0_db, ber, 'r', 'LineWidth', 2);
 111   hold on;
 112
 113   theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1);
 114   legend({'Simulated phase noise', 'Without phase noise'}, 'Location', 'southwest');
 115
 116   title(strcat(num2str(M), '-PSK with phase noise and correction'));
 117   grid on;
 118   xlabel('$E_b/N_0$ (dB)');
 119   ylabel('BER');
 120
 121   formatFigure;
 122 end