Added presentation; DE-QPSK; CD with FFT; split-step Fourier

[4yp.git] / chromaticDispersion1Signal.m

diff --git a/chromaticDispersion1Signal.m b/chromaticDispersion1Signal.m
index 7d9e980..dd169ac 100644 (file)
--- a/chromaticDispersion1Signal.m
+++ b/chromaticDispersion1Signal.m
@@ -1,11 +1,11 @@
 M = 4;
-numSymbs = 100000;
+numSymbs = 5e5;
 
 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
+sps = 2; % samples per symbol
 
 fs = Rsym * sps; % sampling freq (Hz)
 Tsamp = 1 / fs;
@@ -16,39 +16,50 @@ data = randi([0 M - 1], numSymbs, 1);
 modData = pskmod(data, M, pi / M, 'gray');
 x = txFilter(modData, rolloff, span, sps);
 
+x = normalizeEnergy(x, numSymbs*sps, 1);
+
 %% Simulate chromatic dispersion
 D = 17; % ps / (nm km)
 lambda = 1550; % nm
-z = 10; % km
+z = 5000 % km
 
-[xCD, xCDkstart] = chromaticDispersion(x, D, lambda, z, Tsamp);
-xCD = normalizeEnergy(xCD, numSymbs, 1);
+[xCD, xCDkstart] = chromaticDispersion_FFT(x, D, lambda, z, Tsamp);
 
 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(y, rolloff, span, sps);
+[rCDComp, CDCompkstart] = CDCompensation(r, D, lambda, z, Tsamp);
+rCDComp = normalizeEnergy(rCDComp, numSymbs*sps, 1);
 
-r = rxFilter(yCDComp, rolloff, span, sps);
-rSampled = r(sps*span/2+1:sps:(numSymbs + span/2) * sps);
+rSampled = rCDComp(sps*span/2+1:sps:(numSymbs+span/2)*sps);
+rNoCompSa = r(sps*span/2+1:sps:(numSymbs+span/2)*sps);
 
 %% if no CD comp, then rotate constellation. Use:
-%{
-theta = angle(-sum(rSampled .^ M)) / M;
+theta = angle(-sum(rNoCompSa .^ 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);
-%}
+rNoCompSa = rNoCompSa .* 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);
+
 
-rAdaptEq = adaptiveCMA(rSampled);
 
+%%rAdaptEq = adaptiveCMA(rSampled);
+%{
 %% Compare original signal and compensated signal
 figure(101);
 clf;
@@ -57,7 +68,7 @@ subplot(211);
 plot(t(1:length(x)), real(normalizeEnergy(x, numSymbs*sps, 1)), 'b');
 hold on
 plot(t(1:length(x)), real(normalizeEnergy(yCDComp(1:length(x)), numSymbs*sps, 1)), 'r');
-plot(tsym, real(rAdaptEq), 'xg');
+plot(tsym, real(rAdaptEq), 'x', 'Color', [0, 0.6, 0], 'LineWidth', 2);
 hold off;
 title('Real part');
 legend('original', 'dispersion compensated', 'CMA equalized samples');
@@ -66,20 +77,25 @@ subplot(212);
 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');
+plot(tsym, imag(rAdaptEq), 'x', 'Color', [0, 0.6, 0], 'LineWidth', 2);
 hold off;
 title('Imag part');
 axis([t(6000*sps+1) t(6000*sps+150) -Inf +Inf]);
 
 scatterplot(modData);
-title('Constellation of original modulation');
-scatterplot(rSampled);
-title('Constellation of matched filter output');
+formatFigure;
+%title('Constellation of original modulation', 'interpreter', 'latex');
+xlabel('In-Phase', 'interpreter', 'latex');
+%scatterplot(rSampled);
+%title('Constellation of matched filter output');
+scatterplot(rNoCompSa);
+title('Constellation of dispersed signal', 'interpreter', 'latex');
 scatterplot(rAdaptEq);
 title('Constellation of adaptive filter output');
-
+%}
 demodData = pskdemod(rSampled, M, pi / M, 'gray');
-demodAdapt = pskdemod(rAdaptEq, M, pi / M, 'gray');
+%%demodAdapt = pskdemod(rAdaptEq, M, pi / M, 'gray');
 
 [~, ber] = biterr(data, demodData)
-[~, ber] = biterr(data, demodAdapt)
+%[~, berNoComp] = biterr(data, pskdemod(rNoCompSa, M, pi/M, 'gray'))
+%[~, ber] = biterr(data, demodAdapt)