X-Git-Url: https://adrianiainlam.tk/git/?a=blobdiff_plain;f=CD_AWGN.m;h=241d22a66839a16e93452f377593fd4aeedeaca7;hb=d9a5ae9d4e5369f11357d1f5825acc379fbe7e52;hp=be7de8cc9d3a1a5e868a92af5f878e64ce3c3633;hpb=5e9be3c421c4d52da9df842548f421751fa294d4;p=4yp.git

diff --git a/CD_AWGN.m b/CD_AWGN.m
index be7de8c..241d22a 100644
--- a/CD_AWGN.m
+++ b/CD_AWGN.m
@@ -1,18 +1,18 @@
-numSymbs = 10000;
+numSymbs = 2^16;
 M = 4;
 
 Rsym = 2.5e10; % symbol rate (sym/sec)
 
 rolloff = 0.25;
 span = 6; % filter span
-sps = 4; % samples per symbol
+sps = 8; % samples per symbol
 
 fs = Rsym * sps; % sampling freq (Hz)
 Tsamp = 1 / fs;
 
 t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).';
 
-EbN0_db = 0:0.2:14;
+EbN0_db = 0:0.5:14;
 EbN0 = 10 .^ (EbN0_db ./ 10);
 
 Es = 1;
@@ -36,26 +36,30 @@ x = txFilter(modData, rolloff, span, sps);
 %% Simulate chromatic dispersion
 D = 17; % ps / (nm km)
 lambda = 1550; % nm
-z = 10; % km
+z = 3000; % km
 
-xCD = chromaticDispersion(x, D, lambda, z, Tsamp);
-xCD = normalizeEnergy(xCD, numSymbs, 1);
+
+[xCD, xCDkstart] = chromaticDispersion(x, D, lambda, z, Tsamp);
+
+TsampOrig = Tsamp;
 
 for i = 1:plotlen
+  sps = 8;
+
   snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps);
-  noiseEnergy = 10 ^ (-snr / 10);
 
   y = awgn(xCD, snr, 'measured');
 
-  yCDComp = CDCompensation(y, D, lambda, z, Tsamp);
+  r = rxFilter(y, rolloff, span, sps);
 
-  r = rxFilter(yCDComp, rolloff, span, sps);
-  rNoComp = rxFilter(y, rolloff, span, sps);
-  %% normalize energy
-  %r = normalizeEnergy(r, numSymbs, 1 + noiseEnergy);
+  sps = 2;
+  Tsamp = TsampOrig * 4;
 
-  rSampled = r(sps*span/2+1:sps:(numSymbs + span/2) * sps);
-  rNoCompSampled = rNoComp(sps*span/2+1:sps:(numSymbs+span/2)*sps);
+  [rCDComp, CDCompkstart] = CDCompensation(r, D, lambda, z, Tsamp);
+  rCDComp = normalizeEnergy(rCDComp, numSymbs*sps, 1);
+
+  rSampled = rCDComp(2:2:end);
+  rNoCompSampled = r(2:2:end);
 
   %% rotate rNoCompSampled to match original data
   theta = angle(-sum(rNoCompSampled .^ M)) / M;
@@ -63,9 +67,20 @@ for i = 1:plotlen
   if abs(theta - pi / M) / (pi / M) < 0.1
     theta = -pi / M;
   end
-  theta
   rNoCompSampled = rNoCompSampled .* exp(-j * theta);
 
+
+  %% Not entirely sure why, but after using FFT instead of time-domain
+  %% convolution for simulating CD, we now need to do the same rotation
+  %% for rSampled as well, but this time with a positive rotation.
+  theta = angle(-sum(rSampled .^ M)) / M;
+  if abs(theta + pi / M) / (pi / M) < 0.1
+    theta = +pi / M;
+  end
+  rSampled = rSampled .* exp(-1j * theta);
+
+
+
   %% adaptive filter
   adaptFilterOut = adaptiveCMA(rSampled);
 
@@ -79,25 +94,33 @@ for i = 1:plotlen
   [~, berAdapt(i)] = biterr(data, demodAdapt);
   %%[~, berMatlabAdapt(i)] = biterr(data, demodMatlabAdapt);
 
-
-  if EbN0_db(i) == 12
+%{
+  if EbN0_db(i) == 14
     figure(1);
-    scatterplot(rSampled);
-    title('Constellation after CD compensation');
+    scatterplot(normalizeEnergy(rSampled, numSymbs, 1));
+    formatFigure;
+    title('Constellation after CD comp.', 'interpreter', 'latex');
+    xlabel('In-Phase', 'interpreter', 'latex');
+    ylabel('Quadrature', 'interpreter', 'latex');
+    set(gca, 'FontSize', 18);
     %%scatterplot(modData);
     %%title('Original constellation');
-    scatterplot(rNoCompSampled);
-    title('Constellation without CD compensation');
-    scatterplot(adaptFilterOut);
-    title('Constellation with CD compensation and adaptive filter');
+    scatterplot(normalizeEnergy(rNoCompSampled, numSymbs, 1));
+    formatFigure;
+    title('Constellation without CD comp.', 'interpreter', 'latex');
+    xlabel('In-Phase', 'interpreter', 'latex');
+    ylabel('Quadrature', 'interpreter', 'latex');
+    set(gca, 'FontSize', 18);
+    %scatterplot(adaptFilterOut);
+    %title('Constellation with CD compensation and adaptive filter');
     %scatterplot(matlabEq);
     %title('Matlab equalizer');
     ber(i)
     %berNoComp(i)
-    berAdapt(i)
-    berMatlabAdapt(i)
+    %berAdapt(i)
+    %berMatlabAdapt(i)
   end
-
+%}
 end
 
 figure(1);
@@ -106,15 +129,22 @@ clf;
 %% Plot simulated results
 semilogy(EbN0_db, ber, 'r', 'LineWidth', 2);
 hold on;
-%%semilogy(EbN0_db, berNoComp, 'g', 'LineWidth', 2);
-semilogy(EbN0_db, berAdapt, 'm', 'LineWidth', 1.4);
+semilogy(EbN0_db, berNoComp, 'm', 'LineWidth', 2);
+semilogy(EbN0_db, berAdapt, 'Color', [0, 0.6, 0], 'LineWidth', 2);
 %%%semilogy(EbN0_db, berMatlabAdapt, 'c', 'LineWidth', 1.4);
 
 theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1);
-legend({'CD + AWGN + CD comp.', 'CD + AWGN + CD comp.~+ CMA', ...
-        'Theoretical AWGN'}, 'Location', 'southwest');
-
-title(strcat(num2str(M), '-PSK with chromatic dispersion and compensation'));
+%%legend({'CD + AWGN + CD comp.', 'CD + AWGN + CD comp.~+ CMA', ...
+%%        'Theoretical AWGN'}, 'Location', 'southwest');
+%%legend({'CD + AWGN + CD comp.', 'CD + AWGN', 'Theoretical AWGN'}, ...
+%%       'Location', 'southwest');
+legend({'CD + AWGN + CD comp.', 'CD + AWGN', ...
+        'CD + AWGN + CD comp.~+ CMA', 'Theoretical AWGN'}, 'Location', ...
+       'Southwest');
+
+%%title(strcat(num2str(M), '-PSK with chromatic dispersion and compensation'));
+title({'QPSK with chromatic dispersion and compensation', ...
+       strcat(['$D = 17$ ps/(nm km), $z = ', num2str(z), '$ km'])});
 grid on;
 xlabel('$E_b/N_0$ (dB)');
 ylabel('BER');