From 5fae00773184080617ac022c07495d365975e0e1 Mon Sep 17 00:00:00 2001 From: Adrian Iain Lam Date: Tue, 12 Feb 2019 01:37:28 +0000 Subject: [PATCH] Working Kerr effect; PDM; speedups; removed unused files Kerr effect is now working, showing the expected curve shape as Tx power varies. Currently this uses a fixed step size of 1km, with attenuation, and amplification every 50km. Amplification will be removed later to simulate a true passive network, but now kept to verify results. Since non-linearity is now introduced, we change the Tx signal to having 8 samples/sec (previously 2). The Rx then downsamples at the Rx filter. "channel.m" is now added to simulate the overall channel with all effects considered. Currently this includes Kerr, CD and phase noise. PDM was then implemented, with "pdmchannel.m" being the PDM analogue of "channel.m", with corresponding new files for the PDM version of the CMA equalizer and the split-step solver. Simulations are now significantly faster by doing all filtering in the frequency domain (so multiplication instead of convolution). The behaviour is changed slightly: instead of zero initial conditions, we now take periodic (circular) boundary conditions. Simulation lengths were also changed to an integer power of 2 to have the most efficient FFT. Some initial testing files that are no longer needed are removed. --- CDCompensation.m | 17 ++++-- CD_AWGN.m | 30 ++++++---- adaptiveCMA.m | 9 +-- agrawalAppendixB.m | 58 ------------------ baseband.m | 69 ---------------------- channel.m | 123 ++++++++++++++++++++++++++++++++++++++ chromaticDispersion.m | 17 +++++- chromaticDispersion1Signal.m | 16 +++-- chromaticDispersionTest.m | 35 ----------- chromaticDispersion_FFT.m | 43 -------------- discretePSK_BER_SNR.m | 57 ------------------ formatFigure.m | 7 +++ kerr.m | 86 +++++++++++++-------------- kerr1Signal.m | 110 ---------------------------------- passband.m | 79 ------------------------- pdm_adaptiveCMA.m | 72 +++++++++++++++++++++++ pdmchannel.m | 137 +++++++++++++++++++++++++++++++++++++++++++ phasenoise1signal.m | 15 ++--- phasenoise_AWGN.m | 22 ++++--- rxFilter.m | 18 ++++-- splitstepfourier.m | 46 +++++++-------- ssf_pdm.m | 61 +++++++++++++++++++ txFilter.m | 10 +++- 23 files changed, 567 insertions(+), 570 deletions(-) delete mode 100644 agrawalAppendixB.m delete mode 100644 baseband.m create mode 100644 channel.m delete mode 100644 chromaticDispersionTest.m delete mode 100644 chromaticDispersion_FFT.m delete mode 100644 discretePSK_BER_SNR.m delete mode 100644 kerr1Signal.m delete mode 100644 passband.m create mode 100644 pdm_adaptiveCMA.m create mode 100644 pdmchannel.m create mode 100644 ssf_pdm.m diff --git a/CDCompensation.m b/CDCompensation.m index 1f3cc2c..9f01950 100644 --- a/CDCompensation.m +++ b/CDCompensation.m @@ -21,11 +21,20 @@ function [yCDComp, kstart] = CDCompensation(y, D, lambda, z, Tsamp) 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. diff --git a/CD_AWGN.m b/CD_AWGN.m index 4db7aad..241d22a 100644 --- a/CD_AWGN.m +++ b/CD_AWGN.m @@ -1,18 +1,18 @@ -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; @@ -36,26 +36,30 @@ x = txFilter(modData, rolloff, span, sps); %% 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; diff --git a/adaptiveCMA.m b/adaptiveCMA.m index 084eab4..41b1e48 100644 --- a/adaptiveCMA.m +++ b/adaptiveCMA.m @@ -1,4 +1,4 @@ -function adaptFilterOut = adaptiveCMA(rSampled) +function [adaptFilterOut, convergeIdx] = adaptiveCMA(rSampled) %% adaptive filter %% CMA taps = 19; % ODD taps @@ -17,6 +17,7 @@ function adaptFilterOut = adaptiveCMA(rSampled) adaptFilterOut = zeros(numSymbs, 1); converged = 0; convergeCount = 0; + convergeIdx = Inf; for it = 1:numSymbs if it <= (taps - 1) / 2; @@ -30,14 +31,14 @@ function adaptFilterOut = adaptiveCMA(rSampled) 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; diff --git a/agrawalAppendixB.m b/agrawalAppendixB.m deleted file mode 100644 index 295e508..0000000 --- a/agrawalAppendixB.m +++ /dev/null @@ -1,58 +0,0 @@ -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'); diff --git a/baseband.m b/baseband.m deleted file mode 100644 index d6f0f1b..0000000 --- a/baseband.m +++ /dev/null @@ -1,69 +0,0 @@ -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 diff --git a/channel.m b/channel.m new file mode 100644 index 0000000..260d21f --- /dev/null +++ b/channel.m @@ -0,0 +1,123 @@ +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; diff --git a/chromaticDispersion.m b/chromaticDispersion.m index 49a49a3..bc04e59 100644 --- a/chromaticDispersion.m +++ b/chromaticDispersion.m @@ -26,10 +26,21 @@ function [xCD, kstart] = chromaticDispersion(x, D, lambda, z, Tsamp) 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 diff --git a/chromaticDispersion1Signal.m b/chromaticDispersion1Signal.m index dd169ac..6441dd5 100644 --- a/chromaticDispersion1Signal.m +++ b/chromaticDispersion1Signal.m @@ -5,7 +5,7 @@ Rsym = 2.5e10; % symbol rate (sym/sec) 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; @@ -21,9 +21,9 @@ x = normalizeEnergy(x, numSymbs*sps, 1); %% 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); @@ -32,11 +32,16 @@ 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; @@ -57,7 +62,6 @@ end rSampled = rSampled .* exp(-1j * theta); - %%rAdaptEq = adaptiveCMA(rSampled); %{ %% Compare original signal and compensated signal diff --git a/chromaticDispersionTest.m b/chromaticDispersionTest.m deleted file mode 100644 index 513738f..0000000 --- a/chromaticDispersionTest.m +++ /dev/null @@ -1,35 +0,0 @@ -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'); diff --git a/chromaticDispersion_FFT.m b/chromaticDispersion_FFT.m deleted file mode 100644 index b091485..0000000 --- a/chromaticDispersion_FFT.m +++ /dev/null @@ -1,43 +0,0 @@ -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. diff --git a/discretePSK_BER_SNR.m b/discretePSK_BER_SNR.m deleted file mode 100644 index cdff01f..0000000 --- a/discretePSK_BER_SNR.m +++ /dev/null @@ -1,57 +0,0 @@ -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 diff --git a/formatFigure.m b/formatFigure.m index a09586d..c3d607f 100644 --- a/formatFigure.m +++ b/formatFigure.m @@ -4,3 +4,10 @@ set(gca, 'TickLabelInterpreter', 'Latex'); 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'); diff --git a/kerr.m b/kerr.m index c67bbb2..7e534dd 100644 --- a/kerr.m +++ b/kerr.m @@ -1,4 +1,4 @@ -numSymbs = 5e4; +numSymbs = 2^16; M = 4; Rsym = 2.5e10; % symbol rate (sym/sec) @@ -6,40 +6,46 @@ Tsym = 1 / Rsym; % symbol period (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); @@ -50,58 +56,52 @@ for i = 1:plotlen %% 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; diff --git a/kerr1Signal.m b/kerr1Signal.m deleted file mode 100644 index f4ea292..0000000 --- a/kerr1Signal.m +++ /dev/null @@ -1,110 +0,0 @@ -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; diff --git a/passband.m b/passband.m deleted file mode 100644 index 9e88af6..0000000 --- a/passband.m +++ /dev/null @@ -1,79 +0,0 @@ -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 diff --git a/pdm_adaptiveCMA.m b/pdm_adaptiveCMA.m new file mode 100644 index 0000000..7e6c9eb --- /dev/null +++ b/pdm_adaptiveCMA.m @@ -0,0 +1,72 @@ +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 diff --git a/pdmchannel.m b/pdmchannel.m new file mode 100644 index 0000000..997d646 --- /dev/null +++ b/pdmchannel.m @@ -0,0 +1,137 @@ +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; diff --git a/phasenoise1signal.m b/phasenoise1signal.m index f819e65..e2240b9 100644 --- a/phasenoise1signal.m +++ b/phasenoise1signal.m @@ -4,7 +4,7 @@ M = 4; 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; @@ -39,18 +39,19 @@ linewidthLO = 10e6; % Hz 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'); @@ -60,9 +61,9 @@ 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; diff --git a/phasenoise_AWGN.m b/phasenoise_AWGN.m index f055bcc..621a33d 100644 --- a/phasenoise_AWGN.m +++ b/phasenoise_AWGN.m @@ -1,11 +1,11 @@ -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; @@ -48,15 +48,20 @@ linewidthTx = 0; % Hz 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); @@ -68,9 +73,12 @@ for it = 1 : iterations 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); @@ -99,7 +107,7 @@ for it = 1 : iterations if EbN0_db(i) == 8 && it == 1 figure(1234); - plot(repelem(-phiestsPSK, sps)); + plot(repelem(-phiestsPSK, 8)); hold on; plot(pTxLoPSK); legend('estimate', 'actual'); diff --git a/rxFilter.m b/rxFilter.m index eb25a32..e8e572a 100644 --- a/rxFilter.m +++ b/rxFilter.m @@ -8,13 +8,23 @@ function r = rxFilter(y, rolloff, span, sps) %% 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 diff --git a/splitstepfourier.m b/splitstepfourier.m index f71828d..f0ec5fb 100644 --- a/splitstepfourier.m +++ b/splitstepfourier.m @@ -1,5 +1,6 @@ -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)) @@ -7,8 +8,7 @@ function [xCD, kstart] = splitstepfourier(x, D, lambda, z, Tsamp,gamma,stepnum) %% - 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 @@ -17,40 +17,34 @@ function [xCD, kstart] = splitstepfourier(x, D, lambda, z, Tsamp,gamma,stepnum) 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. diff --git a/ssf_pdm.m b/ssf_pdm.m new file mode 100644 index 0000000..55511fe --- /dev/null +++ b/ssf_pdm.m @@ -0,0 +1,61 @@ +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. diff --git a/txFilter.m b/txFilter.m index 507ffe1..10c387d 100644 --- a/txFilter.m +++ b/txFilter.m @@ -8,11 +8,17 @@ function x = txFilter(modData, rolloff, span, sps) %% 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 -- 2.7.4