chromaticDispersion1Signal.m

   1 M = 4;
   2 numSymbs = 5e5;
   3
   4 Rsym = 2.5e10; % symbol rate (sym/sec)
   5
   6 span = 6; % Tx/Rx filter span
   7 rolloff = 0.25; % Tx/Rx RRC rolloff
   8 sps = 2; % 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 data = randi([0 M - 1], numSymbs, 1);
  16 modData = pskmod(data, M, pi / M, 'gray');
  17 x = txFilter(modData, rolloff, span, sps);
  18
  19 x = normalizeEnergy(x, numSymbs*sps, 1);
  20
  21 %% Simulate chromatic dispersion
  22 D = 17; % ps / (nm km)
  23 lambda = 1550; % nm
  24 z = 5000 % km
  25
  26 [xCD, xCDkstart] = chromaticDispersion_FFT(x, D, lambda, z, Tsamp);
  27
  28 EbN0_db = 8;
  29 snr = EbN0_db + 10 * log10(log2(M)) - 10 * log10(sps);
  30
  31 %%y = awgn(xCD, snr, 'measured');
  32 y = xCD;
  33
  34 r = rxFilter(y, rolloff, span, sps);
  35 [rCDComp, CDCompkstart] = CDCompensation(r, D, lambda, z, Tsamp);
  36 rCDComp = normalizeEnergy(rCDComp, numSymbs*sps, 1);
  37
  38 rSampled = rCDComp(sps*span/2+1:sps:(numSymbs+span/2)*sps);
  39 rNoCompSa = r(sps*span/2+1:sps:(numSymbs+span/2)*sps);
  40
  41 %% if no CD comp, then rotate constellation. Use:
  42 theta = angle(-sum(rNoCompSa .^ M)) / M;
  43 %% if theta approx +pi/M, wrap to -pi/M
  44 if abs(theta - pi / M) / (pi / M) < 0.1
  45   theta = -pi / M;
  46 end
  47 rNoCompSa = rNoCompSa .* exp(-j * theta);
  48
  49
  50 %% Not entirely sure why, but after using FFT instead of time-domain
  51 %% convolution for simulating CD, we now need to do the same rotation
  52 %% for rSampled as well, but this time with a positive rotation.
  53 theta = angle(-sum(rSampled .^ M)) / M;
  54 if abs(theta + pi / M) / (pi / M) < 0.1
  55   theta = +pi / M;
  56 end
  57 rSampled = rSampled .* exp(-1j * theta);
  58
  59
  60
  61 %%rAdaptEq = adaptiveCMA(rSampled);
  62 %{
  63 %% Compare original signal and compensated signal
  64 figure(101);
  65 clf;
  66 tsym = t(sps*span/2+1:sps:(numSymbs+span/2)*sps);
  67 subplot(211);
  68 plot(t(1:length(x)), real(normalizeEnergy(x, numSymbs*sps, 1)), 'b');
  69 hold on
  70 plot(t(1:length(x)), real(normalizeEnergy(yCDComp(1:length(x)), numSymbs*sps, 1)), 'r');
  71 plot(tsym, real(rAdaptEq), 'x', 'Color', [0, 0.6, 0], 'LineWidth', 2);
  72 hold off;
  73 title('Real part');
  74 legend('original', 'dispersion compensated', 'CMA equalized samples');
  75 axis([t(6000*sps+1) t(6000*sps+150) -Inf +Inf]);
  76 subplot(212);
  77 plot(t(1:length(x)), imag(normalizeEnergy(x, numSymbs*sps, 1)), 'b');
  78 hold on;
  79 plot(t(1:length(x)), imag(normalizeEnergy(yCDComp(1:length(x)), numSymbs*sps, 1)), 'r');
  80 plot(tsym, imag(rAdaptEq), 'x', 'Color', [0, 0.6, 0], 'LineWidth', 2);
  81 hold off;
  82 title('Imag part');
  83 axis([t(6000*sps+1) t(6000*sps+150) -Inf +Inf]);
  84
  85 scatterplot(modData);
  86 formatFigure;
  87 %title('Constellation of original modulation', 'interpreter', 'latex');
  88 xlabel('In-Phase', 'interpreter', 'latex');
  89 %scatterplot(rSampled);
  90 %title('Constellation of matched filter output');
  91 scatterplot(rNoCompSa);
  92 title('Constellation of dispersed signal', 'interpreter', 'latex');
  93 scatterplot(rAdaptEq);
  94 title('Constellation of adaptive filter output');
  95 %}
  96 demodData = pskdemod(rSampled, M, pi / M, 'gray');
  97 %%demodAdapt = pskdemod(rAdaptEq, M, pi / M, 'gray');
  98
  99 [~, ber] = biterr(data, demodData)
 100 %[~, berNoComp] = biterr(data, pskdemod(rNoCompSa, M, pi/M, 'gray'))
 101 %[~, ber] = biterr(data, demodAdapt)