| 1 | numSymbs = 2^16; |
| 2 | M = 4; % QPSK |
| 3 | |
| 4 | Rsym = 28e9; % symbol rate (sym/sec) |
| 5 | |
| 6 | %% zs: array of distances z to be simulated |
| 7 | % Example: zs = 42; zs = 40:10:100; zs = [300, 500, 1000]; |
| 8 | zs = 42; |
| 9 | plotlen = length(zs); |
| 10 | |
| 11 | %% Tx RRC filter properties |
| 12 | rolloff = 0.25; |
| 13 | span = 6; % filter span |
| 14 | sps = 16; % samples per symbol |
| 15 | |
| 16 | %% Sampling frequency |
| 17 | fs = Rsym * sps; % Hz |
| 18 | Tsamp = 1 / fs; % s |
| 19 | % t: time vector, s |
| 20 | t = (0 : 1 / fs : numSymbs / Rsym - 1 / fs).'; |
| 21 | |
| 22 | |
| 23 | %% Chromatic dispersion |
| 24 | D = 17; % ps / (nm km) |
| 25 | lambda = 1550; % nm |
| 26 | |
| 27 | %% Laser phase noise |
| 28 | linewidthTx = 0; % Hz |
| 29 | linewidthLO = 1e6; % Hz |
| 30 | |
| 31 | %% Kerr effect / SSFS parameters |
| 32 | gamma = 1.2; % watt^-1 / km |
| 33 | alpha = 0.2; % dB/km |
| 34 | dz = 2; % Step size, km |
| 35 | |
| 36 | %% Polarization state rotation parameters |
| 37 | rot_omega = 1e3; % rad/s |
| 38 | rot_phi = 2; % rad |
| 39 | |
| 40 | %% Launch power, per wavelength channel |
| 41 | power_dBm = 0; |
| 42 | power = 10 .^ (power_dBm / 10) * 1e-3; % watts |
| 43 | |
| 44 | %% WDM properties |
| 45 | wavelength_channels = 3; |
| 46 | dw = 2 * pi * 50e9; % channel spacing (rad/s) |
| 47 | |
| 48 | %% Shot noise |
| 49 | hc = 6.62607015e-34 * 299792458; % J m |
| 50 | Eperphoton = hc / (lambda * 1e-9); % J |
| 51 | |
| 52 | |
| 53 | %% Stores result to be plotted |
| 54 | ber = zeros(plotlen, 1); |
| 55 | if plotlen > 1 |
| 56 | fig = figure; |
| 57 | end |
| 58 | |
| 59 | %% sps and Tsamp change at Tx/Rx, save these for later. |
| 60 | spsOrig = sps; |
| 61 | TsampOrig = Tsamp; |
| 62 | |
| 63 | %% Generate random data for both polarizations |
| 64 | data_x = randi([0, M - 1], numSymbs, wavelength_channels, 'uint8'); |
| 65 | data_y = randi([0, M - 1], numSymbs, wavelength_channels, 'uint8'); |
| 66 | |
| 67 | %% DE-QPSK modulation |
| 68 | modData_x = deqpskmod(data_x); |
| 69 | modData_y = deqpskmod(data_y); |
| 70 | |
| 71 | %% Construct waveforms for each channel separately |
| 72 | A_x_wdm = zeros(numSymbs * sps, wavelength_channels); |
| 73 | A_y_wdm = zeros(numSymbs * sps, wavelength_channels); |
| 74 | carriers = zeros(numSymbs * sps, wavelength_channels); |
| 75 | |
| 76 | for w = 1 : wavelength_channels |
| 77 | %% Compute frequency offsets: |
| 78 | % ___ ___ ___ ___ ___ |
| 79 | % Spectrum | | | | | | | | | | |
| 80 | % | | | | | | | | | | |
| 81 | % ____| |___| |___| |___| |___| |____ --> freq |
| 82 | % channel # 5 3 1 2 4 |
| 83 | % ang freq offset -2dw -dw 0 +dw +2dw |
| 84 | |
| 85 | if mod(w, 2) == 0 |
| 86 | ndw = w / 2 * dw; |
| 87 | else |
| 88 | ndw = (1-w) / 2 * dw; |
| 89 | end |
| 90 | carriers(:, w) = exp(1j * ndw * t); |
| 91 | A_x_wdm(:, w) = txFilter(modData_x(:, w), rolloff, span, sps); |
| 92 | A_y_wdm(:, w) = txFilter(modData_y(:, w), rolloff, span, sps); |
| 93 | end |
| 94 | |
| 95 | %% Sum the WDM waveforms with their frequency offsets |
| 96 | A_x = sum(A_x_wdm .* carriers, 2); |
| 97 | A_y = sum(A_y_wdm .* carriers, 2); |
| 98 | |
| 99 | %% Clear variables no longer needed to reduce memory usage |
| 100 | clear modData_x modData_y A_x_wdm A_y_wdm; |
| 101 | |
| 102 | %% Set launch power. Divide by 2 because half power for each polarization. |
| 103 | A_x = sqrt(power / 2) * A_x; |
| 104 | A_y = sqrt(power / 2) * A_y; |
| 105 | |
| 106 | %% Rotate polarization states |
| 107 | A_x = A_x .* cos(rot_omega * t) + ... |
| 108 | A_y .* sin(rot_omega * t) * exp(-1j * rot_phi); |
| 109 | A_y = A_x .* -sin(rot_omega * t) * exp(1j * rot_phi) + ... |
| 110 | A_y .* cos(rot_omega * t); |
| 111 | |
| 112 | %% Now loop through each z |
| 113 | for i = 1 : plotlen |
| 114 | z = zs(i); |
| 115 | |
| 116 | sps = spsOrig; |
| 117 | Tsamp = TsampOrig; |
| 118 | |
| 119 | %% Split-step Fourier |
| 120 | [Al_x, Al_y] = ssfs(A_x, A_y, D, lambda, z, dz, Tsamp, gamma, alpha); |
| 121 | |
| 122 | %% Phase noise |
| 123 | Al_x = phaseNoise(Al_x, linewidthTx, linewidthLO, Tsamp); |
| 124 | Al_y = phaseNoise(Al_y, linewidthTx, linewidthLO, Tsamp); |
| 125 | |
| 126 | %% Here, only receive the central channel 1. |
| 127 | % For channel n: Al_x .* conj(carriers(:, n)); etc. |
| 128 | r_x = rxFilter(Al_x, rolloff, span, sps); |
| 129 | r_y = rxFilter(Al_y, rolloff, span, sps); |
| 130 | % Rx filter performs downsampling as well, keep track of this |
| 131 | sps = 2; |
| 132 | Tsamp = Tsamp * spsOrig / sps; |
| 133 | |
| 134 | %% Rx shot noise |
| 135 | photonpersym = mean(abs(r_x) .^ 2) / Rsym / Eperphoton; |
| 136 | snr = photonpersym; |
| 137 | r_x = awgn(r_x, snr, 'measured', 'linear'); |
| 138 | r_y = awgn(r_y, snr, 'measured', 'linear'); |
| 139 | |
| 140 | %% -- Begin DSP channel equalization -- |
| 141 | %% Chromatic dispersion compensation |
| 142 | r_x = CDCompensation(r_x, D, lambda, z, Tsamp); |
| 143 | r_y = CDCompensation(r_y, D, lambda, z, Tsamp); |
| 144 | r_x = r_x(2:2:end); |
| 145 | r_y = r_y(2:2:end); |
| 146 | |
| 147 | %% Adaptive filter |
| 148 | [r_x, r_y] = pdm_adaptiveCMA(r_x, r_y); |
| 149 | |
| 150 | %% Phase noise correction |
| 151 | r_x = phaseNoiseCorr(r_x, M, 0, 40).'; |
| 152 | r_y = phaseNoiseCorr(r_y, M, 0, 40).'; |
| 153 | |
| 154 | %% Demodulate DE-QPSK |
| 155 | demod_x = deqpskdemod(r_x); |
| 156 | demod_y = deqpskdemod(r_y); |
| 157 | |
| 158 | %% Calculate and store BER |
| 159 | [~, ber(i)] = biterr([data_x(:, 1); data_y(:, 1)], [demod_x; demod_y]); |
| 160 | |
| 161 | q = 20 * log10(erfcinv(2*ber)*sqrt(2)); |
| 162 | if i > 1 |
| 163 | figure(fig); |
| 164 | plot(zs, q); |
| 165 | end |
| 166 | end |
| 167 | |
| 168 | ber |
| 169 | q |