4 Rsym = 28e9; % symbol rate (sym/sec)
6 %% zs: array of distances z to be simulated
7 % Example: zs = 42; zs = 40:10:100; zs = [300, 500, 1000];
12 %% Tx RRC filter properties
14 span = 6; % filter span
15 sps = 16; % samples per symbol
21 t = (0 : 1 / fs : numSymbs / Rsym - 1 / fs).';
24 %% Chromatic dispersion
25 D = 17; % ps / (nm km)
30 linewidthLO = 1e6; % Hz
32 %% Kerr effect / SSFS parameters
33 gamma = 1.2; % watt^-1 / km
35 dz = 2; % Step size, km
37 %% Polarization state rotation parameters
38 rot_omega = 1e3; % rad/s
41 %% Launch power, per wavelength channel
43 power = 10 .^ (power_dBm / 10) * 1e-3; % watts
46 wavelength_channels = 3;
47 dw = 2 * pi * 50e9; % channel spacing (rad/s)
50 hc = 6.62607015e-34 * 299792458; % J m
51 Eperphoton = hc / (lambda * 1e-9); % J
54 %% Stores result to be plotted
55 ber = zeros(plotlen, 1);
57 fig = figure; hold on;
60 %% sps and Tsamp change at Tx/Rx, save these for later.
64 %% Generate random data for both polarizations
65 data_x = randi([0, M - 1], numSymbs, wavelength_channels, 'uint8');
66 data_y = randi([0, M - 1], numSymbs, wavelength_channels, 'uint8');
69 modData_x = deqpskmod(data_x);
70 modData_y = deqpskmod(data_y);
72 %% Construct waveforms for each channel separately
73 A_x_wdm = zeros(numSymbs * sps, wavelength_channels);
74 A_y_wdm = zeros(numSymbs * sps, wavelength_channels);
75 carriers = zeros(numSymbs * sps, wavelength_channels);
77 for w = 1 : wavelength_channels
78 %% Compute frequency offsets:
80 % Spectrum | | | | | | | | | |
82 % ____| |___| |___| |___| |___| |____ --> freq
84 % ang freq offset -2dw -dw 0 +dw +2dw
91 carriers(:, w) = exp(1j * ndw * t);
92 A_x_wdm(:, w) = txFilter(modData_x(:, w), rolloff, span, sps);
93 A_y_wdm(:, w) = txFilter(modData_y(:, w), rolloff, span, sps);
96 %% Sum the WDM waveforms with their frequency offsets
97 A_x = sum(A_x_wdm .* carriers, 2);
98 A_y = sum(A_y_wdm .* carriers, 2);
100 %% Clear variables no longer needed to reduce memory usage
101 clear modData_x modData_y A_x_wdm A_y_wdm;
103 %% Set launch power. Divide by 2 because half power for each polarization.
104 A_x = sqrt(power / 2) * A_x;
105 A_y = sqrt(power / 2) * A_y;
107 %% Rotate polarization states
108 A_x = A_x .* cos(rot_omega * t) + ...
109 A_y .* sin(rot_omega * t) * exp(-1j * rot_phi);
110 A_y = A_x .* -sin(rot_omega * t) * exp(1j * rot_phi) + ...
111 A_y .* cos(rot_omega * t);
113 %% Now loop through each z
120 %% Split-step Fourier
121 [A_x, A_y] = ssfs(A_x, A_y, D, lambda, z, dz, Tsamp, gamma, alpha);
124 A_x = phaseNoise(A_x, linewidthTx, linewidthLO, Tsamp);
125 A_y = phaseNoise(A_y, linewidthTx, linewidthLO, Tsamp);
127 %% Here, only receive the central channel 1.
128 % For channel n: A_x .* conj(carriers(:, n)); etc.
129 r_x = rxFilter(A_x, rolloff, span, sps);
130 r_y = rxFilter(A_y, rolloff, span, sps);
131 % Rx filter performs downsampling as well, keep track of this
133 Tsamp = Tsamp * spsOrig / sps;
136 photonpersym = mean(abs(r_x) .^ 2) / Rsym / Eperphoton;
138 r_x = awgn(r_x, snr, 'measured', 'linear');
139 r_y = awgn(r_y, snr, 'measured', 'linear');
141 %% -- Begin DSP channel equalization --
142 %% Chromatic dispersion compensation
143 r_x = CDCompensation(r_x, D, lambda, z, Tsamp);
144 r_y = CDCompensation(r_y, D, lambda, z, Tsamp);
149 [r_x, r_y] = pdm_adaptiveCMA(r_x, r_y);
151 %% Phase noise correction
152 r_x = phaseNoiseCorr(r_x, M, 0, 40).';
153 r_y = phaseNoiseCorr(r_y, M, 0, 40).';
155 %% Demodulate DE-QPSK
156 demod_x = deqpskdemod(r_x);
157 demod_y = deqpskdemod(r_y);
159 %% Calculate and store BER
160 [~, ber(i)] = biterr([data_x(:, 1); data_y(:, 1)], [demod_x; demod_y]);
162 q = 20 * log10(erfcinv(2*ber)*sqrt(2));