From 5117fc5819545753d1c8ebbef03bf7deea0f3ae9 Mon Sep 17 00:00:00 2001 From: Adrian Iain Lam Date: Thu, 23 May 2019 03:16:04 +0100 Subject: [PATCH] Final simulation model --- CDCompensation.m | 37 ++++++++++++ deqpskdemod.m | 19 ++++++ deqpskmod.m | 10 ++++ main.m | 170 ++++++++++++++++++++++++++++++++++++++++++++++++++++++ pdm_adaptiveCMA.m | 76 ++++++++++++++++++++++++ phaseNoise.m | 18 ++++++ phaseNoiseCorr.m | 33 +++++++++++ rxFilter.m | 32 ++++++++++ ssfs.m | 70 ++++++++++++++++++++++ txFilter.m | 28 +++++++++ 10 files changed, 493 insertions(+) create mode 100644 CDCompensation.m create mode 100644 deqpskdemod.m create mode 100644 deqpskmod.m create mode 100644 main.m create mode 100644 pdm_adaptiveCMA.m create mode 100644 phaseNoise.m create mode 100644 phaseNoiseCorr.m create mode 100644 rxFilter.m create mode 100644 ssfs.m create mode 100644 txFilter.m diff --git a/CDCompensation.m b/CDCompensation.m new file mode 100644 index 0000000..b220bcf --- /dev/null +++ b/CDCompensation.m @@ -0,0 +1,37 @@ +function x = CDCompensation(xCD, D, lambda, z, Tsamp) + %% Chromatic dispersion compensation. + %% Params: + %% - xCD: received waveform with CD + %% - D: dispersion coefficient (ps / (nm km)) + %% - lambda: wavelength (nm) + %% - z: length of fibre (km) + %% - Tsamp: sampling time (s) + %% Output: + %% - x: xCD after performing CD compensation + + %% Convert everything to SI base units + c = 299792458; % m/s + D = D * 1e-6; % s/m^2 + lambda = lambda * 1e-9; % m + z = z * 1e3; % m + + %% Discrete FIR filter: + %% N: filter length; k: filter index; h: filter coefficients. + N = 2 * floor(abs(D) * lambda^2 * z / (2 * c * Tsamp^2)) + 1; + k = -floor(N / 2) : floor(N / 2); + h = exp(-1j * pi * c * Tsamp^2 * k .^ 2 / (D * lambda^2 * z)); + + %% Perform filtering in frequency domain + len_fft = max(length(xCD), length(h)); + H = fft(h, len_fft); + XCD = fft(xCD, len_fft); + x = ifft(H.' .* XCD); + + %% Re-order due to circular convolution + l = (N - 1) / 2; + if l > 0 + x = [x(l:end); x(1:l-1)]; + else + x = [x(end); x(1:end-1)]; + end +end diff --git a/deqpskdemod.m b/deqpskdemod.m new file mode 100644 index 0000000..fbcef9c --- /dev/null +++ b/deqpskdemod.m @@ -0,0 +1,19 @@ +function demodData = deqpskdemod(modData) + %% Input: Rx samples from a DE-QPSK constellation. + %% Output: Data as integers between 0 to 3 inclusive. + %% Parameter formats are the same as MATLAB's "pskdemod" function. + + %% MATLAB has no easy way to perform just the decision stage, + % so here we demodulate and then remodulate to get the closest symbols. + demod1 = pskdemod(modData, 4, 0, 'gray'); + remod = pskmod(demod1, 4, 0, 'gray'); + + clear demod1; % save some memory + + [~, wavelength_channels] = size(modData); + %% Instead of looping through the symbols, it will be faster + % to vectorize the operation. So we create an array that is + % delayed by 1 sample. + delayed = [ones(1, wavelength_channels); remod(1:end-1, :)]; + demodData = uint8(pskdemod(remod .* conj(delayed), 4, 0, 'gray')); +end diff --git a/deqpskmod.m b/deqpskmod.m new file mode 100644 index 0000000..13a8459 --- /dev/null +++ b/deqpskmod.m @@ -0,0 +1,10 @@ +function modData = deqpskmod(data) + %% Input: Unmodulated data as integers between 0 to 3 inclusive. + %% Output: Differentially encoded QPSK symbols. + %% Parameter formats are the same as MATLAB's "pskmod" function. + modData = pskmod(data, 4, 0, 'gray'); + numSymbs = size(data, 1); + for i = 2 : numSymbs + modData(i, :) = modData(i, :) .* modData(i - 1, :); + end +end diff --git a/main.m b/main.m new file mode 100644 index 0000000..38dcccc --- /dev/null +++ b/main.m @@ -0,0 +1,170 @@ +numSymbs = 2^16; +M = 4; % QPSK + +Rsym = 28e9; % symbol rate (sym/sec) + +%% zs: array of distances z to be simulated +% Example: zs = 42; zs = 40:10:100; zs = [300, 500, 1000]; +zs = 42; +plotlen = length(zs); + + +%% Tx RRC filter properties +rolloff = 0.25; +span = 6; % filter span +sps = 16; % samples per symbol + +%% Sampling frequency +fs = Rsym * sps; % Hz +Tsamp = 1 / fs; % s +% t: time vector, s +t = (0 : 1 / fs : numSymbs / Rsym - 1 / fs).'; + + +%% Chromatic dispersion +D = 17; % ps / (nm km) +lambda = 1550; % nm + +%% Laser phase noise +linewidthTx = 0; % Hz +linewidthLO = 1e6; % Hz + +%% Kerr effect / SSFS parameters +gamma = 1.2; % watt^-1 / km +alpha = 0.2; % dB/km +dz = 2; % Step size, km + +%% Polarization state rotation parameters +rot_omega = 1e3; % rad/s +rot_phi = 2; % rad + +%% Launch power, per wavelength channel +power_dBm = 0; +power = 10 .^ (power_dBm / 10) * 1e-3; % watts + +%% WDM properties +wavelength_channels = 3; +dw = 2 * pi * 50e9; % channel spacing (rad/s) + +%% Shot noise +hc = 6.62607015e-34 * 299792458; % J m +Eperphoton = hc / (lambda * 1e-9); % J + + +%% Stores result to be plotted +ber = zeros(plotlen, 1); +if plotlen > 1 + fig = figure; hold on; +end + +%% sps and Tsamp change at Tx/Rx, save these for later. +spsOrig = sps; +TsampOrig = Tsamp; + +%% Generate random data for both polarizations +data_x = randi([0, M - 1], numSymbs, wavelength_channels, 'uint8'); +data_y = randi([0, M - 1], numSymbs, wavelength_channels, 'uint8'); + +%% DE-QPSK modulation +modData_x = deqpskmod(data_x); +modData_y = deqpskmod(data_y); + +%% Construct waveforms for each channel separately +A_x_wdm = zeros(numSymbs * sps, wavelength_channels); +A_y_wdm = zeros(numSymbs * sps, wavelength_channels); +carriers = zeros(numSymbs * sps, wavelength_channels); + +for w = 1 : wavelength_channels + %% Compute frequency offsets: + % ___ ___ ___ ___ ___ + % Spectrum | | | | | | | | | | + % | | | | | | | | | | + % ____| |___| |___| |___| |___| |____ --> freq + % channel # 5 3 1 2 4 + % ang freq offset -2dw -dw 0 +dw +2dw + + if mod(w, 2) == 0 + ndw = w / 2 * dw; + else + ndw = (1-w) / 2 * dw; + end + carriers(:, w) = exp(1j * ndw * t); + A_x_wdm(:, w) = txFilter(modData_x(:, w), rolloff, span, sps); + A_y_wdm(:, w) = txFilter(modData_y(:, w), rolloff, span, sps); +end + +%% Sum the WDM waveforms with their frequency offsets +A_x = sum(A_x_wdm .* carriers, 2); +A_y = sum(A_y_wdm .* carriers, 2); + +%% Clear variables no longer needed to reduce memory usage +clear modData_x modData_y A_x_wdm A_y_wdm; + +%% Set launch power. Divide by 2 because half power for each polarization. +A_x = sqrt(power / 2) * A_x; +A_y = sqrt(power / 2) * A_y; + +%% Rotate polarization states +A_x = A_x .* cos(rot_omega * t) + ... + A_y .* sin(rot_omega * t) * exp(-1j * rot_phi); +A_y = A_x .* -sin(rot_omega * t) * exp(1j * rot_phi) + ... + A_y .* cos(rot_omega * t); + +%% Now loop through each z +for i = 1 : plotlen + z = zs(i); + + sps = spsOrig; + Tsamp = TsampOrig; + + %% Split-step Fourier + [A_x, A_y] = ssfs(A_x, A_y, D, lambda, z, dz, Tsamp, gamma, alpha); + + %% Phase noise + A_x = phaseNoise(A_x, linewidthTx, linewidthLO, Tsamp); + A_y = phaseNoise(A_y, linewidthTx, linewidthLO, Tsamp); + + %% Here, only receive the central channel 1. + % For channel n: A_x .* conj(carriers(:, n)); etc. + r_x = rxFilter(A_x, rolloff, span, sps); + r_y = rxFilter(A_y, rolloff, span, sps); + % Rx filter performs downsampling as well, keep track of this + sps = 2; + Tsamp = Tsamp * spsOrig / sps; + + %% Rx shot noise + photonpersym = mean(abs(r_x) .^ 2) / Rsym / Eperphoton; + snr = photonpersym; + r_x = awgn(r_x, snr, 'measured', 'linear'); + r_y = awgn(r_y, snr, 'measured', 'linear'); + + %% -- Begin DSP channel equalization -- + %% Chromatic dispersion compensation + r_x = CDCompensation(r_x, D, lambda, z, Tsamp); + r_y = CDCompensation(r_y, D, lambda, z, Tsamp); + r_x = r_x(2:2:end); + r_y = r_y(2:2:end); + + %% Adaptive filter + [r_x, r_y] = pdm_adaptiveCMA(r_x, r_y); + + %% Phase noise correction + r_x = phaseNoiseCorr(r_x, M, 0, 40).'; + r_y = phaseNoiseCorr(r_y, M, 0, 40).'; + + %% Demodulate DE-QPSK + demod_x = deqpskdemod(r_x); + demod_y = deqpskdemod(r_y); + + %% Calculate and store BER + [~, ber(i)] = biterr([data_x(:, 1); data_y(:, 1)], [demod_x; demod_y]); + + q = 20 * log10(erfcinv(2*ber)*sqrt(2)); + if plotlen > 1 + figure(fig); + plot(zs, q); + end +end + +ber +q diff --git a/pdm_adaptiveCMA.m b/pdm_adaptiveCMA.m new file mode 100644 index 0000000..6dca946 --- /dev/null +++ b/pdm_adaptiveCMA.m @@ -0,0 +1,76 @@ +function [x, y] = pdm_adaptiveCMA(rx, ry) + %% Perform adaptive equalization using CMA. + %% Input: rx, ry: Both polarizations of received signal + %% Output: x, y: Equalizaed signal + + taps = 19; % Number of taps. Should be odd. + mu = 1e-3; % Convergence parameter for gradient descent. + + hxx = zeros(taps, 1); + hxy = zeros(taps, 1); + hyx = zeros(taps, 1); + hyy = zeros(taps, 1); + %% hxx: real indices -K, ..., 0, ..., K. K = floor(taps/2) + %% MATLAB indices 1 1+K taps + + %% Initialize hxx, hxx[0] = 1, hxx[k] = hxx[-k] = 0 + hxx(ceil(taps/2)) = 1; + hxy(ceil(taps/2)) = 1; + hyx(ceil(taps/2)) = 1; + hyy(ceil(taps/2)) = 1; + + numSymbs = length(rx); + + %% Normalize to unit power. + rx = rx / sqrt(mean(abs(rx) .^ 2)); + ry = ry / sqrt(mean(abs(ry) .^ 2)); + + x = zeros(numSymbs, 1); + y = zeros(numSymbs, 1); + + %% Run CMA twice so that the first symbols were also equalized + for loops = 1:2 + %% Loop through each symbol + for it = 1:numSymbs + %% Construct block of length equal to filter length (taps) + if it <= (taps - 1) / 2; + %% If near the start, prepend zeros + xp = [zeros((taps - 1) / 2 - it + 1, 1); rx(1:it + (taps - 1) / 2)]; + yp = [zeros((taps - 1) / 2 - it + 1, 1); ry(1:it + (taps - 1) / 2)]; + elseif it + (taps - 1) / 2 > numSymbs + %% If near the end, append zeros + xp = [rx(it - (taps - 1) / 2 : end); zeros(it + (taps - 1) / 2 - numSymbs, 1)]; + yp = [ry(it - (taps - 1) / 2 : end); zeros(it + (taps - 1) / 2 - numSymbs, 1)]; + else + %% Just slice the signal + xp = rx(it - (taps - 1) / 2 : it + (taps - 1) / 2); + yp = ry(it - (taps - 1) / 2 : it + (taps - 1) / 2); + end + + %% Filtering + xout = sum(hxx .* xp) + sum(hxy .* yp); + yout = sum(hyx .* xp) + sum(hyy .* yp); + x(it) = xout; + y(it) = yout; + + %% Caculate error + ex = 1 - abs(xout) ^ 2; + ey = 1 - abs(yout) ^ 2; + + %% Update filter by gradient descent + hxx = hxx + mu * ex * xout * conj(xp); + hxy = hxy + mu * ex * xout * conj(yp); + hyx = hyx + mu * ey * yout * conj(xp); + hyy = hyy + mu * ey * yout * conj(yp); + + %% If both filters converge to the same polarization, + % re-initialize the filters. + if sum(abs(hxx - hyx)) < 0.01 && sum(abs(hxy - hyy)) < 0.01 + hxx = 0.5 * (hxx + flipud(conj(hyy))); + hyy = conj(flipud(hxx)); + hxy = 0.5 * (hxy - conj(flipud(hyx))); + hyx = -conj(flipud(hxy)); + end + end + end +end diff --git a/phaseNoise.m b/phaseNoise.m new file mode 100644 index 0000000..4363c22 --- /dev/null +++ b/phaseNoise.m @@ -0,0 +1,18 @@ +function [xPN, phasenoise] = phaseNoise(x, linewidthTx, linewidthLO, Tsamp) + %% Simulates laser phase noise. + %% Inputs: + %% - x: input waveform + %% - linewidthTx: Tx laser linewidth (Hz) + %% - linewidthLO: Rx LO laser linewidth (Hz) + %% - Tsamp: Sampling period (s) + %% Outputs: + %% - xPN: output waveform + %% - phasenoise: the actual phase noise added (rad) + dphiTx = sqrt(2 * pi * linewidthTx * Tsamp) * randn(length(x), 1); + dphiLO = sqrt(2 * pi * linewidthLO * Tsamp) * randn(length(x), 1); + phiTx = cumsum(dphiTx); + phiLO = cumsum(dphiLO); + + phasenoise = phiTx - phiLO; + xPN = x .* exp(-1j * phasenoise); +end diff --git a/phaseNoiseCorr.m b/phaseNoiseCorr.m new file mode 100644 index 0000000..eead282 --- /dev/null +++ b/phaseNoiseCorr.m @@ -0,0 +1,33 @@ +function [rc, phiests] = phaseNoiseCorr(r, M, phoffset, blocksize) + %% Phase noise correction. + %% Inputs: + %% - r: Received symbols + %% - M: Order of constellation, M=4 for QPSK + %% - phoffset: Phase of the 0th symbol in the constellation + %% - blocksize: Viterbi-Viterbi algorithm block size + %% Outputs: + %% - rc: Symbols with corrected phase + %% - phiests: Estimates of the phase noise + + n = length(r); + phiests = zeros(1, n); + rc = zeros(1, n); + + for l = 1 : blocksize : n + block = r(l : min(l + blocksize - 1, n)); + + sum_M = sum(block .^ M); + phi_est = angle(sum_M .* exp(1j * M * phoffset)) / M; + + if l > 1 + %% Phase unwrapping + phi_prev = phiests(l - 1); + m = floor(0.5 + (phi_prev - phi_est) * M / (2 * pi)); + phi_est = phi_est + m * 2 * pi / M; + end + + block = block .* exp(1j * -phi_est); + rc(l : min(l + blocksize - 1, n)) = block; + phiests(l : min(l + blocksize - 1, n)) = phi_est; + end +end diff --git a/rxFilter.m b/rxFilter.m new file mode 100644 index 0000000..2868fa4 --- /dev/null +++ b/rxFilter.m @@ -0,0 +1,32 @@ +function r = rxFilter(y, rolloff, span, sps) + %% Receiver matched (root raised cosine) filter, + %% downsampling to 2 samples/sym. + %% Inputs: + %% - y: received waveform + %% - rolloff: rolloff factor in root raised cosine filter + %% - span: filter span (number of symbols) + %% - sps: Input samples per symbol + %% Output: + %% - r: filtered signal + + %% Construct filter object + rxfilter = comm.RaisedCosineReceiveFilter... + ('Shape', 'Square root', ... + 'RolloffFactor', rolloff, ... + 'FilterSpanInSymbols', span, ... + 'InputSamplesPerSymbol', sps, ... + 'Gain', 1 / sqrt(sps)); + + %% Perform filtering in frequency domain + coef = coeffs(rxfilter); + filter_fft = fft(coef.Numerator, length(y)); + y_fft = fft(y); + rs = ifft(y_fft .* filter_fft.'); + + %% Re-order signal due to circular convolution + l = (length(coef.Numerator) - 1) / 2; + rr = [rs(l:end); rs(1:l-1)]; + + %% Downsample + r = downsample(rr, sps/2, 2); +end diff --git a/ssfs.m b/ssfs.m new file mode 100644 index 0000000..2d3abef --- /dev/null +++ b/ssfs.m @@ -0,0 +1,70 @@ +function [x, y] = ssfs(xin, yin, D, lambda, z, dz, Tsamp, gamma, alpha) + %% Split-step Fourier solver, simulates chromatic dispersion, + %% attenuation, Kerr effect, and power splitting / amplification. + %% Params: + %% - xin, yin: input waveform (x and y polarizations) + %% - D: dispersion coefficient (ps / (nm km)) + %% - lambda: wavelength (nm) + %% - z: length of fibre (km) + %% - dz: step size (km) + %% - Tsamp: sampling time (s) + %% - gamma: Non-linear coefficient (W^-1 / km) + %% - alpha: attenuation (dB / km) + %% Output: + %% - x, y: output waveform (both polarizations) + + %% Convert everything to SI base units + c = 299792458; % m/s + D = D * 1e-6; % s/m^2 + lambda = lambda * 1e-9; % m + z = z * 1e3; % m + gamma = gamma * 1e-3; % watt^-1 / m + dz = dz * 1e3; % m + alpha = alpha / 10 * log(10); % Np/km + alpha = alpha * 1e-3; % Np/m + + stepnum = z / dz; + + %% Frequency response of CD + n = length(xin); + fs = 1 / Tsamp; + omega = (2*pi * fs / n * [(0 : floor((n-1)/2)), (-ceil((n-1)/2) : -1)]).'; + dispDFT = exp(-1j * omega.^2 * D * lambda^2 * dz / (4 * pi * c)); + + %% Convenient variables to reduce typing + hhz = 1j * (8/9) * gamma * dz; % Kerr phase shift per power + attn = -alpha * dz / 2; % attenuation + + %% Initial Kerr half step + P = abs(xin) .^ 2 + abs(yin) .^ 2; + x = xin .* exp(P .* hhz / 2 + attn / 2); + y = yin .* exp(P .* hhz / 2 + attn / 2); + + for i = 1 : stepnum + %% CD in frequency domain + xDFT = fft(x); + yDFT = fft(y); + x = ifft(xDFT .* dispDFT); + y = ifft(yDFT .* dispDFT); + + %% Kerr effect in time domain + P = abs(x) .^ 2 + abs(y) .^ 2; + x = x .* exp(P .* hhz + attn); + y = y .* exp(P .* hhz + attn); + + %% Power split after 40km + if i * dz == 40e3 + splitnum = 1024; % energy factor, amplitude factor is sqrt of this + x = x ./ sqrt(splitnum); + y = y ./ sqrt(splitnum); + %% Splitter loss - 1:4 coupler * 5 levels, 0.3 dB per level + %% so total loss is 1.5 dB + x = x ./ sqrt(10 ^ 0.15); + y = y ./ sqrt(10 ^ 0.15); + end + end + %% Final Kerr effect has overshot by half step, so cancel this + P = abs(x) .^ 2 + abs(y) .^ 2; + x = x .* exp(-P .* hhz / 2 - attn / 2); + y = y .* exp(-P .* hhz / 2 - attn / 2); +end diff --git a/txFilter.m b/txFilter.m new file mode 100644 index 0000000..6af1ea1 --- /dev/null +++ b/txFilter.m @@ -0,0 +1,28 @@ +function x = txFilter(modData, rolloff, span, sps) + %% Transmitter pulse-shaping (root raised cosine) filter. + %% Inputs: + %% - modData: modulated data + %% - rolloff: rolloff factor in root raised cosine filter. + %% - span: filter span (number of symbols) + %% - sps: samples per symbol + %% Output: + %% - x: pulse-shaped waveform + + %% Construct filter object + txfilter = comm.RaisedCosineTransmitFilter... + ('Shape', 'Square root', ... + 'RolloffFactor', rolloff, ... + 'FilterSpanInSymbols', span, ... + 'OutputSamplesPerSymbol', sps, ... + 'Gain', sqrt(sps)); % so that output has energy 1 + + %% Extract filter coefficients + coef = coeffs(txfilter); + %% Upsample data and perform filtering in frequency domain + filter_fft = fft(coef.Numerator, length(modData) * sps); + modData_fft = fft(upsample(modData, sps)); + x = ifft(modData_fft .* filter_fft.'); + %% Re-order signal due to circular convolution + l = (length(coef.Numerator) - 1) / 2; + x = [x(l:end); x(1:l-1)]; +end -- 2.7.4