kstart = -floor(N/2);
% h: FIR filter
- h = sqrt(j * c * Tsamp^2 / (D * lambda^2 * z)) * ...
- exp(-j * pi * c * Tsamp^2 * k .^ 2 / (D * lambda^2 * z));
+ h = exp(-j * pi * c * Tsamp^2 * k .^ 2 / (D * lambda^2 * z));
+
+ len_fft = max(length(y), length(h));
+ H = fft(h, len_fft);
+ Y = fft(y, len_fft);
+
+ yCDComp = ifft(H.' .* Y);
+ l = (length(h) - 1) / 2;
+ if l > 0
+ yCDComp = [yCDComp(l:end); yCDComp(1:l-1)];
+ else
+ yCDComp = [yCDComp(end); yCDComp(1:end-1)];
+ end
- yCDComp = upfirdn(y, h);
- yCDComp = yCDComp(N:end); % truncate filter transients
end
%% References
%% [1]: S.J. Savory, Digital filters for coherent optical receivers, 2008.
-numSymbs = 5e5;
+numSymbs = 2^16;
M = 4;
Rsym = 2.5e10; % symbol rate (sym/sec)
rolloff = 0.25;
span = 6; % filter span
-sps = 2; % samples per symbol
+sps = 8; % samples per symbol
fs = Rsym * sps; % sampling freq (Hz)
Tsamp = 1 / fs;
t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).';
-EbN0_db = 0:0.2:14;
+EbN0_db = 0:0.5:14;
EbN0 = 10 .^ (EbN0_db ./ 10);
Es = 1;
%% Simulate chromatic dispersion
D = 17; % ps / (nm km)
lambda = 1550; % nm
-z = 60;%000; % km
+z = 3000; % km
-usingFFT = 1
-xCD = chromaticDispersion_FFT(x, D, lambda, z, Tsamp);
-%%xCD = normalizeEnergy(xCD, numSymbs, 1);
-%%xCD = x;
+
+[xCD, xCDkstart] = chromaticDispersion(x, D, lambda, z, Tsamp);
+
+TsampOrig = Tsamp;
for i = 1:plotlen
+ sps = 8;
+
snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps);
- noiseEnergy = 10 ^ (-snr / 10);
y = awgn(xCD, snr, 'measured');
- %%y = xCD;
r = rxFilter(y, rolloff, span, sps);
- rCDComp = CDCompensation(r, D, lambda, z, Tsamp);
+
+ sps = 2;
+ Tsamp = TsampOrig * 4;
+
+ [rCDComp, CDCompkstart] = CDCompensation(r, D, lambda, z, Tsamp);
rCDComp = normalizeEnergy(rCDComp, numSymbs*sps, 1);
- rSampled = rCDComp(sps*span/2+1:sps:(numSymbs+span/2)*sps);
- rNoCompSampled = r(sps*span/2+1:sps:(numSymbs+span/2)*sps);
+ rSampled = rCDComp(2:2:end);
+ rNoCompSampled = r(2:2:end);
%% rotate rNoCompSampled to match original data
theta = angle(-sum(rNoCompSampled .^ M)) / M;
-function adaptFilterOut = adaptiveCMA(rSampled)
+function [adaptFilterOut, convergeIdx] = adaptiveCMA(rSampled)
%% adaptive filter
%% CMA
taps = 19; % ODD taps
adaptFilterOut = zeros(numSymbs, 1);
converged = 0;
convergeCount = 0;
+ convergeIdx = Inf;
for it = 1:numSymbs
if it <= (taps - 1) / 2;
xout = sum(hxx .* xp);
ex = 1 - abs(xout) ^ 2;
- if abs(ex) < 1e-3
+ if abs(ex) < 0.01
convergeCount = convergeCount + 1;
else
convergeCount = 0;
end
if ~converged && convergeCount >= 10
- converged = 1
- it
+ converged = 1;
+ convergeIdx = it;
end
adaptFilterOut(it) = xout;
+++ /dev/null
-distance = input('Enter fiber length in L_D ');
-beta2 = input('dispersion: 1 for normal, -1 for anomalous ');
-N = input('Nonlinear parameter N = '); % Soliton order
-mshape = input('m = 0 for sech, m > 0 for super-Gaussian ');
-chirp0 = 0;
-
-% set simulation parameters
-nt = 1024; Tmax = 32; % FFT points and window size
-step_num = round(20 * distance * N^2); % No. of z steps
-deltaz = distance/step_num; % step size in z
-dtau = (2*Tmax) / nt; % step size in tau
-
-%% tau and omega arrays
-tau = (-nt/2 : nt/2-1) * dtau; % temporal grid
-omega = (pi/Tmax) * [(0:nt/2-1) (-nt/2:-1)]; % freq grid
-
-if mshape == 0
- uu = sech(tau) .* exp(-0.5j * chirp0 * tau.^2);
-else
- uu = exp(-0.5 * (1 + 1j * chirp0) .* tau.^(2 * mshape));
-end
-
-%% plot input pulse shape and spectrum
-temp = fftshift(ifft(uu)) .* (nt * dtau) / sqrt(2 * pi); % spectrum
-figure(1); clf; subplot(2,1,1);
-plot(tau, abs(uu).^2, '--k'); hold on;
-axis([-5 5 0 inf]);
-xlabel('Normalized Time');
-ylabel('Normalized Power');
-title('Input and Output pulse shape and spectrum');
-
-subplot(2, 1, 2);
-plot(fftshift(omega)/(2*pi), abs(temp) .^ 2, '--k'); hold on;
-axis([-.5 .5 0 inf]);
-xlabel('Normaized freq');
-ylabel('spectral power');
-
-%% store dispersive phase shifts to speed up code
-dispersion = exp(0.5j * beta2 * omega.^2 * deltaz); % [hase factor
-hhz = 1j * N^2 * deltaz;
-
-%% begin main loop
-%% N/2 -> D -> N/2 first half step nonlinear
-temp = uu .* exp(abs(uu) .^ 2 .* hhz / 2);
-for n = 1 : step_num
- f_temp = ifft(temp) .* dispersion;
- uu = fft(f_temp);
- temp = uu .* exp(abs(uu) .^ 2 .* hhz);
-end
-uu = temp .* exp(-abs(uu) .^ 2 .* hhz); % final field
-temp = fftshift(ifft(uu)) .* (nt*dtau) / sqrt(2 * pi); % final spectrum
-%% end of main loop
-
-%% plot output pulse shape and spectrum
-subplot(2,1,1);
-plot(tau, abs(uu).^2, '-k');
-subplot(2,1,2);
-plot(fftshift(omega) / (2*pi), abs(temp).^2, '-k');
+++ /dev/null
-function baseband(rolloff, M, numSymbs)
- %% Set defaults for inputs
- if nargin < 3
- numSymbs = 1000;
- end
- if nargin < 2
- M = 2;
- end
- if nargin < 1
- rolloff = 0.5;
- end
-
- %% https://www.mathworks.com/help/comm/examples/passband-modulation-with-adjacent-channel-interference.html
- Rsym = 1e6; % symbol rate (sym/sec)
-
- span = 6; % filter span
- sps = 2; % samples per symbol
-
- fs = Rsym * sps; % sampling freq (Hz)
-
- t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs)';
-
- EbN0_db = 0:0.2:10;
- EbN0 = 10 .^ (EbN0_db ./ 10);
- Es = 1;
- Eb = Es / log2(M);
- N0 = Eb ./ EbN0;
-
- EsN0 = EbN0 .* log2(M);
- EsN0_db = 10 .* log10(EsN0);
-
- plotlen = length(EbN0);
- ber = zeros(1, plotlen);
-
- data = randi([0 M - 1], numSymbs, 1);
- modData = pskmod(data, M, 0, 'gray');
-
- xBaseband = txFilter(modData, rolloff, span, sps);
-
- for i = 1:plotlen
- snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps);
- noiseEnergy = 10 ^ (-snr / 10);
-
- yBaseband = awgn(xBaseband, snr, 'measured');
-
- rBaseband = rxFilter(yBaseband, rolloff, span, sps);
-
- rSampled = rBaseband(sps*span/2+1:sps:(numSymbs+span/2)*sps);
- demodData = pskdemod(rSampled, M, 0, 'gray');
- [bitErrors, ber(i)] = biterr(data, demodData);
- end
-
- figure(1);
- clf;
-
- %% Plot simulated results
- semilogy(EbN0_db, ber, 'r', 'LineWidth', 2);
- hold on;
-
- theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1);
- legend('Simulated', 'Discrete');
-
- title(strcat(num2str(M), '-PSK with Gray code'));
- grid on;
- xlabel('$E_b/N_0$ (dB)');
- ylabel('BER');
-
- formatFigure;
-end
--- /dev/null
+numSymbs = 2^16;
+M = 4;
+
+Rsym = 2.5e10; % symbol rate (sym/sec)
+Tsym = 1 / Rsym; % symbol period (sec)
+
+rolloff = 0.25;
+span = 6; % filter span
+sps = 8; % samples per symbol
+
+fs = Rsym * sps; % sampling freq (Hz)
+Tsamp = 1 / fs;
+
+t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).';
+
+power_dBm = -6:1:4;
+%%power_dBm = 0;
+power = 10 .^ (power_dBm / 10) * 1e-3; % watts
+
+Es = power * Tsym; % joules
+Eb = Es / log2(M); % joules
+
+N0ref_db = 10; % Eb/N0 at power = 1mW
+%% Fix N0, such that Eb/N0 = N0ref_db at power = 1mW
+N0 = 1e-3 * Tsym / (log2(M) * 10 ^ (N0ref_db / 10)); % joules
+%% At current settings, N0 = 0.002 pJ
+
+plotlen = length(power);
+
+ber = zeros(1, plotlen);
+
+data = randi([0 M - 1], numSymbs, 1);
+%%modData = dpskmod(data, M, 0, 'gray');
+modData = pskmod(data, M, 0, 'gray');
+for i = 2:numSymbs
+ modData(i) = modData(i) * modData(i-1);
+end
+
+
+%% Chromatic dispersion
+D = 17; % ps / (nm km)
+lambda = 1550; % nm
+z = 100; % km
+
+
+linewidthTx = 0; % Hz
+linewidthLO = 1e6; % Hz
+
+
+TsampOrig = Tsamp;
+
+x_P1 = txFilter(modData, rolloff, span, sps);
+
+
+for i = 1:plotlen
+ sps = 8;
+ Tsamp = TsampOrig;
+
+ snr = Es(i) / sps / N0;
+ snr_dB = 10 * log10(snr);
+
+ %%x = txFilter(modData, rolloff, span, sps);
+ %% Now, sum(abs(x) .^ 2) / length(x) should be 1.
+ %% We can set its power simply by multiplying.
+ x = sqrt(power(i)) * x_P1;
+
+ %% We can now do split-step Fourier.
+ gamma = 1.2; % watt^-1 / km
+
+
+ xCDKerr = splitstepfourier(x, D, lambda, z, Tsamp, gamma);
+
+ xpn = phaseNoise(xCDKerr, linewidthTx, linewidthLO, Tsamp);
+
+ y = awgn(xpn, snr_dB, 'measured', 'db');
+ %y = xCDKerr;
+
+ r = rxFilter(y, rolloff, span, sps);
+ sps = 2;
+ Tsamp = Tsamp * 4;
+
+ rCDComp = CDCompensation(r, D, lambda, z, Tsamp);
+ rCDComp = normalizeEnergy(rCDComp, numSymbs * sps, 1);
+
+ rSampled = rCDComp(2:2:end);
+
+ %% adaptive filter
+ [adaptFilterOut, convergeIdx] = adaptiveCMA(rSampled);
+
+ pncorr = phaseNoiseCorr(adaptFilterOut, M, 0, 40).';
+
+ demodAdapt = pskdemod(pncorr, M, 0, 'gray');
+ remod = pskmod(demodAdapt, M, 0, 'gray');
+ delayed = [1; remod(1:end-1)];
+ demod = pskdemod(remod .* conj(delayed), M, 0, 'gray');
+
+ if convergeIdx < Inf
+ [~, ber(i)] = biterr(data(convergeIdx:end), demod(convergeIdx:end));
+ else
+ [~, ber(i)] = biterr...
+ (data(ceil(0.8*numSymbs):end), ...
+ demod(ceil(0.8*numSymbs):end));
+ end
+end
+
+ber
+
+
+figure(1);
+clf;
+
+%% Plot simulated results
+qp = 20 * log10(erfcinv(2*ber)*sqrt(2));
+plot(power_dBm, qp, 'Color', [0, 0.6, 0], 'LineWidth', 2);
+hold on;
+
+title({'CD + Kerr + CD compensation', ...
+ strcat(['$D = 17$ ps/(nm km), $z = ', num2str(z), '$ km'])});
+grid on;
+xlabel('Optical power (dBm)');
+ylabel('$20 \log_{10}\left(\sqrt{2}\mathrm{erfc}^{-1}(2 BER)\right)$');
+
+formatFigure;
t = k * Tsamp;
% Impulse response
- g = sqrt(c / (j * D * lambda^2 * z)) * ...
- exp(j * pi * c / (D * lambda^2 * z) * t .^ 2);
+ g = exp(j * pi * c / (D * lambda^2 * z) * t .^ 2);
- xCD = conv(x, g);
+ lenx = length(x);
+ leng = length(g);
+
+ len_fft = max(lenx, leng);
+
+ G = fft(g, len_fft);
+ X = fft(x, len_fft);
+
+ xCD = ifft(G.' .* X);
+ l = (leng - 1) / 2;
+ if l > 0
+ xCD = [xCD(l:end); xCD(1:l-1)];
+ end
kstart = 1 - kmax;
end
span = 6; % Tx/Rx filter span
rolloff = 0.25; % Tx/Rx RRC rolloff
-sps = 2; % samples per symbol
+sps = 8; % samples per symbol
fs = Rsym * sps; % sampling freq (Hz)
Tsamp = 1 / fs;
%% Simulate chromatic dispersion
D = 17; % ps / (nm km)
lambda = 1550; % nm
-z = 5000 % km
+z = 500 % km
-[xCD, xCDkstart] = chromaticDispersion_FFT(x, D, lambda, z, Tsamp);
+[xCD, xCDkstart] = chromaticDispersion(x, D, lambda, z, Tsamp);
EbN0_db = 8;
snr = EbN0_db + 10 * log10(log2(M)) - 10 * log10(sps);
y = xCD;
r = rxFilter(y, rolloff, span, sps);
+
+sps = 2;
+Tsamp = Tsamp * 4;
+
+
[rCDComp, CDCompkstart] = CDCompensation(r, D, lambda, z, Tsamp);
rCDComp = normalizeEnergy(rCDComp, numSymbs*sps, 1);
-rSampled = rCDComp(sps*span/2+1:sps:(numSymbs+span/2)*sps);
-rNoCompSa = r(sps*span/2+1:sps:(numSymbs+span/2)*sps);
+rSampled = rCDComp(2:2:end);
+rNoCompSa = r(2:2:end);
%% if no CD comp, then rotate constellation. Use:
theta = angle(-sum(rNoCompSa .^ M)) / M;
rSampled = rSampled .* exp(-1j * theta);
-
%%rAdaptEq = adaptiveCMA(rSampled);
%{
%% Compare original signal and compensated signal
+++ /dev/null
-Rsym = 2.5e10; % symbol/sec
-sps = 4; % samples/symbol
-fs = Rsym * sps; % sampling freq (Hz)
-Tsamp = 1 / fs;
-
-x = 10 * normpdf(1:30, 10, 3);
-
-xCD = chromaticDispersion(x, 20, 1550, 30, Tsamp);
-
-figure(10001);
-subplot(121);
-plot(abs(x) .^ 2);
-axis([0 30 -Inf Inf]);
-title('Before dispersion');
-ylabel('Power (arb. unit)');
-xlabel('Time (arb. unit)');
-subplot(122);
-plot(abs(xCD) .^ 2)
-axis([20 50 -Inf Inf]);
-title('After dispersion');
-xlabel('Time (arb.~unit)');
-ylabel('Power (arb.~unit)');
-
-subplot(121);
-formatFigure;
-subplot(122)
-formatFigure;
-set(figure(10001), 'Units', 'centimeters', ...
- 'Position', [0 0 24 8], 'PaperPositionMode', 'auto');
-subplot(121);
-set(gca, 'XTick', 0:5:30);
-subplot(122);
-set(gca, 'XTick', 20:5:50);
-
-%saveas(gcf, 'dispersionTest.png');
+++ /dev/null
-function [xCD, kstart] = chromaticDispersion_FFT(x, D, lambda, z, Tsamp)
- %% Simulate chromatic dispersion.
- %% Params:
- %% - x: input waveform (pulse-shaped)
- %% - D: dispersion coefficient (ps / (nm km))
- %% - lambda: wavelength (nm)
- %% - z: length of fibre (km)
- %% - Tsamp: sampling time (s)
- %% Output:
- %% - xCD: x after being dispersed. Energy of xCD is not normalized.
- %% - kstart: starting index of the discrete signal
-
- %% 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
-
- %% Time domain filter length, needed for compatibility with
- %% time-domain (convolution) chromatic dispersion.
- kmax = floor(abs(D) * lambda^2 * z / (2 * c * Tsamp^2));
- kstart = 1 - kmax;
-
- x = [zeros(kmax, 1); x]; % prepend zeros to allow for transients
-
- xDFT = fft(x);
- n = length(x);
- fs = 1 / Tsamp;
-
- omega = (2*pi * fs / n * [(0 : floor((n-1)/2)), (-ceil((n-1)/2) : -1)]).';
- dispDFT = xDFT .* exp(-1j * omega.^2 * D * lambda^2 * z / (4 * pi * c));
-
- xCD = ifft(dispDFT);
-
- if ceil(kmax/2) > 0
- %% fix the order of the samples due to prepending zeros before FFT
- xCD = [xCD(ceil(kmax/2):end); xCD(1:ceil(kmax/2)-1)];
- end
- %% pad zeros for compatibility with time-domain filter
- xCD = [zeros(floor((kmax-1)/2), 1); xCD; zeros(ceil((kmax+1)/2), 1)];
-end
-%% References
-%% [1]: S.J. Savory, Digital filters for coherent optical receivers, 2008.
+++ /dev/null
-function discretePSK_BER_SNR(M, numSymbs)
- %% Set defaults for inputs
- if nargin < 2
- numSymbs = 1000;
- end
- if nargin < 1
- M = 2;
- end
-
- EbN0_db = 0:0.2:10;
- EbN0 = 10 .^ (EbN0_db ./ 10);
-
- Es = 1;
- Eb = Es / log2(M);
- N0 = Eb ./ EbN0;
-
- EsN0 = EbN0 .* log2(M);
- EsN0_db = 10 .* log10(EsN0);
-
- plotlen = length(EbN0);
-
- ber = zeros(1, plotlen);
-
- data = randi([0, M - 1], 1, numSymbs);
- txsig = pskmod(data, M, 0, 'gray');
-
- for i = 1:plotlen
- rxsig = awgn(txsig, EsN0_db(i));
- demodData = pskdemod(rxsig, M, 0, 'gray');
-
- [bitErrors, ber(i)] = biterr(data, demodData);
- end
-
- figure(1);
- clf;
-
- %% Plot simulated results
- semilogy(EbN0_db, ber, 'r', 'LineWidth', 2);
- hold on;
-
- theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1);
-
- if M == 2 || M == 4
- legend('Simulated', 'Theoretical');
- else
- legend('Simulated', 'Approximation');
- end
-
- title(strcat(num2str(M), '-PSK with Gray code'));
- grid on;
- xlabel('$E_b/N_0$ (dB)');
- ylabel('BER');
-
- formatFigure;
- %saveas(gcf, strcat('BER_SNR_', num2str(M), 'PSK_', num2str(numSymbs), ...
- % '.svg'));
-end
set(gca, 'FontSize', 14);
leg = findobj(gcf, 'Type', 'Legend');
set(leg, 'FontSize', 14);
+set(leg, 'Interpreter', 'Latex');
+tit = get(gca, 'title');
+set(tit, 'Interpreter', 'Latex');
+xl = get(gca, 'XLabel');
+set(xl, 'Interpreter', 'Latex');
+yl = get(gca, 'YLabel');
+set(yl, 'Interpreter', 'Latex');
-numSymbs = 5e4;
+numSymbs = 2^16;
M = 4;
Rsym = 2.5e10; % symbol rate (sym/sec)
rolloff = 0.25;
span = 6; % filter span
-sps = 2; % samples per symbol
+sps = 8; % samples per symbol
fs = Rsym * sps; % sampling freq (Hz)
Tsamp = 1 / fs;
t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).';
-
-power_dBm = -3:0.25:4;
+power_dBm = -6:1:4;
+%%power_dBm = 0;
power = 10 .^ (power_dBm / 10) * 1e-3; % watts
Es = power * Tsym; % joules
Eb = Es / log2(M); % joules
-N0ref_dB = 10; % Eb/N0 at power = 1mW
-%% Fix N0, such that Eb/N0 = N0ref_dB at power = 1mW
-N0 = 1e-3 * Tsym / (log2(M) * 10 ^ (N0ref_dB / 10)); % joules
-
+N0ref_db = 10; % Eb/N0 at power = 1mW
+%% Fix N0, such that Eb/N0 = N0ref_db at power = 1mW
+N0 = 1e-3 * Tsym / (log2(M) * 10 ^ (N0ref_db / 10)); % joules
+%% At current settings, N0 = 0.002 pJ
plotlen = length(power);
ber = zeros(1, plotlen);
data = randi([0 M - 1], numSymbs, 1);
-modData = pskmod(data, M, pi / M, 'gray');
+modData = dpskmod(data, M, 0, 'gray');
+%%modData = pskmod(data, M, pi/4, 'gray');
%% Chromatic dispersion
D = 17; % ps / (nm km)
lambda = 1550; % nm
-z = 600; % km
+z = 100; % km
+TsampOrig = Tsamp;
+
for i = 1:plotlen
+ sps = 8;
+ Tsamp = TsampOrig;
+
snr = Es(i) / sps / N0;
snr_dB = 10 * log10(snr);
%% We can now do split-step Fourier.
gamma = 1.2; % watt^-1 / km
- %%stepnum = round(40 * z * gamma); % Nonlinear Fiber optics, App B
- stepnum = 100;
- xCD = splitstepfourier(x, D, lambda, z, Tsamp, gamma, stepnum);
- y = awgn(xCD, snr, power(i), 'linear');
- %%y = xCD;
+
+ xCDKerr = splitstepfourier(x, D, lambda, z, Tsamp, gamma);
+
+ y = awgn(xCDKerr, snr_dB, 'measured', 'db');
+ %y = xCDKerr;
r = rxFilter(y, rolloff, span, sps);
+ sps = 2;
+ Tsamp = Tsamp * 4;
+
rCDComp = CDCompensation(r, D, lambda, z, Tsamp);
- rCDComp = normalizeEnergy(rCDComp, numSymbs*sps, 1);
+ rCDComp = normalizeEnergy(rCDComp, numSymbs * sps, 1);
- rSampled = rCDComp(sps*span/2+1:sps:(numSymbs+span/2)*sps);
- rNoCompSampled = r(sps*span/2+1:sps:(numSymbs+span/2)*sps);
+ rSampled = rCDComp(2:2:end);
- %% rotate rNoCompSampled to match original data
- theta = angle(-sum(rNoCompSampled .^ M)) / M;
- %% if theta approx +pi/M, wrap to -pi/M
- if abs(theta - pi / M) / (pi / M) < 0.1
- theta = -pi / M;
- end
- rNoCompSampled = rNoCompSampled .* exp(-j * theta);
+ %% adaptive filter
+ [adaptFilterOut, convergeIdx] = adaptiveCMA(rSampled);
+ demod = dpskdemod(adaptFilterOut, M, 0, 'gray');
+ %%demod = pskdemod(adaptFilterOut, M, pi/4, 'gray');
- %% 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;
+ if convergeIdx < Inf
+ [~, ber(i)] = biterr(data(convergeIdx:end), demod(convergeIdx:end));
+ else
+ [~, ber(i)] = biterr...
+ (data(ceil(0.8*numSymbs):end), ...
+ demod(ceil(0.8*numSymbs):end));
end
- rSampled = rSampled .* exp(-1j * theta);
-
+end
- %% adaptive filter
- adaptFilterOut = adaptiveCMA(rSampled);
+ber
- demodAdapt = pskdemod(adaptFilterOut, M, pi / M, 'gray');
- [~, ber(i)] = biterr(data, demodAdapt);
-end
-figure(1);
+figure;
clf;
%% Plot simulated results
-semilogy(power_dBm, ber, 'Color', [0, 0.6, 0], 'LineWidth', 2);
+qp = 20 * log10(erfcinv(2*ber)*sqrt(2));
+plot(power_dBm, qp, 'Color', [0, 0.6, 0], 'LineWidth', 2);
hold on;
title({'CD + Kerr + CD compensation', ...
- strcat(['$D = 17$ ps/(nm km), $z = ', num2str(z), '$ km']), ...
- strcat(['$E_b/N_0 = ', num2str(N0ref_dB), '$ dB at 1 mW'])});
+ strcat(['$D = 17$ ps/(nm km), $z = ', num2str(z), '$ km'])});
grid on;
xlabel('Optical power (dBm)');
-ylabel('BER');
+ylabel('$20 \log_{10}\left(\sqrt{2}\mathrm{erfc}^{-1}(2 BER)\right)$');
formatFigure;
+++ /dev/null
-numSymbs = 5e5;
-M = 4;
-
-Rsym = 2.5e10; % symbol rate (sym/sec)
-Tsym = 1 / Rsym; % symbol period (sec)
-
-rolloff = 0.25;
-span = 6; % filter span
-sps = 2; % samples per symbol
-
-fs = Rsym * sps; % sampling freq (Hz)
-Tsamp = 1 / fs;
-
-t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs).';
-
-
-%%power_dBm = -3:0.2:4;
-power_dBm = [0];
-power = 10 .^ (power_dBm / 10) * 1e-3; % watts
-
-Es = power * Tsym; % joules
-Eb = Es / log2(M); % joules
-
-N0ref_db = 10; % Eb/N0 at power = 1mW
-%% Fix N0, such that Eb/N0 = N0ref_db at power = 1mW
-N0 = 1e-3 * Tsym / (log2(M) * 10 ^ (N0ref_db / 10)); % joules
-
-
-plotlen = length(power);
-
-ber = zeros(1, plotlen);
-
-data = randi([0 M - 1], numSymbs, 1);
-modData = pskmod(data, M, pi / M, 'gray');
-
-
-%% Chromatic dispersion
-D = 17; % ps / (nm km)
-lambda = 1550; % nm
-z = 600; % km
-
-
-for i = 1:plotlen
- snr = Es(i) / sps / N0;
- snr_dB = 10 * log10(snr);
-
- x = txFilter(modData, rolloff, span, sps);
- %% Now, sum(abs(x) .^ 2) / length(x) should be 1.
- %% We can set its power simply by multiplying.
- x = sqrt(power(i)) * x;
-
- %% We can now do split-step Fourier.
- gamma = 1.2; % watt^-1 / km
- %%stepnum = round(40 * z * gamma); % Nonlinear Fiber optics, App B
- stepnum = 100;
- xCD = splitstepfourier(x, D, lambda, z, Tsamp, gamma, stepnum);
-
- y = awgn(xCD, snr, power(i), 'linear');
- %%y = xCD;
-
- r = rxFilter(y, rolloff, span, sps);
- rCDComp = CDCompensation(r, D, lambda, z, Tsamp);
- rCDComp = normalizeEnergy(rCDComp, numSymbs*sps, 1);
-
- rSampled = rCDComp(sps*span/2+1:sps:(numSymbs+span/2)*sps);
- rNoCompSampled = r(sps*span/2+1:sps:(numSymbs+span/2)*sps);
-
- %% rotate rNoCompSampled to match original data
- theta = angle(-sum(rNoCompSampled .^ M)) / M;
- %% if theta approx +pi/M, wrap to -pi/M
- if abs(theta - pi / M) / (pi / M) < 0.1
- theta = -pi / M;
- end
- rNoCompSampled = rNoCompSampled .* 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);
-
- %% adaptive filter
- adaptFilterOut = adaptiveCMA(rSampled);
-
- demodAdapt = pskdemod(adaptFilterOut, M, pi / M, 'gray');
- [~, ber(i)] = biterr(data, demodAdapt)
-end
-
-return
-
-
-figure(1);
-clf;
-
-%% Plot simulated results
-semilogy(power_dBm, ber, 'Color', [0, 0.6, 0], 'LineWidth', 2);
-hold on;
-
-title({'CD + Kerr + CD compensation', ...
- strcat(['$D = 17$ ps/(nm km), $z = ', num2str(z), '$ km'])});
-grid on;
-%%xlabel('$E_b/N_0$ (dB)');
-xlabel('Optical power (dBm)');
-ylabel('BER');
-
-formatFigure;
+++ /dev/null
-function passband(rolloff, M, numSymbs)
- %% Set defaults for inputs
- if nargin < 3
- numSymbs = 1000;
- end
- if nargin < 2
- M = 2;
- end
- if nargin < 1
- rolloff = 0.5;
- end
-
- %% https://www.mathworks.com/help/comm/examples/passband-modulation-with-adjacent-channel-interference.html
- Rsym = 2.5e10; % symbol rate (sym/sec)
-
- span = 6; % filter span
- sps = 4; % samples per symbol
-
- fs = Rsym * sps; % sampling freq (Hz)
-
- t = (0 : 1 / fs : numSymbs / Rsym + (1.5 * span * sps - 1) / fs)';
-
-
- EbN0_db = 0:0.2:10;
- EbN0 = 10 .^ (EbN0_db ./ 10);
-
- Es = 1;
- Eb = Es / log2(M);
- N0 = Eb ./ EbN0;
-
- EsN0 = EbN0 .* log2(M);
- EsN0_db = 10 .* log10(EsN0);
-
- plotlen = length(EbN0);
-
- ber = zeros(1, plotlen);
-
-
- data = randi([0 M - 1], numSymbs, 1);
- modData = pskmod(data, M, 0, 'gray');
-
- xBaseband = txFilter(modData, rolloff, span, sps);
-
- fc = 1e12; % Carrier freq (Hz)
- carrier = sqrt(2) * exp(j * 2 * pi * fc * t);
-
- xPassband = xBaseband .* carrier(1:length(xBaseband));
-
- for i = 1:plotlen
- snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps);
- noiseEnergy = 10 ^ (-snr / 10);
-
- yPassband = awgn(xPassband, snr, 'measured');
- yLO = yPassband .* conj(carrier(1:length(yPassband)));
- rBaseband = rxFilter(yLO, rolloff, span, sps);
-
- rSampled = rBaseband(sps*span/2+1:sps:(numSymbs + span/2) * sps);
-
- demodData = pskdemod(rSampled, M, 0, 'gray');
- [bitErrors, ber(i)] = biterr(data, demodData);
- end
-
- figure(1);
- clf;
-
- %% Plot simulated results
- semilogy(EbN0_db, ber, 'r', 'LineWidth', 2);
- hold on;
-
- theoreticalPSK(EbN0_db, M, 'b', 'LineWidth', 1);
- legend('Simulated RRC', 'Discrete');
-
- title(strcat(num2str(M), '-PSK RRC with Gray code'));
- grid on;
- xlabel('$E_b/N_0$ (dB)');
- ylabel('BER');
-
- formatFigure;
-end
--- /dev/null
+function [x, y] = pdm_adaptiveCMA(rx, ry)
+ taps = 19; % ODD taps
+ 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;
+
+ mu = 1e-3;
+ numSymbs = length(rx);
+
+ %% Check average energy of symbols
+ rx = rx / sqrt(mean(abs(rx) .^ 2));
+ ry = ry / sqrt(mean(abs(ry) .^ 2));
+
+ x = zeros(numSymbs, 1);
+ y = zeros(numSymbs, 1);
+ converged = 0;
+ convergeCount = 0;
+ convergeIdx = Inf;
+
+ for it = 1:numSymbs
+ if it <= (taps - 1) / 2;
+ 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
+ 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
+ xp = rx(it - (taps - 1) / 2 : it + (taps - 1) / 2);
+ yp = ry(it - (taps - 1) / 2 : it + (taps - 1) / 2);
+ end
+
+ xout = sum(hxx .* xp) + sum(hxy .* yp);
+ yout = sum(hyx .* xp) + sum(hyy .* yp);
+ ex = 1 - abs(xout) ^ 2;
+ ey = 1 - abs(yout) ^ 2;
+
+ if abs(ex) < 0.01
+ convergeCount = convergeCount + 1;
+ else
+ convergeCount = 0;
+ end
+ if ~converged && convergeCount >= 10
+ converged = 1;
+ convergeIdx = it;
+ end
+
+ x(it) = xout;
+ y(it) = yout;
+
+ 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 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
--- /dev/null
+numSymbs = 2^14;
+M = 4;
+
+Rsym = 2.5e10; % symbol rate (sym/sec)
+Tsym = 1 / Rsym; % symbol period (sec)
+
+rolloff = 0.25;
+span = 6; % filter span
+sps = 8; % samples per symbol
+
+fs = Rsym * sps; % sampling freq (Hz)
+Tsamp = 1 / fs;
+
+t = (0 : 1 / fs : numSymbs/2 / Rsym - 1/fs).';
+
+power_dBm = -6:1:4;
+%%power_dBm = 0;
+power = 10 .^ (power_dBm / 10) * 1e-3; % watts
+
+Es = power * Tsym; % joules
+Eb = Es / log2(M); % joules
+
+N0ref_db = 10; % Eb/N0 at power = 1mW
+%% Fix N0, such that Eb/N0 = N0ref_db at power = 1mW
+N0 = 1e-3 * Tsym / (log2(M) * 10 ^ (N0ref_db / 10)); % joules
+%% At current settings, N0 = 0.002 pJ
+
+plotlen = length(power);
+
+ber = zeros(1, plotlen);
+
+data = randi([0 M - 1], numSymbs, 1);
+%%modData = dpskmod(data, M, 0, 'gray');
+modData = pskmod(data, M, 0, 'gray');
+for i = 2:numSymbs
+ modData(i) = modData(i) * modData(i-1);
+end
+
+
+%% Chromatic dispersion
+D = 17; % ps / (nm km)
+lambda = 1550; % nm
+z = 100; % km
+
+linewidthTx = 0; % Hz
+linewidthLO = 1e6; % Hz
+
+TsampOrig = Tsamp;
+
+sig_x = txFilter(modData(1:numSymbs/2), rolloff, span, sps);
+sig_y = txFilter(modData(numSymbs/2+1:end), rolloff, span, sps);
+
+for i = 1:plotlen
+ sps = 8;
+ Tsamp = TsampOrig;
+
+ snr = Es(i) / sps / N0;
+ snr_dB = 10 * log10(snr);
+
+ %%x = txFilter(modData, rolloff, span, sps);
+ %% Now, sum(abs(x) .^ 2) / length(x) should be 1.
+ %% We can set its power simply by multiplying.
+ %%x = sqrt(power(i)) * x;
+ txx = sig_x * sqrt(power(i));
+ txy = sig_y * sqrt(power(i));
+
+ rot_omega = 1e3; % rad/s
+ rot_phi = 2;
+ rot_x = txx .* cos(rot_omega * t) + ...
+ txy .* sin(rot_omega * t) * exp(-1j * rot_phi);
+ rot_y = txx .* -sin(rot_omega * t) * exp(1j * rot_phi) + ...
+ txy .* cos(rot_omega * t);
+
+ %% We can now do split-step Fourier.
+ gamma = 1.2; % watt^-1 / km
+
+ [xCDKerr, yCDKerr] = ssf_pdm(rot_x, rot_y, ...
+ D, lambda, z, Tsamp, gamma);
+
+ xpn = phaseNoise(xCDKerr, linewidthTx, linewidthLO, Tsamp);
+ ypn = phaseNoise(yCDKerr, linewidthTx, linewidthLO, Tsamp);
+
+ xout = awgn(xpn, snr_dB, 'measured', 'db');
+ yout = awgn(ypn, snr_dB, 'measured', 'db');
+
+ rx = rxFilter(xout, rolloff, span, sps);
+ ry = rxFilter(yout, rolloff, span, sps);
+ sps = 2;
+ Tsamp = Tsamp * 4;
+
+ rxCDComp = CDCompensation(rx, D, lambda, z, Tsamp);
+ ryCDComp = CDCompensation(ry, D, lambda, z, Tsamp);
+
+ rxSampled = rxCDComp(1:2:end);
+ rySampled = ryCDComp(1:2:end);
+
+ %% adaptive filter
+ [xCMA, yCMA] = pdm_adaptiveCMA(rxSampled, rySampled);
+
+ xpncorr = phaseNoiseCorr(xCMA, M, 0, 40).';
+ ypncorr = phaseNoiseCorr(yCMA, M, 0, 40).';
+
+ demodx = pskdemod(xpncorr, M, 0, 'gray');
+ remodx = pskmod(demodx, M, 0, 'gray');
+ delayx = [1; remodx(1:end-1)];
+ demodx = pskdemod(remodx .* conj(delayx), M, 0, 'gray');
+ clear remodx
+ clear delayx
+
+ demody = pskdemod(ypncorr, M, 0, 'gray');
+ remody = pskmod(demody, M, 0, 'gray');
+ delayy = [1; remody(1:end-1)];
+ demody = pskdemod(remody .* conj(delayy), M, 0, 'gray');
+ clear remody
+ clear delayy
+
+
+ [~, ber(i)] = biterr(data, [demodx; demody]);
+end
+
+ber
+
+figure;
+clf;
+
+%% Plot simulated results
+qp = 20 * log10(erfcinv(2*ber)*sqrt(2));
+plot(power_dBm, qp, 'Color', [0, 0.6, 0], 'LineWidth', 2);
+hold on;
+
+title({'CD + Kerr + CD compensation', ...
+ strcat(['$D = 17$ ps/(nm km), $z = ', num2str(z), '$ km'])});
+grid on;
+xlabel('Optical power (dBm)');
+ylabel('$20 \log_{10}\left(\sqrt{2}\mathrm{erfc}^{-1}(2 BER)\right)$');
+
+formatFigure;
Rsym = 2.5e10; % symbol rate (sym/sec)
rolloff = 0.5;
span = 6; % filter span
-sps = 4; % samples per symbol
+sps = 8; % samples per symbol
fs = Rsym * sps; % sampling freq (Hz)
Tsamp = 1 / fs;
snr = EbN0_db + 10 * log10(log2(M)) - 10 * log10(sps);
-%%y = awgn(xPN, snr, 'measured');
yPSK = awgn(xPSKpn, snr, 'measured');
yDPSK = xDPSKpn;
rPSK = rxFilter(yPSK, rolloff, span, sps);
rDPSK = rxFilter(yDPSK, rolloff, span, sps);
+sps = 2;
+Tsamp = Tsamp * 4;
-rPSKSa = rPSK(sps*span/2+1:sps:(numSymbs+span/2)*sps);
-rDPSKSa = rDPSK(sps*span/2+1:sps:(numSymbs+span/2)*sps);
+rPSKSa = rPSK(1:2:end);
+rDPSKSa = rDPSK(1:2:end);
-[rPSKSaPhEq, phiestsPSK] = phaseNoiseCorr(rPSKSa, M, pi/M, 40);
+[rPSKSaPhEq, phiestsPSK] = phaseNoiseCorr(rPSKSa, M, pi/M, 80);
demodPSK = pskdemod(rPSKSaPhEq, M, pi/M, 'gray');
demodDPSK = dpskdemod(rDPSKSa, M, pi/M, 'gray');
figure(2);
-plot(t(1:40000), repelem(-phiestsPSK, sps));
+plot(t(1:80000), repelem(-phiestsPSK, 8));
hold on;
-plot(t(1:40000), pTxLoPSK(1:40000));
+plot(t(1:80000), pTxLoPSK(1:80000));
legend('estimate', 'actual');
title('Phase noise estimation');
hold off;
-numSymbs = 5e3;
+numSymbs = 5e4;
M = 4;
rolloff = 0.5;
Rsym = 2.5e10; % symbol rate (sym/sec)
span = 6; % filter span
-sps = 4; % samples per symbol
+sps = 8; % samples per symbol
fs = Rsym * sps; % sampling freq (Hz)
Tsamp = 1 / fs;
linewidthLO = 5e6; % Hz
%linewidthLO = Rsym * 1e-3;
-iterations = 25;
+iterations = 1;
avgSa = 40;
+TsampOrig = Tsamp;
+
for it = 1 : iterations
[xPSKpn, pTxLoPSK] = phaseNoise(xPSK, linewidthTx, linewidthLO, Tsamp);
[xDEPSKpn, pTxLoDEPSK] = phaseNoise(xDEPSK, linewidthTx, linewidthLO, Tsamp);
[xDPSKpn, pTxLoDPSK] = phaseNoise(xDPSK, linewidthTx, linewidthLO, Tsamp);
for i = 1:plotlen
+ Tsamp = TsampOrig;
+ sps = 8;
+
snr = EbN0_db(i) + 10 * log10(log2(M)) - 10 * log10(sps);
noiseEnergy = 10 ^ (-snr / 10);
rDEPSK = rxFilter(yDEPSK, rolloff, span, sps);
rDPSK = rxFilter(yDPSK, rolloff, span, sps);
- rPSKSamp = rPSK(sps*span/2+1:sps:(numSymbs+span/2)*sps);
- rDEPSKSamp = rDEPSK(sps*span/2+1:sps:(numSymbs+span/2)*sps);
- rDPSKSamp = rDPSK(sps*span/2+1:sps:(numSymbs+span/2)*sps);
+ sps = 2;
+ Tsamp = TsampOrig * 4;
+
+ rPSKSamp = rPSK(1:2:end);
+ rDEPSKSamp = rDEPSK(1:2:end);
+ rDPSKSamp = rDPSK(1:2:end);
[rPSKSampEq, phiestsPSK] = phaseNoiseCorr(rPSKSamp, M, pi/M, avgSa);
[rDEPSKSampEq, phiestsDEPSK] = phaseNoiseCorr(rDEPSKSamp, M, 0, avgSa);
if EbN0_db(i) == 8 && it == 1
figure(1234);
- plot(repelem(-phiestsPSK, sps));
+ plot(repelem(-phiestsPSK, 8));
hold on;
plot(pTxLoPSK);
legend('estimate', 'actual');
%% Output:
%% - r: filtered signal (energy not normalized),
%% normalize and then sample.
- filter = comm.RaisedCosineReceiveFilter...
+ rxfilter = comm.RaisedCosineReceiveFilter...
('Shape', 'Square root', ...
'RolloffFactor', rolloff, ...
'FilterSpanInSymbols', span, ...
'InputSamplesPerSymbol', sps, ...
- 'DecimationFactor', 1);
+ 'DecimationFactor', 4, ...
+ 'DecimationOffset', mod(span+1, 4));
- r = filter([y; zeros(span * sps, 1)]);
- r = r(span * sps / 2 + 1 : end); % truncate filter transients
+ 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 conv
+ l = (length(coef.Numerator) - 1) / 2;
+ rr = [rs(l:end); rs(1:l-1)];
+
+ r = downsample(rr, 4, 3);
end
-function [xCD, kstart] = splitstepfourier(x, D, lambda, z, Tsamp,gamma,stepnum)
- %% Simulate chromatic dispersion.
+function xCDKerr = splitstepfourier(x, D, lambda, z, Tsamp, gamma)
+ %% Simulate chromatic dispersion and Kerr effect,
+ %% with attenuation and amplification.
%% Params:
%% - x: input waveform (pulse-shaped)
%% - D: dispersion coefficient (ps / (nm km))
%% - z: length of fibre (km)
%% - Tsamp: sampling time (s)
%% Output:
- %% - xCD: x after being dispersed. Energy of xCD is not normalized.
- %% - kstart: starting index of the discrete signal
+ %% - xCDKerr: x after being dispersed.
%% Convert everything to SI base units
c = 299792458; % m/s
z = z * 1e3; % m
gamma = gamma * 1e-3; % watt^-1 / m
- %stepnum = 1
- dz = z / stepnum;
+ dz = 1; % km
+ dz = dz * 1e3; % m
+ stepnum = z / dz;
- hhz = 1j * gamma * dz;
+ alpha = 0.2; % dB/km
+ alpha = alpha / 10 * log(10); % /km
+ alpha = alpha * 1e-3; % /m
+ hhz = 1j * gamma * dz;
- %% Time domain filter length, needed for compatibility with
- %% time-domain (convolution) chromatic dispersion.
- kmax = floor(abs(D) * lambda^2 * z / (2 * c * Tsamp^2));
- kstart = 1 - kmax;
- x = [zeros(kmax, 1); x]; % prepend zeros to allow for transients
+ xCDKerr = x .* exp(abs(x) .^ 2 .* hhz / 2 - alpha * dz / 4);
- xCD = x .* exp(abs(x) .^ 2 .* hhz / 2);
for i = 1 : stepnum
- xDFT = fft(xCD);
- n = length(xCD);
+ xDFT = fft(xCDKerr);
+ n = length(xCDKerr);
fs = 1 / Tsamp;
omega = (2*pi * fs / n * [(0 : floor((n-1)/2)), (-ceil((n-1)/2) : -1)]).';
dispDFT = xDFT .* exp(-1j * omega.^2 * D * lambda^2 * dz / (4 * pi * c));
+ xCDKerr = ifft(dispDFT);
- xCD = ifft(dispDFT);
-
- xCD = xCD .* exp(abs(xCD) .^ 2 .* hhz / 2);
- end
- xCD = xCD .* exp(-abs(xCD) .^ 2 .* hhz / 2);
-
- if ceil(kmax/2) > 0
- %% fix the order of the samples due to prepending zeros before FFT
- xCD = [xCD(ceil(kmax/2):end); xCD(1:ceil(kmax/2)-1)];
+ xCDKerr = xCDKerr .* exp(abs(xCDKerr) .^ 2 .* hhz - alpha * dz / 2);
+ if mod(i, 50) == 0
+ xCDKerr = xCDKerr * 10; % amplification
+ end
end
- %% pad zeros for compatibility with time-domain filter
- xCD = [zeros(floor((kmax-1)/2), 1); xCD; zeros(ceil((kmax+1)/2), 1)];
+ xCDKerr = xCDKerr .* exp(-abs(xCDKerr) .^ 2 .* hhz / 2 + alpha * dz / 4);
end
%% References
%% [1]: S.J. Savory, Digital filters for coherent optical receivers, 2008.
--- /dev/null
+function [xCDKerr, yCDKerr] = ...
+ ssf_pdm(x, y, D, lambda, z, Tsamp, gamma)
+ %% Simulate chromatic dispersion and Kerr effect,
+ %% with attenuation and amplification.
+ %% Params:
+ %% - x: input waveform (pulse-shaped)
+ %% - D: dispersion coefficient (ps / (nm km))
+ %% - lambda: wavelength (nm)
+ %% - z: length of fibre (km)
+ %% - Tsamp: sampling time (s)
+ %% Output:
+ %% - xCDKerr: x after being dispersed.
+
+ %% 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 = 1; % km
+ dz = dz * 1e3; % m
+ stepnum = z / dz;
+
+ alpha = 0.2; % dB/km
+ alpha = alpha / 10 * log(10); % /km
+ alpha = alpha * 1e-3; % /m
+
+ hhz = 1j * (8/9) * gamma * dz;
+
+ P = abs(x) .^ 2 + abs(y) .^ 2;
+
+ xCDKerr = x .* exp(P .* hhz / 2 - alpha * dz / 4);
+ yCDKerr = y .* exp(P .* hhz / 2 - alpha * dz / 4);
+
+ for i = 1 : stepnum
+ xDFT = fft(xCDKerr);
+ yDFT = fft(yCDKerr);
+ n = length(xCDKerr);
+ 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));
+
+ xCDKerr = ifft(xDFT .* dispDFT);
+ yCDKerr = ifft(yDFT .* dispDFT);
+
+ P = abs(xCDKerr) .^ 2 + abs(yCDKerr) .^ 2;
+ xCDKerr = xCDKerr .* exp(P .* hhz - alpha * dz / 2);
+ yCDKerr = yCDKerr .* exp(P .* hhz - alpha * dz / 2);
+ if mod(i, 50) == 0
+ xCDKerr = xCDKerr * 10; % amplification
+ yCDKerr = yCDKerr * 10;
+ end
+ end
+ P = abs(xCDKerr) .^ 2 + abs(yCDKerr) .^ 2;
+ xCDKerr = xCDKerr .* exp(-P .* hhz / 2 + alpha * dz / 4);
+ yCDKerr = yCDKerr .* exp(-P .* hhz / 2 + alpha * dz / 4);
+end
+%% References
+%% [1]: S.J. Savory, Digital filters for coherent optical receivers, 2008.
%% Output:
%% - x: pulse-shaped waveform
- filter = comm.RaisedCosineTransmitFilter...
+ txfilter = comm.RaisedCosineTransmitFilter...
('Shape', 'Square root', ...
'RolloffFactor', rolloff, ...
'FilterSpanInSymbols', span, ...
'OutputSamplesPerSymbol', sps, ...
'Gain', sqrt(sps)); % so that output has energy 1
- x = filter([modData; zeros(span, 1)]);
+
+ coef = coeffs(txfilter);
+ filter_fft = fft(coef.Numerator, length(modData) * sps);
+ modData_fft = fft(upsample(modData, sps));
+ x = ifft(modData_fft .* filter_fft.');
+ l = (length(coef.Numerator) - 1) / 2;
+ x = [x(l:end); x(1:l-1)];
end
projects / 4yp.git / commitdiff