From 1ccbdf85c430d80efbc6884c433f595149b36251 Mon Sep 17 00:00:00 2001 From: Sourcery AI <> Date: Tue, 1 Nov 2022 19:34:42 +0000 Subject: [PATCH] 'Refactored by Sourcery' --- .../CommpyTest_BERvsEbN0_wPulseShaping.py | 66 ++++++----- ...e sinais \303\263pticos e ru\303\255do.py" | 4 +- ...\247\303\243o de sinais \303\263pticos.py" | 4 +- jupyter notebooks/Sistemas coerentes.py | 28 ++--- .../Split-Step Fourier Method (SSFM).py | 53 +++++---- .../.ipynb_checkpoints/dsp-checkpoint.py | 2 +- jupyter notebooks/utils/dsp.py | 89 +++++++-------- jupyter notebooks/utils/models.py | 104 +++++++++--------- jupyter notebooks/utils/tx.py | 43 ++++---- 9 files changed, 188 insertions(+), 205 deletions(-) diff --git a/jupyter notebooks/CommpyTest_BERvsEbN0_wPulseShaping.py b/jupyter notebooks/CommpyTest_BERvsEbN0_wPulseShaping.py index 3fda061..cec58a3 100644 --- a/jupyter notebooks/CommpyTest_BERvsEbN0_wPulseShaping.py +++ b/jupyter notebooks/CommpyTest_BERvsEbN0_wPulseShaping.py @@ -68,10 +68,10 @@ def GrayMapping(M, constType): L = int(np.sqrt(M)-1) bitsSymb = int(np.log2(M)) - + code = GrayCode(bitsSymb) a = list(code.generate_gray()) - + if constType == 'qam': PAM = np.arange(-L, L+1, 2) PAM = np.array([PAM]) @@ -79,27 +79,27 @@ def GrayMapping(M, constType): # generate complex square M-QAM constellation const = repmat(PAM, L+1, 1) + 1j*repmat(np.flip(PAM.T,0), 1, L+1) const = const.T - + for ind in np.arange(1,L+1,2): const[ind] = np.flip(const[ind],0) - + elif constType == 'psk': pskPhases = np.arange(0,2*np.pi,2*np.pi/M) - + # generate complex M-PSK constellation const = np.exp(1j*pskPhases) - + const = const.reshape(M,1) const_ = np.zeros((M,2),dtype=complex) - - for ind in range(0,M): + + for ind in range(M): const_[ind,0] = const[ind,0] # complex constellation symbol const_[ind,1] = int(a[ind],2) # mapped bit sequence (as integer decimal) - + # sort complex symbols column according to their mapped bit sequence (as integer decimal) const = const_[const_[:,1].real.argsort()] - + return const @@ -121,12 +121,10 @@ def modulateGray(bits, M, constType): def demodulateGray(symb, M ,constType): const = GrayMapping(M, constType) - - hdDecision_vec = np.vectorize(hdDecision, excluded = [1]) - index_list = hdDecision_vec(symb, const[:,0]) - demod_bits = dec2bitarray(index_list, int(np.log2(M))) - - return demod_bits + + hdDecision_vec = np.vectorize(hdDecision, excluded = [1]) + index_list = hdDecision_vec(symb, const[:,0]) + return dec2bitarray(index_list, int(np.log2(M))) def pulseShape(pulseType, SpS=2, N=1024, alpha=0.1, Ts=1): @@ -158,7 +156,7 @@ def pulseShape(pulseType, SpS=2, N=1024, alpha=0.1, Ts=1): alpha = 0.05 # Rolloff do filtro RRC N = 1024 # Número de coeficientes do filtro RRC EbN0dB = 20 - + # generate random bits bitsTx = np.random.randint(2, size=3*2**12) @@ -176,7 +174,7 @@ def pulseShape(pulseType, SpS=2, N=1024, alpha=0.1, Ts=1): # pulse shaping #pulseFilter = pulseShape('rrc', SpS, N, alpha, Ts) tindex, rrcFilter = rrcosfilter(N, alpha, Ts, Fa) -sigTx = filterNoDelay(pulseFilter, symbolsUp) +sigTx = filterNoDelay(pulseFilter, symbolsUp) sigTx = sigTx/np.sqrt(signal_power(sigTx)) # plot eye diagrams @@ -211,7 +209,7 @@ def pulseShape(pulseType, SpS=2, N=1024, alpha=0.1, Ts=1): # + # matched filter -sigRx = filterNoDelay(rrcFilter, sigRx) +sigRx = filterNoDelay(rrcFilter, sigRx) sigRx = sigRx/SpS # plot eye diagrams @@ -240,7 +238,7 @@ def pulseShape(pulseType, SpS=2, N=1024, alpha=0.1, Ts=1): BERtheory = theoryBER(M, EbN0dB,'qam') # print results -print('EbN0: %3.2f dB, EbN0_est: %3.2f dB,\nBERtheory: %3.1e, BER: %3.1e ' %(EbN0dB, EbN0dB_est, BERtheory, BER)) +print('EbN0: %3.2f dB, EbN0_est: %3.2f dB,\nBERtheory: %3.1e, BER: %3.1e ' %(EbN0dB, EbN0dB_est, BERtheory, BER)) print('Total of bits counted: ', ERR.size) # plot constellations @@ -268,10 +266,11 @@ def pulseShape(pulseType, SpS=2, N=1024, alpha=0.1, Ts=1): BER = np.zeros(EbN0dB_.shape) EbN0dB_est = np.zeros(EbN0dB_.shape) +discard = 100 for indSNR in range(EbN0dB_.size): EbN0dB = EbN0dB_[indSNR] - + # generate random bits bitsTx = np.random.randint(2, size=3*2**18) @@ -279,18 +278,18 @@ def pulseShape(pulseType, SpS=2, N=1024, alpha=0.1, Ts=1): mod = QAMModem(m=M) symbTx = mod.modulate(bitsTx) Es = mod.Es; - + # Normalize symbols energy to 1 symbTx = symbTx/np.sqrt(Es) - + # Aumenta a taxa de amostragem do sinal para SpS amostras/símbolo symbolsUp = upsample(symbTx, SpS) - + # filtro formatador de pulso tindex, rrcFilter = rrcosfilter(N, alpha, Ts, Fa) - symbolsUp = filterNoDelay(rrcFilter, symbolsUp) + symbolsUp = filterNoDelay(rrcFilter, symbolsUp) symbolsUp = symbolsUp/np.sqrt(signal_power(symbolsUp)) - + # AWGN channel snrdB = EbN0dB + 10*np.log10(np.log2(M)) noiseVar = 1/(10**(snrdB/10)) @@ -300,24 +299,23 @@ def pulseShape(pulseType, SpS=2, N=1024, alpha=0.1, Ts=1): noise = 1/np.sqrt(2)*noise symbRx = symbolsUp + noise - + # filtro casado - symbRx = filterNoDelay(rrcFilter, symbRx) + symbRx = filterNoDelay(rrcFilter, symbRx) symbRx = symbRx/SpS - + # Decimação para uma amostra por símbolo symbRx = symbRx[0::SpS] - + # Demodulate received symbols - bitsRx = mod.demodulate(np.sqrt(Es)*symbRx, demod_type = 'hard') - discard = 100 + bitsRx = mod.demodulate(np.sqrt(Es)*symbRx, demod_type = 'hard') numBits = bitsTx.size; - + # BER calculation, EbN0 estimation ERR = np.logical_xor(bitsRx[discard:numBits-discard], bitsTx[discard:numBits-discard]) BER[indSNR] = ERR.sum()/ERR.size EbN0dB_est[indSNR] = 10*np.log10(1/(signal_power(symbRx-symbTx)*np.log2(M))) - + print('EbN0: %3.2f dB, EbN0_est: %3.2f dB, BER: %3.1e ' %(EbN0dB, EbN0dB_est[indSNR], BER[indSNR])) diff --git "a/jupyter notebooks/Detec\303\247\303\243o de sinais \303\263pticos e ru\303\255do.py" "b/jupyter notebooks/Detec\303\247\303\243o de sinais \303\263pticos e ru\303\255do.py" index 555e529..39be16d 100644 --- "a/jupyter notebooks/Detec\303\247\303\243o de sinais \303\263pticos e ru\303\255do.py" +++ "b/jupyter notebooks/Detec\303\247\303\243o de sinais \303\263pticos e ru\303\255do.py" @@ -145,11 +145,11 @@ σ2_s = 2*q*(Ip + Id)*B # variância μ = 0 # média -σ = np.sqrt(σ2_s) +σ = np.sqrt(σ2_s) Is = normal(μ, σ, Namostras) # plotas as primeiras 1000 amostras -plt.plot(Is[0:1000],linewidth = 0.8); +plt.plot(Is[:1000], linewidth = 0.8); plt.xlim(0,1000) plt.ylabel('Is') plt.xlabel('amostra') diff --git "a/jupyter notebooks/Gera\303\247\303\243o de sinais \303\263pticos.py" "b/jupyter notebooks/Gera\303\247\303\243o de sinais \303\263pticos.py" index 4260272..7b233ae 100644 --- "a/jupyter notebooks/Gera\303\247\303\243o de sinais \303\263pticos.py" +++ "b/jupyter notebooks/Gera\303\247\303\243o de sinais \303\263pticos.py" @@ -301,9 +301,7 @@ def mzm(Ai, Vπ, u, Vb): :return Ao: amplitude da portadora modulada """ π = np.pi - Ao = Ai*np.cos(0.5/Vπ*(u+Vb)*π) - - return Ao + return Ai*np.cos(0.5/Vπ*(u+Vb)*π) # ### Transmitindo informação na intensidade (potência) da portadora óptica diff --git a/jupyter notebooks/Sistemas coerentes.py b/jupyter notebooks/Sistemas coerentes.py index f2a245f..e5426b0 100644 --- a/jupyter notebooks/Sistemas coerentes.py +++ b/jupyter notebooks/Sistemas coerentes.py @@ -298,17 +298,15 @@ def hybrid_2x4_90(E1, E2): M = Matrix([[C_3dB, zeros(2)], [zeros(2), C_3dB]]) - + U = Matrix([[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, j]]) - + Ei = Matrix([[E1],[0],[0],[E2]]) # vetor 4x1 - - Eo = M*U*M*Ei - - return Eo + + return M*U*M*Ei # fotodetector balanceado def bpd(E1, E2, R=1): @@ -406,19 +404,17 @@ def hybrid_2x4_90deg(E1, E2): :return: hybrid outputs ''' assert E1.size == E2.size, 'E1 and E2 need to have the same size' - + # optical hybrid transfer matrix T = np.array([[ 1/2, 1j/2, 1j/2, -1/2], [ 1j/2, -1/2, 1/2, 1j/2], [ 1j/2, 1/2, -1j/2, -1/2], [-1/2, 1j/2, -1/2, 1j/2]]) - + Ei = np.array([E1, np.zeros((E1.size,)), np.zeros((E1.size,)), E2]) - - Eo = T@Ei - - return Eo + + return T@Ei def coherentReceiver(Es, Elo, Rd=1): ''' @@ -447,12 +443,12 @@ def coherentReceiver(Es, Elo, Rd=1): def phaseNoise(lw, Nsamples, Ts): - σ2 = 2*np.pi*lw*Ts + σ2 = 2*np.pi*lw*Ts phi = np.zeros(Nsamples) - - for ind in range(0, Nsamples-1): + + for ind in range(Nsamples-1): phi[ind+1] = phi[ind] + normal(0, np.sqrt(σ2)) - + return phi diff --git a/jupyter notebooks/Split-Step Fourier Method (SSFM).py b/jupyter notebooks/Split-Step Fourier Method (SSFM).py index 842e76f..c176a02 100644 --- a/jupyter notebooks/Split-Step Fourier Method (SSFM).py +++ b/jupyter notebooks/Split-Step Fourier Method (SSFM).py @@ -40,29 +40,29 @@ def theoryBER(M, EbN0, constType): return Pb def CDcompensation(Ex, Ey, Ltotal, D, Fc, Fs): - + c = 299792458 # speed of light (vacuum) c_kms = c/1e3 λ = c_kms/Fc β2 = -(D*λ**2)/(2*np.pi*c_kms) γ = gamma - + l = len(Ex) - + Nfft = l; ω = 2*np.pi*Fs*fftfreq(Nfft) - + Ex = fft(Ex) #Pol. X Ey = fft(Ey) #Pol. Y - + # CD compensation Ex = Ex*np.exp(-1j*(β2/2)*(ω**2)*(Ltotal)) Ey = Ey*np.exp(-1j*(β2/2)*(ω**2)*(Ltotal)) - + Ex = ifft(Ex) Ey = ifft(Ey) - + return Ex, Ey @@ -70,29 +70,28 @@ def manakovSSF(Ei, hz, Lspan, Ltotal, alpha, gamma, D, Fc, Fs): Ex = Ei[:,0] Ey = Ei[:,1] - + c = 299792458 # speed of light (vacuum) c_kms = c/1e3 λ = c_kms/Fc α = alpha/(10*np.log10(np.exp(1))) β2 = -(D*λ**2)/(2*np.pi*c_kms) γ = gamma - + Nfft = len(Ex) ω = 2*np.pi*Fs*fftfreq(Nfft) - + Nspans = int(np.floor(Ltotal/Lspan)) Nsteps = int(np.floor(Lspan/hz)) Ex = fft(Ex) #Pol. X Ey = fft(Ey) #Pol. Y - + linOperator = np.exp(-(α/2)*(hz/2) + 1j*(β2/2)*(ω**2)*(hz/2)) - - for spanN in tqdm(range(1, Nspans+1)): - for stepN in range(1, Nsteps+1): - + + for _ in tqdm(range(1, Nspans+1)): + for _ in range(1, Nsteps+1): # First linear step (frequency domain) Ex = Ex*linOperator Ey = Ey*linOperator @@ -100,22 +99,22 @@ def manakovSSF(Ei, hz, Lspan, Ltotal, alpha, gamma, D, Fc, Fs): # Nonlinear step (time domain) Ex = ifft(Ex); Ey = ifft(Ey); - + Ex = Ex*np.exp(1j*γ*8/9*(Ex*np.conj(Ex) + Ey*np.conj(Ey))*hz) Ey = Ey*np.exp(1j*γ*8/9*(Ex*np.conj(Ex) + Ey*np.conj(Ey))*hz) - + # Second linear step (frequency domain) Ex = fft(Ex); Ey = fft(Ey); Ex = Ex*linOperator Ey = Ey*linOperator - + Ex = Ex*np.exp(α*Nsteps*hz) Ey = Ey*np.exp(α*Nsteps*hz) - + Ex = ifft(Ex) Ey = ifft(Ey) - + return Ex, Ey def edfa_lin(signal, gain, nf, fc, fs): @@ -207,7 +206,7 @@ def filterNoDelay(h, x): # generate random bits np.random.seed(33) -bits_x = np.random.randint(2, size=3*2**14) +bits_x = np.random.randint(2, size=3*2**14) np.random.seed(23) bits_y = np.random.randint(2, size=3*2**14) @@ -228,7 +227,7 @@ def filterNoDelay(h, x): # pulse shaping tindex, rrcFilter = rrcosfilter(N, alphaRRC, Ts, Fa) -sig_x = filterNoDelay(rrcFilter, symbolsUp_x) +sig_x = filterNoDelay(rrcFilter, symbolsUp_x) sig_y = filterNoDelay(rrcFilter, symbolsUp_y) sig_x = np.sqrt(P0/2)*sig_x/np.sqrt(signal_power(sig_x)) @@ -259,7 +258,7 @@ def filterNoDelay(h, x): sigRx = sigRx/np.sqrt(signal_power(sigRx)) # matched filter -sigRx = filterNoDelay(rrcFilter, sigRx) +sigRx = filterNoDelay(rrcFilter, sigRx) sigRx = sigRx/SpS # downsampling to one sample per symbol @@ -290,7 +289,7 @@ def filterNoDelay(h, x): BERtheory = theoryBER(M, EbN0dB,'qam') # print results -print('EbN0: %3.2f dB, EbN0_est: %3.2f dB,\nBERtheory: %3.1e, BER: %3.1e ' %(EbN0dB, EbN0dB_est, BERtheory, BER)) +print('EbN0: %3.2f dB, EbN0_est: %3.2f dB,\nBERtheory: %3.1e, BER: %3.1e ' %(EbN0dB, EbN0dB_est, BERtheory, BER)) print('Total of bits counted: ', ERR.size) # plot constellations @@ -306,7 +305,7 @@ def filterNoDelay(h, x): # + from scipy.stats.kde import gaussian_kde -y = (sigRx[0:30000]).real +y = sigRx[:30000].real x = np.arange(0,y.size,1) % 24 k = gaussian_kde(np.vstack([x, y])) @@ -319,8 +318,8 @@ def filterNoDelay(h, x): # + from scipy.stats.kde import gaussian_kde -y = (symbRx[0:30000]).real -x = (symbRx[0:30000]).imag +y = symbRx[:30000].real +x = symbRx[:30000].imag k = gaussian_kde(np.vstack([x, y])) k.set_bandwidth(bw_method=k.factor/2) diff --git a/jupyter notebooks/utils/.ipynb_checkpoints/dsp-checkpoint.py b/jupyter notebooks/utils/.ipynb_checkpoints/dsp-checkpoint.py index fe76609..25ec513 100644 --- a/jupyter notebooks/utils/.ipynb_checkpoints/dsp-checkpoint.py +++ b/jupyter notebooks/utils/.ipynb_checkpoints/dsp-checkpoint.py @@ -44,7 +44,7 @@ def pulseShape(pulseType, SpS=2, N=1024, alpha=0.1, Ts=1): def eyediagram(sig, Nsamples, SpS, n=3): - y = sig[0:Nsamples].real + y = sig[:Nsamples].real x = np.arange(0,y.size,1) % (n*SpS) k = gaussian_kde(np.vstack([x, y])) diff --git a/jupyter notebooks/utils/dsp.py b/jupyter notebooks/utils/dsp.py index 327e1f2..6bda216 100644 --- a/jupyter notebooks/utils/dsp.py +++ b/jupyter notebooks/utils/dsp.py @@ -66,19 +66,19 @@ def eyediagram(sig, Nsamples, SpS, n=3, ptype='fast', plotlabel=None): if np.iscomplex(sig).any(): d = 1 - plotlabel_ = plotlabel+' [real]' + plotlabel_ = f'{plotlabel} [real]' else: d = 0 plotlabel_ = plotlabel - - for ind in range(0, d+1): + + for ind in range(d+1): if ind == 0: - y = sig[0:Nsamples].real - x = np.arange(0,y.size,1) % (n*SpS) + y = sig[:Nsamples].real + x = np.arange(0,y.size,1) % (n*SpS) else: - y = sig[0:Nsamples].imag - plotlabel_ = plotlabel+' [imag]' - + y = sig[:Nsamples].imag + plotlabel_ = f'{plotlabel} [imag]' + plt.figure(); if ptype == 'fancy': k = gaussian_kde(np.vstack([x, y])) @@ -91,16 +91,16 @@ def eyediagram(sig, Nsamples, SpS, n=3, ptype='fast', plotlabel=None): elif ptype == 'fast': y[x == n*SpS] = np.nan; y[x == 0] = np.nan; - + plt.plot(x/SpS, y, color='blue', alpha=0.8, label=plotlabel_); plt.xlim(min(x/SpS), max(x/SpS)) plt.xlabel('symbol period (Ts)') plt.ylabel('amplitude') plt.title('eye diagram') - + if plotlabel != None: plt.legend(loc='upper left') - + plt.grid() plt.show(); return None @@ -111,22 +111,22 @@ def sincInterp(x, fa): Ta_sinc = 1/fa_sinc Ta = 1/fa t = np.arange(0, x.size*32)*Ta_sinc - - plt.figure() + + plt.figure() y = upsample(x,32) y[y==0] = np.nan plt.plot(t,y.real,'ko', label='x[k]') - + x_sum = 0 - for k in range(0, x.size): + for k in range(x.size): xk_interp = x[k]*np.sinc((t-k*Ta)/Ta) x_sum += xk_interp plt.plot(t, xk_interp) - + plt.legend(loc="upper right") plt.xlim(min(t), max(t)) plt.grid() - + return x_sum, t def lowPassFIR(fc, fa, N, typeF = 'rect'): @@ -164,34 +164,28 @@ def edc(Ei, L, D, Fc, Fs): :return Eo: CD compensated signal """ - Eo = linFiberCh(Ei, L, 0, -D, Fc, Fs) - - return Eo + return linFiberCh(Ei, L, 0, -D, Fc, Fs) -def cpr(Ei, N, constSymb, symbTx): +def cpr(Ei, N, constSymb, symbTx): """ Carrier phase recovery (CPR) """ - ϕ = np.zeros(Ei.shape) + ϕ = np.zeros(Ei.shape) θ = np.zeros(Ei.shape) - - for k in range(0,len(Ei)): + + for k in range(len(Ei)): decided = np.argmin(np.abs(Ei[k]*np.exp(1j*θ[k-1]) - constSymb)) # find closest constellation symbol - + if k % 50 == 0: ϕ[k] = np.angle(symbTx[k]/(Ei[k])) # phase estimation with pilot symbol else: ϕ[k] = np.angle(constSymb[decided]/(Ei[k])) # phase estimation after symbol decision - - if k > N: - θ[k] = np.mean(ϕ[k-N:k]) # moving average filter - else: - θ[k] = np.angle(symbTx[k]/(Ei[k])) - + + θ[k] = np.mean(ϕ[k-N:k]) if k > N else np.angle(symbTx[k]/(Ei[k])) Eo = Ei*np.exp(1j*θ) # compensate phase rotation - + return Eo, ϕ, θ def fourthPowerFOE(Ei, Ts, plotSpec=False): @@ -218,7 +212,7 @@ def fourthPowerFOE(Ei, Ts, plotSpec=False): return f[indFO]/4 -def dbp(Ei, Fs, Ltotal, Lspan, hz=0.5, alpha=0.2, gamma=1.3, D=16, Fc=193.1e12): +def dbp(Ei, Fs, Ltotal, Lspan, hz=0.5, alpha=0.2, gamma=1.3, D=16, Fc=193.1e12): """ Digital backpropagation (symmetric, single-pol.) @@ -240,35 +234,34 @@ def dbp(Ei, Fs, Ltotal, Lspan, hz=0.5, alpha=0.2, gamma=1.3, D=16, Fc=193.1e12): α = -alpha/(10*np.log10(np.exp(1))) β2 = (D*λ**2)/(2*np.pi*c_kms) γ = -gamma - + Nfft = len(Ei) ω = 2*np.pi*Fs*fftfreq(Nfft) - + Nspans = int(np.floor(Ltotal/Lspan)) Nsteps = int(np.floor(Lspan/hz)) - - Ech = Ei.reshape(len(Ei),) + + Ech = Ei.reshape(len(Ei),) Ech = fft(Ech) #single-polarization field - + linOperator = np.exp(-(α/2)*(hz/2) + 1j*(β2/2)*(ω**2)*(hz/2)) - - for spanN in tqdm(range(0, Nspans)): - + + for _ in tqdm(range(Nspans)): Ech = Ech*np.exp((α/2)*Nsteps*hz) - - for stepN in range(0, Nsteps): + + for _ in range(Nsteps): # First linear step (frequency domain) Ech = Ech*linOperator - + # Nonlinear step (time domain) Ech = ifft(Ech) Ech = Ech*np.exp(1j*γ*(Ech*np.conj(Ech))*hz) - + # Second linear step (frequency domain) - Ech = fft(Ech) + Ech = fft(Ech) Ech = Ech*linOperator - + Ech = ifft(Ech) - + return Ech.reshape(len(Ech),) diff --git a/jupyter notebooks/utils/models.py b/jupyter notebooks/utils/models.py index e703af3..b97b9d0 100644 --- a/jupyter notebooks/utils/models.py +++ b/jupyter notebooks/utils/models.py @@ -17,9 +17,7 @@ def mzm(Ai, Vπ, u, Vb): :return Ao: output optical signal """ π = np.pi - Ao = Ai*np.cos(0.5/Vπ*(u+Vb)*π) - - return Ao + return Ai*np.cos(0.5/Vπ*(u+Vb)*π) def iqm(Ai, u, Vπ, VbI, VbQ): """ @@ -33,9 +31,9 @@ def iqm(Ai, u, Vπ, VbI, VbQ): :return Ao: output optical signal """ - Ao = mzm(Ai/np.sqrt(2), Vπ, u.real, VbI) + 1j*mzm(Ai/np.sqrt(2), Vπ, u.imag, VbQ) - - return Ao + return mzm(Ai / np.sqrt(2), Vπ, u.real, VbI) + 1j * mzm( + Ai / np.sqrt(2), Vπ, u.imag, VbQ + ) def linFiberCh(Ei, L, alpha, D, Fc, Fs): """ @@ -103,19 +101,17 @@ def hybrid_2x4_90deg(E1, E2): :return: hybrid outputs """ assert E1.size == E2.size, 'E1 and E2 need to have the same size' - + # optical hybrid transfer matrix T = np.array([[ 1/2, 1j/2, 1j/2, -1/2], [ 1j/2, -1/2, 1/2, 1j/2], [ 1j/2, 1/2, -1j/2, -1/2], [-1/2, 1j/2, -1/2, 1j/2]]) - + Ei = np.array([E1, np.zeros((E1.size,)), np.zeros((E1.size,)), E2]) - - Eo = T@Ei - - return Eo + + return T@Ei def coherentReceiver(Es, Elo, Rd=1): """ @@ -163,7 +159,7 @@ def edfa(Ei, Fs, G=20, NF=4.5, Fc=193.1e12): noise = normal(0, np.sqrt(p_noise), Ei.shape) + 1j*normal(0, np.sqrt(p_noise), Ei.shape) return Ei*np.sqrt(G_lin) + noise -def ssfm(Ei, Fs, paramCh): +def ssfm(Ei, Fs, paramCh): """ Split-step Fourier method (symmetric, single-pol.) @@ -194,14 +190,14 @@ def ssfm(Ei, Fs, paramCh): paramCh.amp = getattr(paramCh, 'amp', 'edfa') paramCh.NF = getattr(paramCh, 'NF', 4.5) - Ltotal = paramCh.Ltotal + Ltotal = paramCh.Ltotal Lspan = paramCh.Lspan hz = paramCh.hz - alpha = paramCh.alpha - D = paramCh.D - gamma = paramCh.gamma - Fc = paramCh.Fc - amp = paramCh.amp + alpha = paramCh.alpha + D = paramCh.D + gamma = paramCh.gamma + Fc = paramCh.Fc + amp = paramCh.amp NF = paramCh.NF # channel parameters @@ -214,20 +210,20 @@ def ssfm(Ei, Fs, paramCh): # generate frequency axis Nfft = len(Ei) ω = 2*np.pi*Fs*fftfreq(Nfft) - + Nspans = int(np.floor(Ltotal/Lspan)) Nsteps = int(np.floor(Lspan/hz)) - + Ech = Ei.reshape(len(Ei),) # define linear operator linOperator = np.exp(-(α/2)*(hz/2) + 1j*(β2/2)*(ω**2)*(hz/2)) - - for spanN in tqdm(range(1, Nspans+1)): + + for _ in tqdm(range(1, Nspans+1)): Ech = fft(Ech) #single-polarization field - + # fiber propagation step - for stepN in range(1, Nsteps+1): + for _ in range(1, Nsteps+1): # First linear step (frequency domain) Ech = Ech*linOperator @@ -236,22 +232,22 @@ def ssfm(Ei, Fs, paramCh): Ech = Ech*np.exp(1j*γ*(Ech*np.conj(Ech))*hz) # Second linear step (frequency domain) - Ech = fft(Ech) + Ech = fft(Ech) Ech = Ech*linOperator # amplification step Ech = ifft(Ech) - if amp =='edfa': + if amp == 'edfa': Ech = edfa(Ech, Fs, alpha*Lspan, NF, Fc) - elif amp =='ideal': + elif amp == 'ideal': Ech = Ech*np.exp(α/2*Nsteps*hz) - elif amp == None: + elif amp is None: Ech = Ech*np.exp(0); - + return Ech.reshape(len(Ech),), paramCh -def manakov_ssf(Ei, Fs, paramCh): +def manakov_ssf(Ei, Fs, paramCh): """ Manakov model split-step Fourier (symmetric, dual-pol.) @@ -282,14 +278,14 @@ def manakov_ssf(Ei, Fs, paramCh): paramCh.amp = getattr(paramCh, 'amp', 'edfa') paramCh.NF = getattr(paramCh, 'NF', 4.5) - Ltotal = paramCh.Ltotal + Ltotal = paramCh.Ltotal Lspan = paramCh.Lspan hz = paramCh.hz - alpha = paramCh.alpha - D = paramCh.D - gamma = paramCh.gamma - Fc = paramCh.Fc - amp = paramCh.amp + alpha = paramCh.alpha + D = paramCh.D + gamma = paramCh.gamma + Fc = paramCh.Fc + amp = paramCh.amp NF = paramCh.NF # channel parameters @@ -302,22 +298,22 @@ def manakov_ssf(Ei, Fs, paramCh): # generate frequency axis Nfft = len(Ei) ω = 2*np.pi*Fs*fftfreq(Nfft) - + Nspans = int(np.floor(Ltotal/Lspan)) Nsteps = int(np.floor(Lspan/hz)) - + Ech_x = Ei[:,0].reshape(len(Ei),) Ech_y = Ei[:,1].reshape(len(Ei),) # define linear operator linOperator = np.exp(-(α/2)*(hz/2) + 1j*(β2/2)*(ω**2)*(hz/2)) - - for spanN in tqdm(range(1, Nspans+1)): + + for _ in tqdm(range(1, Nspans+1)): Ech_x = fft(Ech_x) #polarization x field Ech_y = fft(Ech_y) #polarization y field - + # fiber propagation step - for stepN in range(1, Nsteps+1): + for _ in range(1, Nsteps+1): # First linear step (frequency domain) Ech_x = Ech_x*linOperator Ech_y = Ech_y*linOperator @@ -331,35 +327,35 @@ def manakov_ssf(Ei, Fs, paramCh): # Second linear step (frequency domain) Ech_x = fft(Ech_x) Ech_y = fft(Ech_y) - + Ech_x = Ech_x*linOperator Ech_y = Ech_y*linOperator - + # amplification step Ech_x = ifft(Ech_x) Ech_y = ifft(Ech_y) - - if amp =='edfa': + + if amp == 'edfa': Ech_x = edfa(Ech_x, Fs, alpha*Lspan, NF, Fc) Ech_y = edfa(Ech_y, Fs, alpha*Lspan, NF, Fc) - elif amp =='ideal': + elif amp == 'ideal': Ech_x = Ech_x*np.exp(α/2*Nsteps*hz) Ech_y = Ech_y*np.exp(α/2*Nsteps*hz) - elif amp == None: + elif amp is None: Ech_x = Ech_x*np.exp(0); Ech_y = Ech_y*np.exp(0); Ech = np.array([Ech_x.reshape(len(Ei),), Ech_y.reshape(len(Ei),)]).T - + return Ech, paramCh def phaseNoise(lw, Nsamples, Ts): - σ2 = 2*np.pi*lw*Ts + σ2 = 2*np.pi*lw*Ts phi = np.zeros(Nsamples) - - for ind in range(0, Nsamples-1): + + for ind in range(Nsamples-1): phi[ind+1] = phi[ind] + normal(0, np.sqrt(σ2)) - + return phi diff --git a/jupyter notebooks/utils/tx.py b/jupyter notebooks/utils/tx.py index 49ce2c5..0b5ce3c 100644 --- a/jupyter notebooks/utils/tx.py +++ b/jupyter notebooks/utils/tx.py @@ -37,38 +37,41 @@ def simpleWDMTx(param): param.Fc = getattr(param, 'Fc', 193.1e12) param.freqSpac = getattr(param, 'freqSpac', 50e9) param.Nmodes = getattr(param, 'Nmodes', 1) - + # transmitter parameters Ts = 1/param.Rs # symbol period [s] Fa = 1/(Ts/param.SpS) # sampling frequency [samples/s] Ta = 1/Fa # sampling period [s] - + # central frequencies of the WDM channels freqGrid = np.arange(-np.floor(param.Nch/2), np.floor(param.Nch/2)+1,1)*param.freqSpac - + if (param.Nch % 2) == 0: freqGrid += param.freqSpac/2 - + # IQM parameters Ai = 1 Vπ = 2 Vb = -Vπ Pch = 10**(param.Pch_dBm/10)*1e-3 # optical signal power per WDM channel - + π = np.pi - + t = np.arange(0, int(((param.Nbits)/np.log2(param.M))*param.SpS)) - + # allocate array sigTxWDM = np.zeros((len(t), param.Nmodes), dtype='complex') - symbTxWDM = np.zeros((int(len(t)/param.SpS), param.Nmodes, param.Nch), dtype='complex') - + symbTxWDM = np.zeros( + (len(t) // param.SpS, param.Nmodes, param.Nch), dtype='complex' + ) + + Psig = 0 - - for indMode in range(0, param.Nmodes): + + for indMode in range(param.Nmodes): print('Mode #%d'%(indMode)) - - for indCh in range(0, param.Nch): + + for indCh in range(param.Nch): # generate random bits bitsTx = np.random.randint(2, size=param.Nbits) @@ -79,9 +82,9 @@ def simpleWDMTx(param): # normalize symbols energy to 1 symbTx = symbTx/np.sqrt(Es) - + symbTxWDM[:,indMode,indCh] = symbTx - + # upsampling symbolsUp = upsample(symbTx, param.SpS) @@ -97,17 +100,17 @@ def simpleWDMTx(param): # optical modulation sigTxCh = iqm(Ai, 0.5*sigTx, Vπ, Vb, Vb) sigTxCh = np.sqrt(Pch/param.Nmodes)*sigTxCh/np.sqrt(signal_power(sigTxCh)) - + print('channel %d power : %.2f dBm, fc : %3.4f THz' %(indCh+1, 10*np.log10(signal_power(sigTxCh)/1e-3), (param.Fc+freqGrid[indCh])/1e12)) sigTxWDM[:,indMode] += sigTxCh*np.exp(1j*2*π*(freqGrid[indCh]/Fa)*t) - + Psig += signal_power(sigTxWDM[:,indMode]) - + print('total WDM signal power: %.2f dBm'%(10*np.log10(Psig/1e-3))) - + param.freqGrid = freqGrid - + return sigTxWDM, symbTxWDM, param