<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 OFDM-Synchronisation nach Schmidl u. Cox
Beschrieben in
Step1:  Präambel generieren (BPSK)
Step2:  IFFT
Step3:  Cyclic Prefix einfügen
Step4:  Kanal
Step5:  Synchronisation
Step6:  Entzerrung
Anmerkung
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import random
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
matplotlib.rcParams.update({'font.size': 14})
L = 1024    # Anzahl Punkte in der DFT
K = 3*L//4    # Anzahl belegter Unterträger, muss gerade sein
n_CP = L//4  # Länge des Cyclic Prefix
threshold = 0.8  # Threshold for presence detection
t_off = 1000  # Zeitoffset
f_off = 11.3  # Frequenzoffset in Vielfachen des Unterträgerabstands
h_awgn = [1.,]   # Kanalimpulsantwort für AWGN
h_multipath = np.linspace(1, 0, n_CP//2) * np.exp(1j*np.random.uniform(-np.pi, np.pi, n_CP//2))  # Mehrwegekanal
h_multipath /= np.sqrt(np.mean(np.abs(h_multipath)**2))  # normieren
h = h_awgn  # Kanal wählen
SNR = 10     # SNR in dB
Explanation: OFDM-Synchronisation nach Schmidl u. Cox
Beschrieben in: T. Schmidl u. D. Cox, Robust Frequency and Timing Synchronization for OFDM, IEEE Transactions on Communications, December 1997
Modell: Präambelsymbole generieren -> IFFT -> CP einfügen -> Kanal -> Synchronisation ( -> Entzerrung)
End of explanation
pn1 = np.random.choice([-np.sqrt(2), np.sqrt(2)], K//2) # erste PN Sequenz
pn2 = np.random.choice([-1, 1], K)  # zweite PN Sequenz
Explanation: Präambel generieren (BPSK)
End of explanation
c1 = np.zeros((L,))
c1[L//2-K//2:L//2+K//2:2] = pn1  # nur gerade Unterträger werden belegt
s1 = np.fft.ifft(np.fft.fftshift(c1))
c2 = np.zeros((L,))
c2[L//2-K//2:L//2+K//2] = pn2  # alle Unterträger werden belegt
s2 = np.fft.ifft(np.fft.fftshift(c2))
v = np.sqrt(2)*pn2[::2]/pn1
Explanation: IFFT
End of explanation
s = np.concatenate((s1[-n_CP:], s1, s2[-n_CP:], s2))
f = np.linspace(-0.5, 0.5, L, endpoint=False)
S1 = np.abs(np.fft.fftshift(np.fft.fft(s[:L])))**2
plt.plot(f, np.abs(S1))
plt.title("|S1|"); plt.xlabel("Normalisierte Frequenz"); plt.xlim([-0.5, 0.5])
plt.show()
S2 = np.abs(np.fft.fftshift(np.fft.fft(s[L+2*n_CP:])))**2
plt.plot(f, np.abs(S2))
plt.title("|S2|"); plt.xlabel("Normalisierte Frequenz"); plt.xlim([-0.5, 0.5])
plt.show()
plt.plot(np.real(s[n_CP:n_CP+L]))
plt.title("Erstes Synchronisationssymbol (ohne CP, nur Realteil)"); 
plt.ylabel("Amplitude"); plt.xlabel("Sample Index");plt.xlim([0, L])
plt.show()
Explanation: Cyclic Prefix einfügen
End of explanation
# Faltung mit Kanalimpulsantwort
r = np.convolve(s, h, mode="full")
filter_delay = (len(h)-1)//2
r = r[filter_delay:]  # Einschwingzeit entfernen
# Frequenzoffset
r *= np.exp(1j*2*np.pi*f_off*np.arange(len(r))/L)
# Zeitoffset
r = np.concatenate((np.zeros((t_off,), dtype=complex), r))
# AWGN
P_s = np.mean(np.abs(r[t_off:])**2)  # mittlere Signalleistung
sigma = 10**(-SNR/20) * np.sqrt(P_s)
n = sigma * (np.sqrt(1/2) * np.random.randn(len(r)) + 1j * np.sqrt(1/2) * np.random.randn(len(r)))
r += n
Explanation: Kanal
End of explanation
# Fixed Lag Correlation
Pd = np.asarray([np.sum(np.conj(r[i:i+L//2])*r[i+L//2:i+L]) for i in range(len(r)-L)])
Rd = np.asarray([np.sum(np.abs(r[i+L//2:i+L])**2) for i in range(len(r) - L)])
M = np.abs(Pd/Rd)**2
plt.plot(np.abs(M)); plt.title("M(d)"); plt.xlabel("Verschiebung"); plt.xlim([0, len(M)]); plt.ylim([0,1.1]); plt.show()
plt.plot(np.angle(Pd)); plt.title("arg(P(d))"); plt.xlabel("Verschiebung"); plt.xlim([0, len(M)]); plt.show()
# Signaldetektion
signal_detected = any(M > threshold)
print("Präambel erkannt!" if signal_detected else "Keine Präambel erkannt.")
# Schätzung Zeitoffset
if signal_detected:
    # t_est = np.argmax(M)  # keine gute Idee wegen des Plateaus
    M_over_threshold = [idx for idx in range(len(M)) if M[idx] > threshold]
    t_est = M_over_threshold[0] + 3*n_CP//4  # "ad-hoc" Ansatz
    t_valid = [t_off+len(h), t_off+n_CP]  # Im realen System ist t_off natürlich unbekannt!
    t_est_valid = t_valid[0] <= t_est < t_valid[1]
    print("Geschätzter Zeitoffset: " + str(t_est) + " Samples. Gültiger Bereich: " + str(t_valid) + ". Schätzung im gültigen Bereich: " + str(t_est_valid))
# Fractional Frequency Offset Schätzung und Korrektur
if signal_detected:
    f_F = np.angle(Pd[t_est])/np.pi
    r *= np.exp(-1j*2*np.pi*f_F*np.arange(len(r))/L)
    print("Geschätzter Fractional Frequency Offset: " + str(f_F))
# Integer Frequency Offset Schätzung und Korrektur
if signal_detected:
    # Präambelsymbole extrahieren, Cyclic Prefix wegschneiden
    s1_est = r[t_est:t_est+L]  
    s2_est = r[t_est+L+n_CP:t_est+2*L+n_CP]
    
    # in den Frequenzbereich transformieren, in intuitive Lage bringen
    x1 = np.fft.fftshift(np.fft.fft(s1_est))  
    x2 = np.fft.fftshift(np.fft.fft(s2_est)) 
    
    # B(g) berechnen
    k = np.concatenate((np.arange(-K//2, 0, 2), np.arange(2, K//2, 2))) + L//2
    v_nodc = np.delete(v, K//2//2)
    E_x2 = np.sum(np.abs(x2[k])**2)
    possible_offsets = range(-(L-K)//2//2, (L-K)//2//2, 1) 
    B = np.zeros((len(possible_offsets,)), dtype=np.complex)
    for i, g in enumerate(possible_offsets):
        B[i] = (np.abs(np.sum(np.conj(x1[k+2*g])*np.conj(v_nodc)*x2[k+2*g]))**2)/(2*E_x2**2)
    # Integer Frequency Offset berechnen
    g_est = np.argmax(B)
    f_I = 2*g_est - (L-K)//2
    print("Geschätzter Integer Frequency Offset: " + str(f_I))
    
    print("Geschätzter Frequenzoffset (integer + fractional): " + str(f_I + f_F))
    plt.plot(possible_offsets, np.abs(B)); plt.title("|B(g)|"); plt.xlabel("g"); plt.xlim([possible_offsets[0], possible_offsets[-1]]); plt.show()    
    plt.plot(f, np.abs(x1)); plt.title("|x1|"); plt.xlabel("Normalisierte Frequenz"); plt.xlim([f[0], f[-1]]); plt.show()
    plt.plot(f, np.abs(x2)); plt.title("|x2|"); plt.xlabel("Normalisierte Frequenz"); plt.xlim([f[0], f[-1]]); plt.show()
Explanation: Synchronisation
End of explanation
if signal_detected:
    band_start_idx = L//2 - K//2 + f_I
    band_stop_idx = band_start_idx + K   
    f_I_valid = band_start_idx >= 0 and band_stop_idx < L
    if(not f_I_valid):
        print("Ungültiger Frequenzoffset!")
    else:
        c1_est = x1[band_start_idx : band_stop_idx]
        c2_est = x2[band_start_idx : band_stop_idx]
        plt.scatter(c1_est.real, c1_est.imag); plt.title("c1 geschätzt, ohne Entzerrung"); plt.xlabel("I"); plt.ylabel("Q"); plt.show()
        plt.scatter(c2_est.real, c2_est.imag); plt.title("c2 geschätzt, ohne Entzerrung"); plt.xlabel("I"); plt.ylabel("Q"); plt.show()
        ZF_eq = 1 / (c2_est / pn2)  # berechne Zero-Forcing Entzerrerkoeffizienten um den verbleibenden Phasenoffset zu korrigieren
        c1_est_eq = c1_est * ZF_eq
        c2_est_eq = c2_est * ZF_eq
        plt.scatter(ZF_eq.real, ZF_eq.imag); plt.title("Entzerrerkoeffizienten")
        plt.xlabel("I"); plt.ylabel("Q"); plt.show()
        plt.plot(np.abs(ZF_eq)); plt.title("Betrag der Entzerrerkoeffizienten")
        plt.xlabel("Index"); plt.show()
        plt.scatter(c1_est_eq.real, c1_est_eq.imag); plt.title("c1 geschätzt, mit Entzerrung")
        plt.xlabel("I"); plt.ylabel("Q"); plt.show()
        plt.scatter(c2_est_eq.real, c2_est_eq.imag); plt.title("c2 geschätzt, mit Entzerrung")
        plt.xlabel("I"); plt.ylabel("Q"); plt.show()    
Explanation: Entzerrung
Anmerkung: Entzerrung ist nicht mehr Teil des beschriebenen Verfahrens von Schmidl und Cox, es kann jedoch auch dafür verwendet werden.
End of explanation
<END_TASK>