CD_AWGN.m

   1 numSymbs = 10000;
   2 M = 4;
   3
   4 Rsym = 2.5e10; % symbol rate (sym/sec)
   5
   6 rolloff = 0.25;
   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 EbN0_db = 0:0.2:14;
  16 EbN0 = 10 .^ (EbN0_db ./ 10);
  17
  18 Es = 1;
  19 Eb = Es / log2(M);
  20 N0 = Eb ./ EbN0;
  21
  22 EsN0 = EbN0 .* log2(M);
  23 EsN0_db = 10 .* log10(EsN0);
  24
  25 plotlen = length(EbN0);
  26
  27 ber = zeros(1, plotlen);
  28 berNoComp = zeros(1, plotlen);
  29 berAdapt = zeros(1, plotlen);
  30 berMatlabAdapt = zeros(1, plotlen);
  31
  32 data = randi([0 M - 1], numSymbs, 1);
  33 modData = pskmod(data, M, pi / M, 'gray');
  34 x = txFilter(modData, rolloff, span, sps);
  35
  36 %% Simulate chromatic dispersion
  37 D = 17; % ps / (nm km)
  38 lambda = 1550; % nm
  39 z = 10; % km
  40
  41 xCD = chromaticDispersion(x, D, lambda, z, Tsamp);
  42 xCD = normalizeEnergy(xCD, numSymbs, 1);
  43
  44 for i = 1:plotlen
  45   snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps);
  46   noiseEnergy = 10 ^ (-snr / 10);
  47
  48   y = awgn(xCD, snr, 'measured');
  49
  50   yCDComp = CDCompensation(y, D, lambda, z, Tsamp);
  51
  52   r = rxFilter(yCDComp, rolloff, span, sps);
  53   rNoComp = rxFilter(y, rolloff, span, sps);
  54   %% normalize energy
  55   %r = normalizeEnergy(r, numSymbs, 1 + noiseEnergy);
  56
  57   rSampled = r(sps*span/2+1:sps:(numSymbs + span/2) * sps);
  58   rNoCompSampled = rNoComp(sps*span/2+1:sps:(numSymbs+span/2)*sps);
  59
  60   %% rotate rNoCompSampled to match original data
  61   theta = angle(-sum(rNoCompSampled .^ M)) / M;
  62   %% if theta approx +pi/M, wrap to -pi/M
  63   if abs(theta - pi / M) / (pi / M) < 0.1
  64     theta = -pi / M;
  65   end
  66   theta
  67   rNoCompSampled = rNoCompSampled .* exp(-j * theta);
  68
  69   %% adaptive filter
  70   adaptFilterOut = adaptiveCMA(rSampled);
  71
  72   demodData = pskdemod(rSampled, M, pi / M, 'gray');
  73   demodNoComp = pskdemod(rNoCompSampled, M, pi / M, 'gray');
  74   demodAdapt = pskdemod(adaptFilterOut, M, pi / M, 'gray');
  75   %%demodMatlabAdapt = pskdemod(matlabEq, M, pi / M, 'gray');
  76
  77   [bitErrors, ber(i)] = biterr(data, demodData);
  78   [bitErrors, berNoComp(i)] = biterr(data, demodNoComp);
  79   [~, berAdapt(i)] = biterr(data, demodAdapt);
  80   %%[~, berMatlabAdapt(i)] = biterr(data, demodMatlabAdapt);
  81
  82
  83   if EbN0_db(i) == 12
  84     figure(1);
  85     scatterplot(rSampled);
  86     title('Constellation after CD compensation');
  87     %%scatterplot(modData);
  88     %%title('Original constellation');
  89     scatterplot(rNoCompSampled);
  90     title('Constellation without CD compensation');
  91     scatterplot(adaptFilterOut);
  92     title('Constellation with CD compensation and adaptive filter');
  93     %scatterplot(matlabEq);
  94     %title('Matlab equalizer');
  95     ber(i)
  96     %berNoComp(i)
  97     berAdapt(i)
  98     berMatlabAdapt(i)
  99   end
 100
 101 end
 102
 103 figure(1);
 104 clf;
 105
 106 %% Plot simulated results
 107 semilogy(EbN0_db, ber, 'r', 'LineWidth', 2);
 108 hold on;
 109 %%semilogy(EbN0_db, berNoComp, 'g', 'LineWidth', 2);
 110 semilogy(EbN0_db, berAdapt, 'm', 'LineWidth', 1.4);
 111 %%%semilogy(EbN0_db, berMatlabAdapt, 'c', 'LineWidth', 1.4);
 112
 113 theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1);
 114 legend({'CD + AWGN + CD comp.', 'CD + AWGN + CD comp.~+ CMA', ...
 115         'Theoretical AWGN'}, 'Location', 'southwest');
 116
 117 title(strcat(num2str(M), '-PSK with chromatic dispersion and compensation'));
 118 grid on;
 119 xlabel('$E_b/N_0$ (dB)');
 120 ylabel('BER');
 121
 122 formatFigure;