Evolution

import numpy as np
import random
import matplotlib.pyplot as plt
from scipy.fft import fft, ifft
from scipy import integrate
from ModulationPy import PSKModem
from functools import partial

def gen(N): #генерация точек из созвездия QAM
    numbers = [-3/10,-1/10,1/10,3/10]
    q = []
    for i in range(0, N):
        Q = random.randint(0, len(numbers)-1)
        I = random.randint(0, len(numbers)-1)
        q.append(complex(numbers[Q], numbers[I]))
    return q



modalpha = []
modbeta = []
M = 20
q = gen(M)
im = complex(0, 1)
Tcontainer = []
Zcontainer = []
condlist = []
zeta = []
#///////////////////апсеймплинг
num_symbols = 1
sps = 4
num_taps = 101
beta = 0.35
Ts = sps
t = np.arange(-50, 51)
h = 1/Ts*np.sinc(t/Ts) * np.cos(np.pi*beta*t/Ts) / (1 - (2*beta*t/Ts)**2)
plt.figure(1)
plt.plot(t, h, '.')
plt.grid(True)
plt.show()
q = np.convolve(q, h) #апсеймплинг
plt.figure(2)
plt.plot(q, '.-')
for i in range(num_symbols):
    plt.plot([i*sps+num_taps//2+1,i*sps+num_taps//2+1], [0, q[i*sps+num_taps//2+1]])
plt.grid(True)
plt.show()

q_re = []
for k in q:
    q_re.append(k.real)
plt.plot(np.real(q), np.imag(q), '.')
plt.grid(True)
plt.show()

plt.plot(q_re)
plt.grid(True)
plt.show()

f = fft(q_re)
plt.plot(f.real)
plt.grid(True)
plt.show()
#///////////////////апсеймплинг

#///////////////////прямая задача
T1 = 0 #начальное значение по x
T2 = 5 #конечное значение по x
M = len(q)
eps = (T2 - T1) / (M-1) #шаг
print(eps)
beta_list = []
alpha_list = []
delta_zeta = -np.pi/eps
delta = (np.abs(-np.pi)+np.abs(np.pi))/M
while delta_zeta < np.pi/eps:
    Z = np.exp(-delta_zeta * im * eps)
    Zcontainer.append(Z)
    T = np.matrix([[1, 0], [0, 1]])
    for m in range(1, M): #алгоритм нахождения коэффициентов рассеяния
        addition = 1 / np.sqrt(1 + eps ** 2 * np.abs(q[m]) ** 2)
        tempT = addition * np.matrix([[Z, eps * q[m]],
                                   [-eps * complex.conjugate(q[m]), Z ** (-1)]])
        T = tempT.dot(T)

    Tcontainer.append(T)
    diagmatrix = np.matrix([[np.exp(im * delta_zeta * (T2 + 0.5 * eps)), 0],
                            [0, np.exp(-1 * im * delta_zeta * (T2 + 0.5 * eps))]])
    result = diagmatrix.dot(T) * np.exp(-1 * delta_zeta * im * (T1 - 0.5 * eps))
    alpha = result[0, 0]
    beta = result[1, 0]
    test = np.abs(alpha)**2 + np.abs(beta)**2
    # print("zeta:", delta_zeta, "result:", alpha, beta)
    # print("zeta:", delta_zeta, "condition:", test)
    delta_zeta += delta
    beta_list.append(beta)
    alpha_list.append(alpha)
    condlist.append(test)
    zeta.append(delta_zeta.real)
    modalpha.append(np.abs(alpha))
    modbeta.append(np.abs(beta))

#///////////////////прямая задача

fig, ax = plt.subplots()
plt.plot(zeta, condlist, 'blue', linewidth=2, label='x1')
ax.set_xlabel('$Re zeta$')
ax.set_ylabel('$|a|^2 + |b|^2 $')
plt.grid()
plt.show()

fig, ax = plt.subplots()
plt.plot(zeta, modalpha, 'blue', linewidth=2, label='x1')
ax.set_xlabel('$Re zeta$')
ax.set_ylabel('$|a|^2 + |b|^2 $')
plt.grid()
plt.show()

fig, ax = plt.subplots()
plt.plot(zeta, modbeta, 'orange', linewidth=2, label='x1')
ax.set_xlabel('$Re zeta$')
ax.set_ylabel('$|a|^2 + |b|^2 $')
plt.grid()
plt.show()


fig, (ax1, ax2) = plt.subplots(2, 1)
fig.subplots_adjust(hspace=0.5)
ax1.plot(zeta, modalpha, zeta, modbeta)
ax1.set_xlim(0, 5)
ax1.set_xlabel('zeta')
ax1.set_ylabel('modules')
ax1.grid(True)
plt.show()


mincond = min(condlist)
maxcond = max(condlist)
print(f"condition minimum value: {mincond}")
print(f"condition maximum value: {maxcond}")

#///////////////////эволюция
evol_beta = []
bar_beta = []
t = 2
for k in range(len(beta_list)): #эволюция коэффициентов рассеяния
    evol_beta.append(beta_list[k] * np.exp(-4j * zeta[k] ** 2 * t))


for k in range(len(evol_beta)):
    bar_beta.append(-complex.conjugate(evol_beta[k]))

bar_alpha = []
for k in range(len(alpha_list)):
    bar_alpha.append(complex.conjugate(alpha_list[k]))

L = T1 + (M)*eps
print(M)
print(L)
phi = []
bar_phi = []
for k in range(len(zeta)):
    phi.append(alpha_list[k] * np.array([np.exp(-1j * zeta[k] * L), 0]) +
               evol_beta[k] * np.array([0, np.exp(1j * zeta[k] * L)]))
    bar_phi.append(-bar_alpha[k] * np.array([0, np.exp(1j * zeta[k] * L)]) +
                   bar_beta[k] * np.array([np.exp(-1j * zeta[k] * L), 0]))

newT = []
for k in range(len(zeta)):
    temp = np.array([[phi[k][0], bar_phi[k][0]], [phi[k][1], bar_phi[k][1]]])
    newT.append(temp)


Tmodif = newT
Inverseq = []

#///////////////////обратная задача для данных после эволюции

Z = []
Z_re = []
dzeta = -np.pi/eps
delta1 = (np.abs(-np.pi)+np.abs(np.pi))/M
newzeta = []
while dzeta < np.pi/eps:
    Z.append(np.exp(-1j * dzeta * eps))
    Z_re.append(np.exp(-1j * dzeta * eps).real)
    dzeta += delta1
    newzeta.append(dzeta)

for j in range(M, 0, -1): #Восстановление q(x,t) по матрице, содержащей коэффициенты рассеяния
    matcoef_container = []
    for k in range(0, j):
        a = np.array(Tmodif)[:, 0, 0] * np.array(Z) ** (-k)
        b = np.array(Tmodif)[:, 0, 1] * np.array(Z) ** (-k)
        c = np.array(Tmodif)[:, 1, 0] * np.array(Z) ** (-k)
        d = np.array(Tmodif)[:, 1, 1] * np.array(Z) ** (-k)
        Integral1 = 1 / (2 * np.pi * 1j) * np.trapz(a, zeta)
        Integral2 = 1 / (2 * np.pi * 1j) * np.trapz(b, zeta)
        Integral3 = 1 / (2 * np.pi * 1j) * np.trapz(c, zeta)
        Integral4 = 1 / (2 * np.pi * 1j) * np.trapz(d, zeta)
        matcoef = np.array([[Integral1, Integral2], [Integral3, Integral4]])
        matcoef_container.append(matcoef)

    coef1 = complex.conjugate(matcoef_container[len(matcoef_container)-1][0, 0])
    coef2 = complex.conjugate(matcoef_container[len(matcoef_container)-2][1, 0])
    res = -eps ** (-1) * coef2 / coef1
    Inverseq.append(res)
    InverseContainer = []
    InverseCoefCont = []
    addition = 1 / np.sqrt(1 + eps ** 2 * np.abs(Inverseq[M-j]) ** 2)
    for k in range(len(Tmodif)):
        tempInv = addition * np.array([[Z[k] ** (-1), -eps*Inverseq[M-j]],
                                       [eps * complex.conjugate(Inverseq[M-j]), Z[k]]])
        InverseContainer.append(tempInv)

    for k in range(len(Tmodif)):
        Tmodif[k] = InverseContainer[k].dot(Tmodif[k])


#///////////////////обратная задача для данных после эволюции

print("------")

for k in Inverseq:
    print(k)
print("------")
print(len(q))
print(len(Inverseq))
qinv = Inverseq[::-1]
plt.plot(q, label="q - generated")
plt.plot(qinv, label="q - ILP")
plt.grid()
plt.legend()
plt.show()

#///////////////////SSFM

t = 0
u = q
dt = 0.001
tfinal = 2
H = round(tfinal / dt)
N = len(q)
k = np.array(np.arange(-N/2, N/2, 1))
w = k * 2 * np.pi / N / eps
for t in range(0, H): #Алгоритм метода Фурье с разделённым шагом
  u = np.exp(2j * np.abs(u) ** 2 * dt)*u
  a = np.fft.fftshift(fft(u))
  a = np.exp(1j * w * w * dt) * a
  u = ifft(np.fft.ifftshift(a))
  t += dt




plt.plot(u, label=f"q(x,{tfinal}) - SSFM")
plt.plot(qinv, label=f"q(x,{tfinal}) - ILP")
#plt.plot(q, label="q(x,0)- generated")
plt.grid()
plt.legend()
plt.show()
#///////////////////SSFM

#///////////////////проверква коэффициентов рассеяния для найденного #q(x,t)(С помощью МОЗР)

modalpha2 = []
modbeta2 = []
beta2 = []
delta_zeta = -np.pi/eps
delta = (np.abs(-np.pi)+np.abs(np.pi))/M
while delta_zeta < np.pi/eps:
    Z = np.exp(-delta_zeta * im * eps)
    Zcontainer.append(Z)
    T = np.matrix([[1, 0], [0, 1]])
    for m in range(1, M):
        addition = 1 / np.sqrt(1 + eps ** 2 * np.abs(qinv[m]) ** 2)
        tempT = addition * np.matrix([[Z, eps * qinv[m]],
                                   [-eps * complex.conjugate(qinv[m]), Z ** (-1)]])
        T = tempT.dot(T)

    Tcontainer.append(T)
    diagmatrix = np.matrix([[np.exp(im * delta_zeta * (T2 + 0.5 * eps)), 0], [0, np.exp(-1 * im * delta_zeta * (T2 + 0.5 * eps))]])
    result = diagmatrix.dot(T) * np.exp(-1 * delta_zeta * im * (T1 - 0.5 * eps))
    alpha = result[0, 0]
    beta = result[1, 0]
    test = np.abs(alpha)**2 + np.abs(beta)**2
    print("zeta:", delta_zeta, "result:", alpha, beta)
    print("zeta:", delta_zeta, "condition:", test)
    delta_zeta += delta
    condlist.append(test)
    modalpha2.append(np.abs(alpha))
    beta2.append(beta.real)
    modbeta2.append(np.abs(beta))



mincond = min(condlist)
maxcond = max(condlist)
print(f"condition minimum value: {mincond}")
print(f"condition maximum value: {maxcond}")

plt.plot(modalpha)
plt.plot(modalpha2)
plt.show()

plt.plot(modbeta)
plt.plot(beta2)
plt.show()

plt.plot(evol_beta)
plt.plot(beta2)
plt.show()

print(len(beta2))
print(len(evol_beta))

q(x,t) и q(x,0):

SSFM(синий), МОЗР(оранжевый), q(x,0)(зеленый)

сравнение коэффициентов q(x,0) и q(x,t) (МОЗР)

синий a(zeta) до эволюции

оранжевый a(zeta) после

Действительные части b(zeta)

синий эволюция b(zeta)

оранжевый получен прямой задачей для q(x,t)