+# %%
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy as stats
+from scipy.fftpack import fft, rfft, irfft, fftfreq, rfftfreq
+import scipy.signal
+from scipy.stats import norm, expon
+import math
+
+def diff_letters(a,b):
+ return sum ( a[i] != b[i] for i in range(len(a)) )
+
+
+# %%
+M = 8
+first = (2 - 1 - M)
+last = (16 - 1 - M)
+PAM8 = np.linspace(first, last, 8)
+
+M = 4
+QPSK = [[math.cos(2* math.pi * m / M), math.sin(2 * math.pi * m / M)] for m in range(M)]
+
+M = 8
+PSK8 = [[math.cos(2 * math.pi * m / M), math.sin(2 * math.pi * m / M)] for m in range(M)]
+
+M = 4
+first = (-3)
+last = (3)
+PAM4 = np.linspace(first, last, 4)
+QAM16 = []
+for i in PAM4:
+ for j in PAM4:
+ QAM16.append([i, j])
+
+
+# %%
+# Assign gray codes:
+
+g_PAM8 = {format(k, '03b'): i for k, i in zip([0, 1, 3, 2, 6, 7, 5, 4], PAM8)}
+
+g_QPSK = {format(k, '02b'): i for k, i in zip([0,1,3,2], QPSK)}
+
+g_PSK8 = {format(k, '03b'): i for k, i in zip([0, 1, 3, 2, 6, 7, 5, 4], PSK8)}
+
+g_QAM16 = {format(k, '04b'): i for k, i in zip([0, 1, 3, 2, 4, 5, 7, 6, 12, 13, 15, 14, 8, 9, 11, 10], QAM16)}
+
+# %%
+# PAM-8
+plt.figure(figsize=(4, 1))
+plt.title("8-PAM modulation, gray coding")
+plt.plot(PAM8, [0]*8, 'o')
+plt.yticks([])
+plt.xticks(PAM8)
+plt.grid(True)
+for (k, s) in g_PAM8.items():
+ plt.text(s, 0, k)
+
+# QPSK
+x = [x for x, y in QPSK]
+y = [y for x, y in QPSK]
+
+circle_x = np.linspace(0, 2*np.pi, 100)
+
+plt.figure(figsize=(4, 4))
+plt.title("QPSK modulation, gray coding")
+plt.plot(np.cos(circle_x), np.sin(circle_x), color="orange")
+plt.plot(x,y,'o', color="blue")
+plt.xticks(x)
+plt.yticks(y)
+plt.grid(True)
+for (k, s) in g_QPSK.items():
+ plt.text(*s, k)
+
+# PSK8
+x = [x for x, y in PSK8]
+y = [y for x, y in PSK8]
+
+plt.figure(figsize=(4, 4))
+plt.title("8-PSK modulation, gray coding")
+plt.plot(np.cos(circle_x), np.sin(circle_x), color="orange")
+plt.plot(x,y,'o', color="blue")
+plt.xticks(x, rotation=45)
+plt.yticks(y)
+plt.grid(True)
+for (k, s) in g_PSK8.items():
+ plt.text(*s, k)
+
+# QAM16
+x = [x for x, y in QAM16]
+y = [y for x, y in QAM16]
+
+plt.figure(figsize=(4, 4))
+plt.title("16-QAM")
+plt.plot(x,y,'o')
+plt.xticks(x)
+plt.yticks(y)
+plt.grid(True)
+for (k, s) in g_QAM16.items():
+ plt.text(*s, k)
+
+
+# %%
+# Plot definitions
+num_symbols = 10**6
+SNR_list = np.linspace(-20, 20, 10)
+
+
+# %%
+# PAM8
+
+symbol_error_rate = []
+bit_error_rate = []
+
+symbols = np.random.randint(0, len(PAM8), size=(num_symbols))
+
+PAM8_symbols = np.zeros(num_symbols, dtype=int)
+for i, s in enumerate(symbols):
+ PAM8_symbols[i] = PAM8[s]
+
+for SNR_dB in SNR_list:
+
+ SNR_linear = 10**(SNR_dB / 10)
+
+ signal_power = np.mean(np.abs(PAM8_symbols)**2)
+ noise_power = signal_power / SNR_linear
+
+ noise = math.sqrt(noise_power) * np.random.normal(0, 1, len(PAM8_symbols))
+ rx_signal = PAM8_symbols + noise
+
+ distances = np.zeros(len(PAM8), dtype=float)
+ detected_symbols = np.zeros(num_symbols, dtype=int)
+ for i, r in enumerate(rx_signal):
+ for j, s in enumerate(PAM8):
+ distances[j] = np.linalg.norm(r - s)
+
+ detected_symbols[i] = PAM8[np.argmin(distances)]
+
+ errors = 0
+ bit_errors = 0
+ for i in range(num_symbols):
+ if PAM8_symbols[i] != detected_symbols[i]:
+ errors += 1
+ code = ""
+ expected_code = ""
+ for k, v in g_PAM8.items():
+ if v == PAM8_symbols[i]:
+ expected_code = k
+ elif v == detected_symbols[i]:
+ code = k
+
+ if code != "" and expected_code != "":
+ break
+
+ bit_errors += diff_letters(code, expected_code)
+
+ SER = errors / num_symbols
+
+ BER = bit_errors / (num_symbols * 3) # 3 bits per num_symbols
+
+ # Make the plot for specific SNR_dB
+ if SNR_dB == 30:
+ plt.figure(figsize=(4, 1))
+ plt.title("Recieved symbols, 8-PAM modulation (SNR_dB = 30)")
+ plt.plot(rx_signal, np.zeros(num_symbols), 'o', color="red")
+ plt.plot(PAM8, [0]*len(PAM8), 'o', color="blue")
+ plt.yticks([])
+ plt.xticks(PAM8)
+ plt.grid(True)
+ for (k, s) in g_PAM8.items():
+ plt.text(s, 0, k)
+
+
+ symbol_error_rate.append(SER)
+ bit_error_rate.append(BER)
+
+# Final curve plot
+plt.figure(figsize=(4, 4))
+plt.title("8-PAM modulation error rate")
+plt.semilogy(SNR_list, symbol_error_rate, 'o-', color="blue")
+plt.semilogy(SNR_list, bit_error_rate, 'o-', color="red")
+plt.xlabel("SNR (dB)")
+plt.ylabel("Error rate")
+plt.legend(["Symbol error rate", "Bit error rate"])
+
+
+# %%
+# QPSK
+
+symbol_error_rate = []
+bit_error_rate = []
+
+symbols = np.random.randint(0, len(QPSK), size=(num_symbols))
+
+QPSK_symbols = np.zeros((num_symbols, 2), dtype=float)
+for i, s in enumerate(symbols):
+ QPSK_symbols[i] = QPSK[s]
+
+for SNR_dB in SNR_list:
+
+ SNR_linear = 10**(SNR_dB / 10)
+
+ signal_power = np.mean(np.abs(QPSK_symbols)**2)
+ noise_power = signal_power / SNR_linear
+
+ noise = math.sqrt(noise_power) * np.random.normal(0, 1, size=(num_symbols, 2))
+ rx_signal = QPSK_symbols + noise
+
+ distances = np.zeros(len(QPSK), dtype=float)
+ detected_symbols = np.zeros((num_symbols, 2), dtype=float)
+ for i, r in enumerate(rx_signal):
+ for j, s in enumerate(QPSK):
+ distances[j] = np.linalg.norm(r - s)
+
+ detected_symbols[i] = QPSK[np.argmin(distances)]
+
+ errors = 0
+ bit_errors = 0
+ for i in range(num_symbols):
+ if QPSK_symbols[i][0] != detected_symbols[i][0] or QPSK_symbols[i][1] != detected_symbols[i][1]:
+ errors += 1
+ code = ""
+ expected_code = ""
+ for k, v in g_QPSK.items():
+ if v[0] == QPSK_symbols[i][0] and v[1] == QPSK_symbols[i][1]:
+ expected_code = k
+ elif v[0] == detected_symbols[i][0] and v[1] == detected_symbols[i][1]:
+ code = k
+
+ if code != "" and expected_code != "":
+ break
+
+ bit_errors += diff_letters(code, expected_code)
+
+ SER = errors / num_symbols
+
+ BER = bit_errors / (num_symbols * 2) # 2 bits per symbol
+
+ # Make the plot for specific SNR_dB
+ if SNR_dB == 20:
+ r_x = [x for x, y in rx_signal]
+ r_y = [y for x, y in rx_signal]
+
+ x = [x for x, y in QPSK]
+ y = [y for x, y in QPSK]
+
+ plt.figure(figsize=(4, 4))
+ plt.title("Recieved symbols, QPSK modulation (SNR_dB = 20)")
+ plt.plot(np.cos(circle_x), np.sin(circle_x), color="orange")
+ plt.plot(r_x, r_y, 'o', color="red")
+ plt.plot(x,y,'o', color="blue")
+ plt.xticks(x, rotation=45)
+ plt.yticks(y)
+ plt.grid(True)
+ for (k, s) in g_QPSK.items():
+ plt.text(*s, k)
+
+ symbol_error_rate.append(SER)
+ bit_error_rate.append(BER)
+
+# Final curve plot
+plt.figure(figsize=(4, 4))
+plt.title("QPSK modulation error rate")
+plt.semilogy(SNR_list, symbol_error_rate, 'o-', color="blue")
+plt.semilogy(SNR_list, bit_error_rate, 'o-', color="red")
+plt.legend(["Symbol error rate", "Bit error rate"])
+plt.xlabel("SNR (dB)")
+plt.ylabel("Error rate")
+
+# %%
+# PSK8
+
+symbol_error_rate = []
+bit_error_rate = []
+
+symbols = np.random.randint(0, len(PSK8), size=(num_symbols))
+
+PSK8_symbols = np.zeros((num_symbols, 2), dtype=float)
+for i, s in enumerate(symbols):
+ PSK8_symbols[i] = PSK8[s]
+
+for SNR_dB in SNR_list:
+
+ SNR_linear = 10**(SNR_dB / 10)
+
+ signal_power = np.mean(np.abs(PSK8_symbols)**2)
+ noise_power = signal_power / SNR_linear
+
+ noise = math.sqrt(noise_power) * np.random.normal(0, 1, size=(num_symbols, 2))
+ rx_signal = PSK8_symbols + noise
+
+ distances = np.zeros(len(PSK8), dtype=float)
+ detected_symbols = np.zeros((num_symbols, 2), dtype=float)
+ for i, r in enumerate(rx_signal):
+ for j, s in enumerate(PSK8):
+ distances[j] = np.linalg.norm(r - s)
+
+ detected_symbols[i] = PSK8[np.argmin(distances)]
+
+ errors = 0
+ bit_errors = 0
+ for i in range(num_symbols):
+ if PSK8_symbols[i][0] != detected_symbols[i][0] or PSK8_symbols[i][1] != detected_symbols[i][1]:
+ errors += 1
+ code = ""
+ expected_code = ""
+ for k, v in g_PSK8.items():
+ if v[0] == PSK8_symbols[i][0] and v[1] == PSK8_symbols[i][1]:
+ expected_code = k
+ elif v[0] == detected_symbols[i][0] and v[1] == detected_symbols[i][1]:
+ code = k
+
+ if code != "" and expected_code != "":
+ break
+
+ bit_errors += diff_letters(code, expected_code)
+
+ SER = errors / num_symbols
+
+ BER = bit_errors / (num_symbols * 3) # 3 bits per symbol
+
+ # Make the plot for specific SNR_dB
+ if SNR_dB == 20:
+ r_x = [x for x, y in rx_signal]
+ r_y = [y for x, y in rx_signal]
+
+ x = [x for x, y in PSK8]
+ y = [y for x, y in PSK8]
+
+ plt.figure(figsize=(4, 4))
+ plt.title("Recieved symbols, 8-PSK modulation (SNR_dB = 20)")
+ plt.plot(np.cos(circle_x), np.sin(circle_x), color="orange")
+ plt.plot(r_x, r_y, 'o', color="red")
+ plt.plot(x,y,'o', color="blue")
+ plt.xticks(x, rotation=45)
+ plt.yticks(y)
+ plt.grid(True)
+ for (k, s) in g_PSK8.items():
+ plt.text(*s, k)
+
+ symbol_error_rate.append(SER)
+ bit_error_rate.append(BER)
+
+# Final curve plot
+plt.figure(figsize=(4, 4))
+plt.title("8-PSK modulation error rate")
+plt.semilogy(SNR_list, symbol_error_rate, 'o-', color="blue")
+plt.semilogy(SNR_list, bit_error_rate, 'o-', color="red")
+plt.legend(["Symbol error rate", "Bit error rate"])
+plt.xlabel("SNR (dB)")
+plt.ylabel("Error rate")
+
+
+# %%
+# QAM16
+
+symbol_error_rate = []
+bit_error_rate = []
+
+symbols = np.random.randint(0, len(QAM16), size=(num_symbols))
+
+QAM16_symbols = np.zeros((num_symbols, 2), dtype=float)
+for i, s in enumerate(symbols):
+ QAM16_symbols[i] = QAM16[s]
+
+for SNR_dB in SNR_list:
+
+ SNR_linear = 10**(SNR_dB / 10)
+
+ signal_power = np.mean(np.abs(QAM16_symbols)**2)
+ noise_power = signal_power / SNR_linear
+
+ noise = math.sqrt(noise_power) * np.random.normal(0, 1, size=(num_symbols, 2))
+ rx_signal = QAM16_symbols + noise
+
+ distances = np.zeros(len(QAM16), dtype=float)
+ detected_symbols = np.zeros((num_symbols, 2), dtype=float)
+ for i, r in enumerate(rx_signal):
+ for j, s in enumerate(QAM16):
+ distances[j] = np.linalg.norm(r - s)
+
+ detected_symbols[i] = QAM16[np.argmin(distances)]
+
+ errors = 0
+ bit_errors = 0
+ for i in range(num_symbols):
+ if QAM16_symbols[i][0] != detected_symbols[i][0] or QAM16_symbols[i][1] != detected_symbols[i][1]:
+ errors += 1
+ code = ""
+ expected_code = ""
+ for k, v in g_QAM16.items():
+ if v[0] == QAM16_symbols[i][0] and v[1] == QAM16_symbols[i][1]:
+ expected_code = k
+ elif v[0] == detected_symbols[i][0] and v[1] == detected_symbols[i][1]:
+ code = k
+
+ if code != "" and expected_code != "":
+ break
+
+ bit_errors += diff_letters(code, expected_code)
+
+ SER = errors / num_symbols
+
+ BER = bit_errors / (num_symbols * 4) # 4 bits per symbol
+
+ # Make the plot for specific SNR_dB
+ if SNR_dB == 20:
+ r_x = [x for x, y in rx_signal]
+ r_y = [y for x, y in rx_signal]
+
+ x = [x for x, y in QAM16]
+ y = [y for x, y in QAM16]
+
+ plt.figure(figsize=(4, 4))
+ plt.title("Recieved symbols, 16-QAM modulation (SNR_dB = 20)")
+ plt.plot(r_x, r_y, 'o', color="red")
+ plt.plot(x,y,'o', color="blue")
+ plt.xticks(x)
+ plt.yticks(y)
+ plt.grid(True)
+ for (k, s) in g_QAM16.items():
+ plt.text(*s, k)
+
+ symbol_error_rate.append(SER)
+ bit_error_rate.append(BER)
+
+# Final curve plot
+plt.figure(figsize=(4, 4))
+plt.title("16-QAM modulation error rate")
+plt.semilogy(SNR_list, symbol_error_rate, 'o-', color="blue")
+plt.semilogy(SNR_list, bit_error_rate, 'o-', color="red")
+plt.legend(["Symbol error rate", "Bit error rate"])
+plt.xlabel("SNR (dB)")
+plt.ylabel("Error rate")
+
+
+# %%
+# Hamming code
+
+def hamming_encode(data):
+ if len(data) % 4 != 0:
+ raise ValueError("input data multiple of 4")
+
+ G = [[1, 0, 0, 1, 0, 1, 1],
+ [0, 1, 0, 1, 0, 1, 0],
+ [0, 0, 1, 1, 0, 0, 1],
+ [0, 0, 0, 0, 1, 1, 1]]
+
+ encoded_data = []
+ for k in range(0, len(data), 4):
+ block = data[k:k+4]
+ for i in range(7):
+ sum = 0
+ for j in range(4):
+ sum += block[j] * G[j][i]
+ encoded_data.append(sum % 2)
+
+ return encoded_data
+
+def hamming_decode(recieved_data):
+ if len(recieved_data) % 7 != 0:
+ raise ValueError("Recieved data length mult of 7")
+
+ H = [[1, 0, 1, 0, 1, 0, 1],
+ [1, 1, 0, 0, 1, 1, 0],
+ [1, 1, 1, 1, 0, 0, 0]]
+
+ decoded_data = []
+ for k in range(0, len(recieved_data), 7):
+ block = recieved_data[k:k+7]
+ syndrome= []
+ for i in range(3):
+ syndrome_sum = 0
+ for j in range(7):
+ syndrome_sum += block[j] * H[i][j]
+ syndrome.append(syndrome_sum % 2)
+ error_pos = sum([2 ** i for i, bit in enumerate(syndrome) if bit])
+ if error_pos > 0:
+ block[7-error_pos] ^= 1
+ decoded_data.extend([block[0], block[1], block[2], block[4]])
+ return decoded_data
+
+num_data_bits_by_4 = 10**6
+
+data = np.zeros(4 * num_data_bits_by_4)
+for i in range(4 * num_data_bits_by_4):
+ data[i] = np.random.randint(0, 2)
+
+encoded = np.array(hamming_encode(data))
+
+ber_encoded = []
+ber_plain = []
+
+QPSK_symbols_encoded = np.array([g_QPSK[str(int(block[0])) + str(int(block[1]))] for block in encoded.reshape(int(len(encoded) / 2), 2)])
+QPSK_symbols_plain = np.array([g_QPSK[str(int(block[0])) + str(int(block[1]))] for block in data.reshape(int(len(data) / 2), 2)])
+
+num_symbols_encoded = len(QPSK_symbols_encoded)
+num_symbols_plain = len(QPSK_symbols_plain)
+
+SNR_list = np.linspace(-20, 20, 10)
+
+# Encoded QPSK
+for SNR_dB in SNR_list:
+
+ SNR_linear = 10**(SNR_dB / 10)
+
+ signal_power = np.mean(np.abs(QPSK_symbols_encoded)**2)
+ noise_power = signal_power / SNR_linear
+
+ noise = math.sqrt(noise_power) * np.random.normal(0, 1, size=(num_symbols_encoded, 2))
+ rx_signal = QPSK_symbols_encoded + noise
+
+ distances = np.zeros(len(QPSK), dtype=float)
+ detected_symbols = np.zeros((num_symbols_encoded, 2), dtype=float)
+ for i, r in enumerate(rx_signal):
+ for j, s in enumerate(QPSK):
+ distances[j] = np.linalg.norm(r - s)
+
+ detected_symbols[i] = QPSK[np.argmin(distances)]
+
+ recv = []
+ for i in range(num_symbols_encoded):
+ code = ""
+ for k, v in g_QPSK.items():
+ if v[0] == detected_symbols[i][0] and v[1] == detected_symbols[i][1]:
+ code = k
+ break
+
+ for char in code:
+ recv.append(int(char))
+
+ decoded = np.array(hamming_decode(recv))
+
+ bit_errors = np.sum(data != decoded)
+ BER = bit_errors / (num_symbols_encoded * 2) # 2 bits per symbol
+
+ ber_encoded.append(BER)
+
+# Plain QPSK
+for SNR_dB in SNR_list:
+
+ SNR_linear = 10**(SNR_dB / 10)
+
+ signal_power = np.mean(np.abs(QPSK_symbols_plain)**2)
+ noise_power = signal_power / SNR_linear
+
+ noise = math.sqrt(noise_power) * np.random.normal(0, 1, size=(num_symbols_plain, 2))
+ rx_signal = QPSK_symbols_plain + noise
+
+ distances = np.zeros(len(QPSK), dtype=float)
+ detected_symbols = np.zeros((num_symbols_plain, 2), dtype=float)
+ for i, r in enumerate(rx_signal):
+ for j, s in enumerate(QPSK):
+ distances[j] = np.linalg.norm(r - s)
+
+ detected_symbols[i] = QPSK[np.argmin(distances)]
+
+ recv = []
+ for i in range(num_symbols_plain):
+ code = ""
+ for k, v in g_QPSK.items():
+ if v[0] == detected_symbols[i][0] and v[1] == detected_symbols[i][1]:
+ code = k
+ break
+
+ for char in code:
+ recv.append(int(char))
+
+ bit_errors = np.sum(data != np.array(recv))
+ BER = bit_errors / (num_symbols_plain * 2) # 2 bits per symbol
+
+ ber_plain.append(BER)
+
+plt.figure(3)
+plt.semilogy(SNR_list, ber_encoded, 'o-', label="With hamming", color="blue")
+plt.semilogy(SNR_list, ber_plain, 'o-', label="Without encoding", color="red")
+plt.title("Bit error rate for QPSK communication with and without encoding")
+plt.xlabel("SNR (dB)")
+plt.ylabel("Bit error rate")
+plt.legend()
+plt.grid(True)
+plt.show()
projects / TSKS35.git / commitdiff
