added:pypho_functions

pypho_functions.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+#  pypho.py
+#
+#  Copyright 2014 Arne Striegler (arne.striegler@hm.edu)
+#
+#
+#
+# Functions
+#
+#
+########################################################################
+
+import numpy as np
+import scipy.fftpack
+#import pyfftw
+import urllib2
+import json
+import os
+#import pycurl
+from StringIO import StringIO
+from pypho_constants import pypho_constants
+import zlib
+import zipfile
+cv = pypho_constants();
+import copy 
+import matplotlib.pyplot as plt
+
+########################################################################
+def dbm2w( power_dbm ):
+   "Transfers power value from dBm to W"
+   return  10**(power_dbm/10.0) / 1000.0
+
+
+########################################################################
+def w2dbm( power_w ):
+   "Transfers power value from W to dBm"
+   return 10*np.log10(power_w*1e3)
+
+
+########################################################################
+def polrot(sig, alpha):
+    "Rotates the polarisation of the signal"
+    sigout = np.cos(alpha)*sig[0,:] - np.sin(alpha)*sig[1,:]
+    sigout = np.vstack((sigout, np.sin(alpha)*sig[0,:] + np.cos(alpha)*sig[1,:]))
+    return sigout
+
+########################################################################
+def db2neper(alpha_dB):
+    "Converts attenuation parameter from dB/km into 1/km"
+    return alpha_dB/4.3429
+
+########################################################################
+def DS2beta(D, S, lambda_ref):
+    "Converts D and S into beta_2 [s**2 / km] and beta_3 [s** 3 /km]"
+    if ( D == 0.0) and (S == 0.0):
+        return (0.0, 0.0)
+
+    D *= 1.0e-6
+    S *= 1.0e3
+    fact = lambda_ref**2 / (2.0 * np.pi * cv.lightspeed)
+
+    return (-fact*D, fact**2 * ( S + (2.0 * D / lambda_ref)))
+
+########################################################################
+def fft(E):
+    "Calculate FFT"
+    #fft_object = pyfftw.FFTW(E, b)
+    return scipy.fftpack.fft(E)
+
+########################################################################
+def ifft(E):
+    "Calculate inverse FFT"
+    #fft_object = pyfftw.FFTW(E, b,'FFTW_BACKWARD')
+    #return b
+    return scipy.fftpack.ifft(E)
+
+########################################################################
+def fftshift(E):
+    "Shift the zero-frequency component to the center of the spectrum"
+    return scipy.fftpack.fftshift(E)
+
+########################################################################
+def supresslowpow(E):
+    "Supress all array alements < x dBm"
+
+    E[ np.where(np.log10(np.abs(E)**2)) < -10 ] = 1e-30 + 0*1j
+
+    return E
+
+########################################################################
+def getpower_dBm(E):
+    "Get mean power in dBm"
+    return (w2dbm(np.mean(np.abs(E[0,:])**2)), w2dbm(np.mean(np.abs(E[1,:])**2)) )
+
+########################################################################
+def getpower_W(E):
+    "Get mean power in W"
+    return (np.mean(np.abs(E[0,:])**2), np.mean(np.abs(E[1,:])**2))
+
+########################################################################
+def cloud_init_E(gp, E):
+    "Adjust E if cloud calculation is done"
+    if gp.cloud :
+        return E, []
+    else :
+        return E, E
+
+########################################################################
+def cloud_SaveFiles(gp, E, link):
+    "Save, zip and send all files for cloude computing"
+    
+    # Send to the magic pycloud
+    # Ask for sim ticket
+    
+    cid ="none"
+    body1 = "none"
+    if gp.cloud :
+        buffer = StringIO()
+        c = pycurl.Curl()
+        c.setopt(c.URL, 'https://optical-fiber-transmission-simulation.com/getid/')
+        c.setopt(c.WRITEFUNCTION, buffer.write)
+        c.setopt(c.USERPWD, gp.key + ':' + gp.pswd)
+    
+        c.perform()
+        c.close()
+    
+        body = buffer.getvalue()
+    
+        print ('I got a sim ticket...')
+        cid =  body
+        print (cid)
+        # Create files
+        import json
+    
+        glova = {'simid': cid, 'nob' : gp.nob, 'spb' : gp.spb, 'f' : gp.f, 'lambda0' : gp.lambda0, 'bitrate' : gp.bitrate, 'polarisation' : gp.polarisation}
+        with open('glova.json', 'w') as fp:
+            json.dump(glova, fp)
+
+        with open('link.json', 'w') as fp:
+            json.dump(link, fp)
+
+        np.savetxt('E.txt', E[0]['E'].view(float) )
+
+
+        try:
+            compression = zipfile.ZIP_DEFLATED
+        except:
+            compression = zipfile.ZIP_STORED
+    
+        print 'Creating archive...'
+        zf = zipfile.ZipFile('pypho.zip', mode='w')
+        try:
+            zf.write('E.txt', compress_type=compression)
+            zf.write('glova.json', compress_type=compression)
+            zf.write('link.json', compress_type=compression)
+        finally:
+            print 'Ready to go!'
+            zf.close()
+        
+        
+        
+        # Send files to cloud
+        print ('Sending files to pyclo!')
+        c = pycurl.Curl()
+        buffer1 = StringIO()
+        c.setopt(c.URL, 'https://optical-fiber-transmission-simulation.com/ul/upload.php')
+        c.setopt(c.POST, 1)
+        c.setopt(c.HTTPPOST, [("sim_ticket", cid), ("filecontents", (c.FORM_FILE, "pypho.zip"))])
+    
+        c.setopt(c.WRITEDATA, buffer1)
+        c.perform()
+        c.close()
+        body1 = buffer1.getvalue()
+    
+    return (cid, body1)
+
+
+########################################################################
+def cloud_GetStatus(sim_id, item="*"):
+    "Ask for status information of a simulation"
+    
+    url = 'https://optical-fiber-transmission-simulation.com/getstatus/?sim_ticket='+sim_id+'&item='+item
+
+    content = urllib2.urlopen(url).read()
+
+    return json.loads(content)
+    
+
+########################################################################
+def cloud_GetFile(sim_id, filename):
+    "Download file from cloud"
+    path = "efiles/"
+    retvalue = "OK"
+    url = 'https://optical-fiber-transmission-simulation.com/getfile/?sim_ticket='+sim_id+'&filename='+filename
+
+    try:
+        f = urllib2.urlopen(url, path + sim_id + '_' + filename)
+        print "Downloading " + sim_id
+
+        # Open our local file for writing
+        #with open(os.path.basename(path + sim_id + '_' + filename), "wb") as local_file:
+        with open( path + sim_id + '_' + filename, "wb") as local_file:
+            local_file.write(f.read())
+            local_file.close()
+            print ( path + sim_id + '_' + filename)
+            zfobj = zipfile.ZipFile(path + sim_id + '_' + filename)
+            for name in zfobj.namelist():           
+                uncompressed = zfobj.read(name)
+            
+                # save uncompressed data to disk
+                outputFilename = path + sim_id + '_' + name
+                print "Saving extracted file to ", outputFilename
+                output = open(outputFilename,'wb')
+                output.write(uncompressed)
+                output.close()  
+
+
+    #handle errors
+    except urllib2.HTTPError as e:
+        print "HTTP Error:", e.code, url
+        retvalue = "HTTP Error:", e.code, url
+    except urllib2.URLError as e:
+        print "URL Error:", e.reason, url
+        
+    return retvalue
+
+########################################################################
+def electrical_filter(glova, esig, B):
+    'Get Gauß filterd electrical Signal'
+    sig = []
+    filfunc = np.exp(-((glova.freqax + (-glova.f0))/(B*2.0e9 /2.3582))**2/2.0)**2   
+    sig = ifft((fft(esig)) * fftshift(filfunc))
+        
+    return sig
+
+########################################################################
+def get_decision_matrix(glova, E, constptsarray, bits, LO, ofil, sigsampler):
+    'Create decision matrix'
+    
+    Esiggi = ofil(E = copy.deepcopy(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
+    Esiggi[0]['E'][1] = I_y_I +1.0j*I_y_Q        # Now as electrial signal
+    
+    Esx, Esy, Nsx, Nsy = sigsampler(E = Esiggi, constpoints = constptsarray, bits = bits)
+
+    N_raster    = 1000.0 # Anzahl der Rasterpunkte
+ 
+    for c in [0, 1]:
+        
+        Dec         = np.zeros([int(N_raster), int(N_raster)])        
+        Z_tmp       = np.zeros([int(N_raster), int(N_raster)])
+        Z_max       = np.zeros([int(N_raster), int(N_raster)])
+    
+        if 0 == c: 
+            Es = Esx
+            Ns = Nsx
+            constpts = constptsarray[0]
+        else:
+            Es = Esy
+            Ns = Nsy
+            constpts = constptsarray[1]
+            
+        Es_re_min = 2.0 * np.min(np.real(Es))
+        Es_re_max = 2.0 * np.max(np.real(Es))
+        Es_im_min = 2.0 * np.min(np.imag(Es))
+        Es_im_max = 2.0 * np.max(np.imag(Es))
+        
+        Es_re_ax = np.arange(Es_re_min, Es_re_max, (Es_re_max - Es_re_min) / N_raster )
+        Es_re_ax = Es_re_ax[0:int(N_raster)]  # to be shure to have the correct number of cells
+        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)])
+        
+            for n_samp in range(0, len(Es)):
+                
+                if (const_num == Ns[n_samp]):
+                    Z_tmp += np.exp( -( (np.real(Es[n_samp])-xx)**2 + (np.imag(Es[n_samp])-yy)**2 ) / 0.050 )
+            
+
+            Z_tmp /= np.max(Z_tmp)
+            Z_dec_tmp = np.where(Z_max < Z_tmp)    
+            Dec[Z_dec_tmp] = const_num    
+            Z_max[Z_dec_tmp] = Z_tmp[Z_dec_tmp]
+        
+        if 0 == c: 
+            Dec_x       = copy.deepcopy(Dec)
+            Esx_re_ax   = copy.deepcopy(Es_re_ax)
+            Esx_im_ax   = copy.deepcopy(Es_im_ax)
+        else:
+            Dec_y       = copy.deepcopy(Dec)
+            Esy_re_ax   = copy.deepcopy(Es_re_ax)
+            Esy_im_ax   = copy.deepcopy(Es_im_ax)
+        
+    return Dec_x, Esx_re_ax, Esx_im_ax, Dec_y, Esy_re_ax, Esy_im_ax
+
+########################################################################
+def ninetydeghybrid(gp, Esig, LO):
+    '90 deg hybrid'
+
+    for pol in [1, 0]:
+        E_1 = 0.5 * ( Esig[0]['E'][pol] +      LO[0]['E'][pol] )
+        E_2 = 0.5 * ( Esig[0]['E'][pol] -      LO[0]['E'][pol] )
+        E_3 = 0.5 * ( Esig[0]['E'][pol] + 1.0j*LO[0]['E'][pol] )
+        E_4 = 0.5 * ( Esig[0]['E'][pol] - 1.0j*LO[0]['E'][pol] )
+
+        # Electrical filter
+        
+        E_1 = electrical_filter(glova = gp, esig = np.abs(E_1)**2, B = gp.symbolrate*2.0e-9)
+        E_2 = electrical_filter(glova = gp, esig = np.abs(E_2)**2, B = gp.symbolrate*2.0e-9)
+        E_3 = electrical_filter(glova = gp, esig = np.abs(E_3)**2, B = gp.symbolrate*2.0e-9)
+        E_4 = electrical_filter(glova = gp, esig = np.abs(E_4)**2, B = gp.symbolrate*2.0e-9)
+
+                
+        I_I =  E_1 - E_2
+        I_Q =  E_3 - E_4
+       
+        
+        if pol == 0:
+            I_x_I = I_I
+            I_x_Q = I_Q
+        else:
+            I_y_I = I_I
+            I_y_Q = I_Q
+            
+    return I_x_I, I_x_Q, I_y_I, I_y_Q
+
+########################################################################
+def create_optnoise (glova, E, OSNR):
+    'Create noise vector'
+    E_N = copy.deepcopy(E[0]['E'])
+    P_sig = np.mean(abs(E[0]['E'][0]**2) + abs(E[0]['E'][1]**2))
+    P_N = P_sig/10**(OSNR/10) * (glova.frange/12.5e9)
+    noisesamples = np.random.randn(4, len(E[0]['E'][0])) * np.sqrt(P_N/4)
+    E_N[0] = noisesamples[0] + 1j*noisesamples[1]
+    E_N[1] = noisesamples[2] + 1j*noisesamples[3]
+    
+    return E_N
+
+########################################################################
+def calc_BER (gp, E, constpts, OSNR, Dec_x, Dec_y, Esx_re_ax, Esx_im_ax, Esy_re_ax, Esy_im_ax, M, LO, ofil, sigsampler):
+    'Calculate BER value of given signal'
+    
+    Esiggi = copy.deepcopy(E)
+    Esx_re_ax_max = Esx_re_ax[1] - Esx_re_ax[0]
+    Esx_im_ax_max = Esx_im_ax[1] - Esx_im_ax[0]
+    Esy_re_ax_max = Esy_re_ax[1] - Esy_re_ax[0]
+    Esy_im_ax_max = Esy_im_ax[1] - Esy_im_ax[0]    
+    
+    BER=[0,0]
+    
+    for t in range(0, M):
+        
+        if (int(t/float(M)*100.0)- int((t-1)/float(M)*100.0)) > 0:
+            print('Progress: ', t/float(M)*100.0, '%')
+        
+        #create noise vectors
+        E_tmp   = copy.deepcopy(Esiggi)
+        E_N = create_optnoise (gp, Esiggi, OSNR) 
+            
+        # add noise
+        E_tmp[0]['E'][0] +=  E_N[0]
+        E_tmp[0]['E'][1] +=  E_N[1]
+        
+        E_tmp = ofil(E = E_tmp )
+    
+        # detect signal with noise
+        I_x_I, I_x_Q, I_y_I, I_y_Q = ninetydeghybrid(gp, E_tmp, LO)
+        E_tmp[0]['E'][0] = I_x_I + 1.0j*I_x_Q
+        E_tmp[0]['E'][1] = I_y_I + 1.0j*I_y_Q
+        
+        Esx, Esy, Nsx, Nsy = sigsampler(E = E_tmp, style = 'NO')
+
+       # plt.figure(1)
+        #plt.subplot(2, 1, 1); plt.plot(np.real(Esx), np.imag(Esx), '.')
+        #plt.subplot(2, 1, 2); plt.plot(np.real(Esy), np.imag(Esy), '.')
+      
+        # Calc BER            
+        Ex_x_pos = np.floor( (np.real(Esx) - Esx_re_ax[0]) / Esx_re_ax_max )
+        Ex_y_pos = np.floor( (np.imag(Esx) - Esx_im_ax[0]) / Esx_im_ax_max )
+        Ey_x_pos = np.floor( (np.real(Esy) - Esy_re_ax[0]) / Esy_re_ax_max )
+        Ey_y_pos = np.floor( (np.imag(Esy) - Esy_im_ax[0]) / Esy_im_ax_max )
+    
+    
+        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 = Dec_x[int(Ex_y_pos[n_samp]), int(Ex_x_pos[n_samp])]
+                BER[0] += np.cumsum( np.ceil(np.add(constpts[0][4][int(const_num_ist)], constpts[0][4][int(Nsx[n_samp])] ) % 2.0) )[-1]        
+             #   if np.cumsum( np.ceil(np.add(constpts[0][4][int(const_num_ist)] ,constpts[0][4][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=constpts[0][2][int(Nsx[n_samp])], marker='x')  
+             #   else:
+             #       plt.subplot(2, 1, 1); plt.plot(np.real(Esx[n_samp]), np.imag(Esx[n_samp]), color=constpts[0][2][int(Nsx[n_samp])], 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][4][int(const_num_ist)], constpts[1][4][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[0] = np.log10(BER[0] / float(len(constpts[0][0]) * M * len(Esx)   ))
+    BER[1] = np.log10(BER[1] / float(len(constpts[1][0]) * M * len(Esy)   ))
+        
+    return BER
\ No newline at end of file