Four wave mixing (FWM) is a non linear fiber effects, which produces mixing products. The mixing products strongly depends on the power, the non linear coefficent $\gamma$, the chromatic dispersion $D(\lambda)$, polarisation and the wavelength difference of the input signals.
As can be seen in the figures below the example shows the input signal spectrum, which consists of 5 cw-signals, where the signal in the middle has a different polarisation than it's neighbours.
At the end of the fiber the spectrum shows additional channels, which result from FWM.
Fig. 1 : Input and output spectrum of the X- (red) and Y- (blue) polarisation with four wave mixing products
Figure 2 shows how the mixng products evolve over the transmssion distance.
Fig. 2 : Power of the four wave mixing products as function of the transmission distance.
The figures are plotted by the script below. Feel free to play around with the polarisation, channel spacing and chromatic dispersion coefficient.
#Import functions and libraries
import sys
sys.path.append('../')
from pypho_setup import pypho_setup
from pypho_fiber import pypho_fiber
from pypho_cwlaser import pypho_cwlaser
from pypho_optfi import pypho_optfi
from pypho_functions import *
import numpy as np
import copy
import matplotlib.patches as mpatches
import matplotlib.pyplot as plt
# Define some paramaters
delta_f = 25 # Frequency difference of the two cw-signals
L = 20e3 # Fiber length in m
l = 0.1e3 # Step fiber length in m
k = 7 # Defines how man channels left and right will be analyzed
chut = np.arange(-k*delta_f, k*delta_f+delta_f, delta_f) # Channels under test
P_FWM_x = np.zeros((chut.size, int(L/l)))
P_FWM_y = np.zeros((chut.size, int(L/l)))
# Define network elements
gp = pypho_setup(nos = 4*16, sps = 1*128, symbolrate = 10e9)
sig_f0 = pypho_cwlaser(glova = gp, power = 3, Df = 0, teta = 1*np.pi/4.0)
sig_f = pypho_cwlaser(glova = gp, power = 3, Df = delta_f, teta = 1*np.pi/2.0)
filter_df = pypho_optfi(glova = gp, Df = 0, B = 10)
SSMF = pypho_fiber(glova = gp, l = l, D = 5.0, S = 0.0, alpha = 0.2, gamma = 1.4, phi_max = .1)
# Simulation
# Define wavelength channel
E_f0 = sig_f0()
E_fm1 = sig_f(Df = +delta_f, teta = np.pi/2)
E_fp1 = sig_f(Df = -delta_f, teta = np.pi/2)
E_fm2 = sig_f(Df = -2*delta_f, teta = np.pi/2)
E_fp2 = sig_f(Df = +2*delta_f, teta = np.pi/2)
E = copy.deepcopy(E_fm1)
E[0]['E'][0] = E_f0[0]['E'][0] + E_fp1[0]['E'][0] + E_fm1[0]['E'][0] + E_fp2[0]['E'][0] + E_fm2[0]['E'][0] # Multiplex all signals X-Pol
E[0]['E'][1] = E_f0[0]['E'][1] + E_fp1[0]['E'][1] + E_fm1[0]['E'][1] + E_fp2[0]['E'][1] + E_fm2[0]['E'][1] # Multiplex all signals Y-Pol
E_Tx = copy.deepcopy(E)
# Plot Spectrum Input signals
plt.figure(1)
plt.subplot(2, 1, 1)
plt.plot((gp.freqax-gp.f0)*1e-9, np.log10(abs(fftshift(fft(E[0]['E'][0] )**2))), 'r')
plt.plot((gp.freqax-gp.f0)*1e-9, np.log10(abs(fftshift(fft(E[0]['E'][1] )**2))), 'g:')
plt.title("Input spectrum", loc='left')
plt.ylabel('Spec. density');
plt.grid(True)
# Fiber transmission
c1 = 0
for z in np.arange(0,L,l):
print(z)
c2 = 0
for ch in chut: # Loop through wavelengths and save mean power value of both pol axis
E_filtered = filter_df(copy.deepcopy(E), Df = ch)
P_FWM_x[c2][c1] = np.mean( np.abs(E_filtered[0]['E'][0])**2 )
P_FWM_y[c2][c1] = np.mean( np.abs(E_filtered[0]['E'][1])**2 )
c2 = c2 + 1
c1 = c1 + 1
E = SSMF(E = E)
# Plot Spectrum Output signals
plt.subplot(2, 1, 2)
plt.plot((gp.freqax-gp.f0)*1e-9, np.log10(abs(fftshift(fft(E[0]['E'][0] )**2))), 'r')
plt.plot((gp.freqax-gp.f0)*1e-9, np.log10(abs(fftshift(fft(E[0]['E'][1] )**2))), 'g:')
plt.title("Output spectrum", loc='left')
plt.ylabel('Spec. density'); plt.xlabel('Frequency deviation [GHz]');
plt.grid(True)
# Plot power of both pol axis as function of transmission distance
plt.figure(2)
for i in np.arange(0,chut.size):
plt.subplot(2, 1, 1)
plt.plot(np.arange(0,L,l), 10*np.log10(P_FWM_x[i]*1e3), label=str(chut[i]) + ' GHz')
plt.subplot(2, 1, 2)
plt.plot(np.arange(0,L,l), 10*np.log10(P_FWM_y[i]*1e3))
plt.subplot(2, 1, 2)
plt.ylabel('$10log |E_y|^2$'); plt.xlabel('Transmission distance [m]');
plt.ylim((-120, 10))
plt.grid(True)
plt.subplot(2, 1, 1)
plt.ylabel('$10log |E_y|^2$');
plt.ylim((-120, 10))
plt.grid(True)
plt.legend()
plt.show()