--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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
projects / 100GbE-PON.git / commitdiff