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

[4yp.git] / chromaticDispersion1Signal.m

diff --git a/chromaticDispersion1Signal.m b/chromaticDispersion1Signal.m
index db25965..7d9e980 100644 (file)
--- a/chromaticDispersion1Signal.m
+++ b/chromaticDispersion1Signal.m
@@ -1,63 +1,85 @@
 M = 4;
-numSymbs = 1000;
+numSymbs = 100000;
 
-%% https://www.mathworks.com/help/comm/examples/passband-modulation-with-adjacent-channel-interference.html
 Rsym = 2.5e10; % symbol rate (sym/sec)
 
 span = 6; % Tx/Rx filter span
 rolloff = 0.25; % Tx/Rx RRC rolloff
 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).';
 
 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 = 1000; % km
+z = 10; % km
 
 [xCD, xCDkstart] = chromaticDispersion(x, D, lambda, z, Tsamp);
 xCD = normalizeEnergy(xCD, numSymbs, 1);
 
+EbN0_db = 8;
+snr = EbN0_db + 10 * log10(log2(M)) - 10 * log10(sps);
+noiseEnergy = 10 ^ (-snr / 10);
 
+%%y = awgn(xCD, snr, 'measured');
 y = xCD;
 
-
 yCDComp = CDCompensation(y, D, lambda, z, Tsamp);
+%%yCDComp = y;
+
+r = rxFilter(yCDComp, rolloff, span, sps);
+rSampled = r(sps*span/2+1:sps:(numSymbs + span/2) * sps);
+
+%% if no CD comp, then rotate constellation. Use:
+%{
+theta = angle(-sum(rSampled .^ M)) / M;
+%% if theta approx +pi/M, wrap to -pi/M
+if abs(theta - pi / M) / (pi / M) < 0.1
+  theta = -pi / M;
+end
+rSampled = rSampled .* exp(-j * theta);
+%}
+
+rAdaptEq = adaptiveCMA(rSampled);
 
 %% Compare original signal and compensated signal
-figure(1);
+figure(101);
+clf;
+tsym = t(sps*span/2+1:sps:(numSymbs+span/2)*sps);
 subplot(211);
-plot(real(x(1:300)));
+plot(t(1:length(x)), real(normalizeEnergy(x, numSymbs*sps, 1)), 'b');
 hold on
-plot(real(yCDComp(1:300)));
-hold off
+plot(t(1:length(x)), real(normalizeEnergy(yCDComp(1:length(x)), numSymbs*sps, 1)), 'r');
+plot(tsym, real(rAdaptEq), 'xg');
+hold off;
 title('Real part');
-legend('original', 'dispersion compensated');
+legend('original', 'dispersion compensated', 'CMA equalized samples');
+axis([t(6000*sps+1) t(6000*sps+150) -Inf +Inf]);
 subplot(212);
-plot(imag(x(1:300)));
-hold on
-plot(imag(yCDComp(1:300)));
-hold off
+plot(t(1:length(x)), imag(normalizeEnergy(x, numSymbs*sps, 1)), 'b');
+hold on;
+plot(t(1:length(x)), imag(normalizeEnergy(yCDComp(1:length(x)), numSymbs*sps, 1)), 'r');
+plot(tsym, imag(rAdaptEq), 'xg');
+hold off;
 title('Imag part');
-
-
-r = rxFilter(yCDComp, rolloff, span, sps);
-r = normalizeEnergy(r, numSymbs, 1); % Add noise energy if needed
-
-rSampled = r(sps*span/2+1:sps:(numSymbs + span/2) * sps);
+axis([t(6000*sps+1) t(6000*sps+150) -Inf +Inf]);
 
 scatterplot(modData);
 title('Constellation of original modulation');
 scatterplot(rSampled);
-title('Constellation of sampled received waveform');
+title('Constellation of matched filter output');
+scatterplot(rAdaptEq);
+title('Constellation of adaptive filter output');
+
+demodData = pskdemod(rSampled, M, pi / M, 'gray');
+demodAdapt = pskdemod(rAdaptEq, M, pi / M, 'gray');
 
-demodData = pskdemod(rSampled, M, 0, 'gray');
+[~, ber] = biterr(data, demodData)
+[~, ber] = biterr(data, demodAdapt)