Commit 35c90120 authored by Arne Striegler's avatar Arne Striegler 😁

Version 3: Init

parent a2265e3f
......@@ -77,21 +77,22 @@ for c in range(0, 0):
############################
E = amp(E = E, Pmean = 0)
#Ein = copy.deepcopy(E)
print('1:', osnr( E = E))
Ein = copy.deepcopy(E)
plt.figure(1)
# Get decision matrix
Dec_x, Esx_re_ax, Esx_im_ax, Dec_y, Esy_re_ax, Esy_im_ax = get_decision_matrix(gp, E, [constpts_x, constpts_y], [symbols_x, symbols_y], filter_f0)
Dec_x, Esx_re_ax, Esx_im_ax, Dec_y, Esy_re_ax, Esy_im_ax = get_decision_matrix(gp, Ein, [constpts_x, constpts_y], [symbols_x, symbols_y], filter_f0)
print('1:', osnr( E = E))
# Plot decision matrix
plt.figure(1)
plt.subplot(2, 1, 1); h = plt.contourf(Esx_re_ax/1000, Esx_im_ax/1000, Dec_x, 32, cmap=plt.cm.bone )
plt.plot(np.real(E[0]['E'][0][int(gp.sps/2)::gp.sps]), np.imag(E[0]['E'][0][int(gp.sps/2)::gp.sps]), '.')
plt.subplot(2, 1, 2); h = plt.contourf(Esy_re_ax, Esy_im_ax, Dec_y, 32, cmap=plt.cm.bone )
plt.subplot(2, 1, 1); h = plt.contourf(Esx_re_ax, Esx_im_ax, Dec_x, 32, cmap=plt.cm.jet )
plt.plot(np.real(E[0]['E'][0][int(gp.sps/2)::gp.sps]), np.imag(E[0]['E'][0][int(gp.sps/2)::gp.sps]), '*')
plt.subplot(2, 1, 2); h = plt.contourf(Esy_re_ax, Esy_im_ax, Dec_y, 32, cmap=plt.cm.jet )
plt.plot(np.real(E[0]['E'][1][int(gp.sps/2)::gp.sps]), np.imag(E[0]['E'][1][int(gp.sps/2)::gp.sps]), '.')
############################
# Calculate BER
############################
#BER = calc_BER (gp, E, [constpts_x, constpts_y], osnr( E = E), Dec_x, Dec_y, Esx_re_ax, Esx_im_ax, Esy_re_ax, Esy_im_ax, 1, filter_f0, [symbols_x, symbols_y])
#print(BER)
BER = calc_BER (gp, E, [constpts_x, constpts_y], osnr( E = E), Dec_x, Dec_y, Esx_re_ax, Esx_im_ax, Esy_re_ax, Esy_im_ax, 10, filter_f0, [symbols_x, symbols_y])
print('1:', osnr( E = E))
print(BER)
plt.show()
\ No newline at end of file
......@@ -9,7 +9,7 @@ from pypho_setup import pypho_setup
from pypho_symbols import pypho_symbols
from pypho_signalsrc import pypho_signalsrc
from pypho_lasmod import pypho_lasmod
from pypho_fiber import pypho_fiber
from pypho_cfiber import pypho_fiber
from pypho_functions import *
import numpy as np
import copy
......
......@@ -352,7 +352,7 @@ class pypho_fiber(object):
y_max = self.glova.frange * 2.0 * np.pi / 1.0e12
x = np.arange(0,self.glova.sps*self.glova.nos)
k = 2.0 * np.pi *self.glova.fres / 1.0e12
self.y = x*k - ( (x*k-y_max/2.0)/np.abs(x*k-y_max/2.0) +1)*y_max/2.0
#self.y = x*k - ( (x*k-y_max/2.0)/np.abs(x*k-y_max/2.0) +1)*y_max/2.0
......
......@@ -109,8 +109,8 @@ def electrical_filter(glova, esig, B):
def get_decision_matrix(glova, E, constptsarray, symbols, ofil):
'Create decision matrix'
Esiggi = copy.deepcopy(E)
Esiggi = ofil(E = Esiggi )
#Esiggi = copy.deepcopy(E)
Esiggi = ofil(E = E )
#I_x_I, I_x_Q, I_y_I, I_y_Q = ninetydeghybrid(glova, Esiggi, LO) # Get photo currents
#Esiggi[0]['E'][0] = I_x_I +1.0j*I_x_Q # Now as electrial signal
......@@ -120,6 +120,8 @@ def get_decision_matrix(glova, E, constptsarray, symbols, ofil):
Esy = Esiggi[0]['E'][1][int(glova.sps/2)::glova.sps]
Nsx = symbols[0]
Nsy = symbols[1]
#print(Esx)
del(symbols)
N_raster = 1000.0 # Anzahl der Rasterpunkte
for c in [0, 1]:
......@@ -147,7 +149,7 @@ def get_decision_matrix(glova, E, constptsarray, symbols, ofil):
Es_im_ax = np.arange(Es_im_min, Es_im_max, (Es_im_max - Es_im_min) / N_raster )
Es_im_ax = Es_im_ax[0:int(N_raster)] # to be shure to have the correct number of cells
xx, yy = np.meshgrid(Es_re_ax, Es_im_ax, sparse = True)
for const_num in range(0, len(constpts[0])):
Z_tmp = np.zeros([int(N_raster), int(N_raster)])
......@@ -208,8 +210,9 @@ def ninetydeghybrid(gp, Esig, LO):
def create_optnoise (gp, P_sig, OSNR):
'Create noise vector'
E_N = np.zeros((2, gp.sps*gp.nos)) +0.0j
E_N = np.zeros((2, gp.sps*gp.nos)) + 0.0j
P_N = P_sig/10.0**(OSNR/10.0) * (gp.frange/12.5e9)
print(gp.frange, OSNR)
noisesamples = np.random.randn(4, gp.sps*gp.nos) * np.sqrt(P_N/4)
E_N[0] = noisesamples[0] + 1.0j*noisesamples[1]
E_N[1] = noisesamples[2] + 1.0j*noisesamples[3]
......@@ -222,6 +225,7 @@ def calc_BER (gp, E, constpts, OSNR, Dec_x, Dec_y, Esx_re_ax, Esx_im_ax, Esy_re_
'Calculate BER value of given signal'
Esiggi = copy.deepcopy(E)
del(E)
Nsx = symbols[0]
Nsy = symbols[1]
del(symbols)
......@@ -269,31 +273,32 @@ def calc_BER (gp, E, constpts, OSNR, Dec_x, Dec_y, Esx_re_ax, Esx_im_ax, Esy_re_
for n_samp in range( 0, len(Esx) ) :
# X-Pol
if int(Ex_x_pos[n_samp]) < 1000 and int(Ex_x_pos[n_samp]) >= 0 and int(Ex_y_pos[n_samp]) < 1000 and int(Ex_y_pos[n_samp] )>= 0:
const_num_ist = int( Dec_x[int(Ex_y_pos[n_samp]), int(Ex_x_pos[n_samp])] )
const_num_ist = int( Dec_x[int(Ex_y_pos[n_samp]), int(Ex_x_pos[n_samp])])
BER[0] += np.cumsum( np.ceil(np.add(constpts[0][1][const_num_ist], constpts[0][1][int(Nsx[n_samp])] ) % 2.0) )[-1]
if t == 0:
if np.cumsum( np.ceil(np.add(constpts[0][1][int(const_num_ist)] ,constpts[0][1][int(Nsx[n_samp])] ) % 2.0) )[-1] > 0:
#print(0, constpts[0][4][int(const_num_ist)] ,constpts[0][4][int(Nsx[n_samp])], BER[1], const_num_ist, Nsx[n_samp], n_samp)
if np.cumsum( np.ceil(np.add(constpts[0][1][int(const_num_ist)] ,constpts[0][1][int(Nsx[n_samp])] ) % 2.0) )[-1] > 0:
#print(0, constpts[0][4][int(const_num_ist)] ,constpts[0][4][int(Nsx[n_samp])], BER[1], const_num_ist, Nsx[n_samp], n_samp)
plt.subplot(2, 1, 1); plt.plot(np.real(Esx[n_samp]), np.imag(Esx[n_samp]), color=str(Nsx[n_samp]*1/len(constpts[0][1])), marker='x')
else:
plt.subplot(2, 1, 1); plt.plot(np.real(Esx[n_samp]), np.imag(Esx[n_samp]), color=str(Nsx[n_samp]*1/len(constpts[0][1])), marker='o', markersize=2)
plt.subplot(2, 1, 1); plt.plot(np.real(Esx[n_samp]), np.imag(Esx[n_samp]), color=str(1-Nsx[n_samp]*1/len(constpts[0][1])), marker='x')
else:
plt.subplot(2, 1, 1); plt.plot(np.real(Esx[n_samp]), np.imag(Esx[n_samp]), color=str(1-Nsx[n_samp]*1/len(constpts[0][1])), marker='o', markersize=2)
# Y-Pol
if int(Ey_x_pos[n_samp]) < 1000 and int(Ey_x_pos[n_samp]) >= 0 and int(Ey_y_pos[n_samp]) < 1000 and int(Ey_y_pos[n_samp]) >= 0:
const_num_ist = Dec_y[int(Ey_y_pos[n_samp]), int(Ey_x_pos[n_samp])]
BER[1] += np.cumsum( np.ceil(np.add(constpts[1][1][int(const_num_ist)], constpts[1][1][int(Nsy[n_samp])] ) % 2.0) )[-1]
# if np.cumsum( np.ceil(np.add(constpts[1][4][int(const_num_ist)] ,constpts[1][4][int(Nsy[n_samp])] ) % 2.0) )[-1] > 0:
# #print(1, constpts[1][4][int(const_num_ist)] ,constpts[1][4][int(Nsy[n_samp])], BER[1], const_num_ist, Nsy[n_samp], n_samp)
# plt.subplot(2, 1, 2); plt.plot(np.real(Esy[n_samp]), np.imag(Esy[n_samp]), color=constpts[1][2][int(Nsy[n_samp])], marker='x')
# else:
# plt.subplot(2, 1, 2); plt.plot(np.real(Esy[n_samp]), np.imag(Esy[n_samp]), color=constpts[1][2][int(Nsy[n_samp])], marker='o', markersize=2)
BER[1] += np.cumsum( np.ceil(np.add(constpts[1][1][int(const_num_ist)], constpts[1][1][int(Nsy[n_samp])] ) % 2.0) )[-1]
if t == 0:
if np.cumsum( np.ceil(np.add(constpts[1][1][int(const_num_ist)] ,constpts[1][1][int(Nsy[n_samp])] ) % 2.0) )[-1] > 0:
#print(1, constpts[1][4][int(const_num_ist)] ,constpts[1][4][int(Nsy[n_samp])], BER[1], const_num_ist, Nsy[n_samp], n_samp)
plt.subplot(2, 1, 2); plt.plot(np.real(Esy[n_samp]), np.imag(Esy[n_samp]), color=str(1-Nsx[n_samp]*1/len(constpts[0][1])), marker='x')
else:
plt.subplot(2, 1, 2); plt.plot(np.real(Esy[n_samp]), np.imag(Esy[n_samp]), color=str(1-Nsx[n_samp]*1/len(constpts[0][1])), marker='o', markersize=2)
BER[0] = np.log10(BER[0] / float(len(constpts[0][0]) * M * len(Esx) ))
......
......@@ -39,6 +39,8 @@ class pypho_optfi(object):
def __call__(self, E = None, Df = None, B = None, filtype = None, alpha = None, loss = None):
n_samples = self.glova.sps * self.glova.nos
if E == None:
print ("ERROR: pypho_optfi: You must define an optical signal")
sys.exit("PyPho stopped!")
......@@ -119,11 +121,11 @@ class pypho_optfi(object):
filfunc *= dbm2w(-1*self.loss)*1000.0
#print("--->"+ self.glova.wisdom_pyfftw)
pyfftw.import_wisdom(self.glova.wisdom_pyfftw)
input_array_f = np.zeros(self.glova.sps * self.glova.nos, dtype=np.complex128)
output_array_f = np.zeros(self.glova.sps * self.glova.nos, dtype=np.complex128)
input_array_f = np.zeros(n_samples, dtype=np.complex128)
output_array_f = np.zeros(n_samples, dtype=np.complex128)
fft_f = pyfftw.FFTW(input_array_f, output_array_f, direction='FFTW_FORWARD', flags=[self.glova.fftw3_plan], threads=self.glova.fftw3_threads)
input_array_b = np.zeros(self.glova.sps * self.glova.nos, dtype=np.complex128)
output_array_b = np.zeros(self.glova.sps * self.glova.nos, dtype=np.complex128)
input_array_b = np.zeros(n_samples, dtype=np.complex128)
output_array_b = np.zeros(n_samples, dtype=np.complex128)
fft_b = pyfftw.FFTW(input_array_b, output_array_b, direction='FFTW_BACKWARD', flags=[self.glova.fftw3_plan], threads=self.glova.fftw3_threads)
for Ei in E:
......@@ -131,27 +133,29 @@ class pypho_optfi(object):
pyfftw.FFTW.execute(fft_f)
input_array_b[:] = output_array_f * fftshift(np.sqrt( filfunc) )
pyfftw.FFTW.execute(fft_b)
E[z]['E'][0] = output_array_b[:]
E[z]['E'][0] = output_array_b[:] / n_samples
input_array_f[:] = E[z]['E'][1]
pyfftw.FFTW.execute(fft_f)
input_array_b[:] = output_array_f * fftshift(np.sqrt( filfunc) )
pyfftw.FFTW.execute(fft_b)
E[z]['E'][1] = output_array_b[:]
E[z]['E'][1] = output_array_b[:] / n_samples
# E[z]['E'][0] = ifft( ( fft(E[z]['E'][0]) ) * fftshift(np.sqrt( filfunc) ))
# E[z]['E'][1] = ifft( ( fft(E[z]['E'][1]) ) * fftshift(np.sqrt( filfunc) ))
E[z]['noise'] = E[z]['noise'] * filfunc
E[z]['noise'] *= filfunc
z += 1
# TODO : Freemem FFTW
# TODO : for Ei in E ändern in z 0 to length.E
# fig=plt.figure('filter')
# ax = fig.add_subplot(111)
# ax.plot((self.glova.freqax()-self.glova.f0)*1e-9, filfunc)
# ax.set_ylim(-0.2, 1.2)
#fig=plt.figure('filter')
#ax = fig.add_subplot(111)
#ax.plot((self.glova.freqax()-self.glova.f0)*1e-9, filfunc)
#ax.plot((self.glova.freqax() - self.glova.f0) * 1e-9, np.abs(input_array_f[:]))
#ax.set_ylim(-0.2, 1.2)
return E
......
......@@ -38,7 +38,7 @@ class pypho_setup(object):
self.pswd = "testuser";
self.key = "testpwd";
self.fftw3_threads = 4
self.fftw3_threads = 8
self.fftw3_plan = 'FFTW_PATIENT' # FFTW_ESTIMATE, FFTW_MEASURE, FFTW_PATIENT, FFTW_EXHAUSTIVE
self.fftw3_path_pyfftw = "../config/wisdom_pyfftw_" + str(nos * sps); #" + hex(uuid.getnode()) + "
self.fftw3_path_fftw3 = "../config/wisdom_fftw3.txt"; #+ hex(uuid.getnode()) + "_"
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment