Added adaptive CMA equalizer, "fixed" CD and phase noise

[4yp.git] / CD_AWGN.m

diff --git a/CD_AWGN.m b/CD_AWGN.m
index 1ea6ad8..be7de8c 100644 (file)
--- a/CD_AWGN.m
+++ b/CD_AWGN.m
@@ -10,7 +10,7 @@ sps = 4; % samples per symbol
 fs = Rsym * sps; % sampling freq (Hz)
 Tsamp = 1 / fs;
 
-t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs)';
+t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).';
 
 EbN0_db = 0:0.2:14;
 EbN0 = 10 .^ (EbN0_db ./ 10);
@@ -25,20 +25,22 @@ EsN0_db = 10 .* log10(EsN0);
 plotlen = length(EbN0);
 
 ber = zeros(1, plotlen);
+berNoComp = zeros(1, plotlen);
+berAdapt = zeros(1, plotlen);
+berMatlabAdapt = zeros(1, plotlen);
 
 data = randi([0 M - 1], numSymbs, 1);
-modData = pskmod(data, M, 0, 'gray');
+modData = pskmod(data, M, pi / M, 'gray');
 x = txFilter(modData, rolloff, span, sps);
 
 %% Simulate chromatic dispersion
-D = 20; % ps / (nm km)
+D = 17; % ps / (nm km)
 lambda = 1550; % nm
 z = 10; % km
 
 xCD = chromaticDispersion(x, D, lambda, z, Tsamp);
 xCD = normalizeEnergy(xCD, numSymbs, 1);
 
-
 for i = 1:plotlen
   snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps);
   noiseEnergy = 10 ^ (-snr / 10);
@@ -48,13 +50,54 @@ for i = 1:plotlen
   yCDComp = CDCompensation(y, D, lambda, z, Tsamp);
 
   r = rxFilter(yCDComp, rolloff, span, sps);
+  rNoComp = rxFilter(y, rolloff, span, sps);
   %% normalize energy
   %r = normalizeEnergy(r, numSymbs, 1 + noiseEnergy);
 
   rSampled = r(sps*span/2+1:sps:(numSymbs + span/2) * sps);
-  demodData = pskdemod(rSampled, M, 0, 'gray');
+  rNoCompSampled = rNoComp(sps*span/2+1:sps:(numSymbs+span/2)*sps);
+
+  %% rotate rNoCompSampled to match original data
+  theta = angle(-sum(rNoCompSampled .^ M)) / M;
+  %% if theta approx +pi/M, wrap to -pi/M
+  if abs(theta - pi / M) / (pi / M) < 0.1
+    theta = -pi / M;
+  end
+  theta
+  rNoCompSampled = rNoCompSampled .* exp(-j * theta);
+
+  %% adaptive filter
+  adaptFilterOut = adaptiveCMA(rSampled);
+
+  demodData = pskdemod(rSampled, M, pi / M, 'gray');
+  demodNoComp = pskdemod(rNoCompSampled, M, pi / M, 'gray');
+  demodAdapt = pskdemod(adaptFilterOut, M, pi / M, 'gray');
+  %%demodMatlabAdapt = pskdemod(matlabEq, M, pi / M, 'gray');
 
   [bitErrors, ber(i)] = biterr(data, demodData);
+  [bitErrors, berNoComp(i)] = biterr(data, demodNoComp);
+  [~, berAdapt(i)] = biterr(data, demodAdapt);
+  %%[~, berMatlabAdapt(i)] = biterr(data, demodMatlabAdapt);
+
+
+  if EbN0_db(i) == 12
+    figure(1);
+    scatterplot(rSampled);
+    title('Constellation after CD compensation');
+    %%scatterplot(modData);
+    %%title('Original constellation');
+    scatterplot(rNoCompSampled);
+    title('Constellation without CD compensation');
+    scatterplot(adaptFilterOut);
+    title('Constellation with CD compensation and adaptive filter');
+    %scatterplot(matlabEq);
+    %title('Matlab equalizer');
+    ber(i)
+    %berNoComp(i)
+    berAdapt(i)
+    berMatlabAdapt(i)
+  end
+
 end
 
 figure(1);
@@ -63,9 +106,13 @@ 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, berMatlabAdapt, 'c', 'LineWidth', 1.4);
 
 theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1);
-legend({'CD + AWGN + CD compensation', 'AWGN only'}, 'Location', 'southwest');
+legend({'CD + AWGN + CD comp.', 'CD + AWGN + CD comp.~+ CMA', ...
+        'Theoretical AWGN'}, 'Location', 'southwest');
 
 title(strcat(num2str(M), '-PSK with chromatic dispersion and compensation'));
 grid on;