rewrite
This commit is contained in:
5
.vscode/settings.json
vendored
Normal file
5
.vscode/settings.json
vendored
Normal file
@@ -0,0 +1,5 @@
|
|||||||
|
{
|
||||||
|
"python.analysis.extraPaths": [
|
||||||
|
".venv/**"
|
||||||
|
]
|
||||||
|
}
|
||||||
@@ -1,3 +1,7 @@
|
|||||||
# ieee802_constellations
|
# ieee802_constellations
|
||||||
|
|
||||||
|
## references
|
||||||
|
|
||||||
|
- <https://www.ieee802.org/3/bn/public/nov13/prodan_3bn_02_1113.pdf>
|
||||||
|
|
||||||
--created with create_repo.py
|
--created with create_repo.py
|
||||||
251
const.py
251
const.py
@@ -1,251 +0,0 @@
|
|||||||
from functools import cache
|
|
||||||
from typing import Tuple
|
|
||||||
import numpy as np
|
|
||||||
import const_utils
|
|
||||||
import matplotlib.pyplot as plt
|
|
||||||
|
|
||||||
class ConstellationPoints():
|
|
||||||
|
|
||||||
def __init__(self, length:int=None, constellation_dict:dict=None, radius:Tuple[int, float]=1):
|
|
||||||
self._radius = radius
|
|
||||||
if constellation_dict:
|
|
||||||
self._constellation = constellation_dict
|
|
||||||
elif length:
|
|
||||||
self._length = length
|
|
||||||
|
|
||||||
|
|
||||||
def _generate_from_length(self, length):
|
|
||||||
const_utils.validate_intpow2(length, 'length')
|
|
||||||
object.__setattr__(self, '_length', length)
|
|
||||||
object.__setattr__(self, '_n', int(self._length/2+0.5))
|
|
||||||
object.__setattr__(self, '_m', self._length // 2)
|
|
||||||
object.__setattr__(self, '_constellation', const_utils.bidict(generate_rectangular_constellation(self._length)))
|
|
||||||
|
|
||||||
|
|
||||||
def _generate_from_dict(self, constellation_dict):
|
|
||||||
const_utils.validate_coords(constellation_dict)
|
|
||||||
const_utils.validate_int(next(iter(constellation_dict.keys())), 'labels')
|
|
||||||
|
|
||||||
# check if constellation is a one-to-one mapping
|
|
||||||
if len(set(constellation_dict.values())) != len(constellation_dict):
|
|
||||||
raise ValueError('constellation must be a one-to-one mapping of labels and coordinates')
|
|
||||||
|
|
||||||
object.__setattr__(self, '_constellation', const_utils.bidict(constellation_dict))
|
|
||||||
|
|
||||||
object.__setattr__(self, '_length', len(constellation_dict))
|
|
||||||
object.__setattr__(self, '_n', int(self._length/2+0.5))
|
|
||||||
object.__setattr__(self, '_m', self._length // 2)
|
|
||||||
|
|
||||||
|
|
||||||
def plot(self):
|
|
||||||
# STUB ConstellationPots.plot()
|
|
||||||
raise NotImplementedError()
|
|
||||||
|
|
||||||
|
|
||||||
def __len__(self):
|
|
||||||
return self._length
|
|
||||||
|
|
||||||
|
|
||||||
@property
|
|
||||||
def symbols(self):
|
|
||||||
return dict(self._constellation).values()
|
|
||||||
|
|
||||||
|
|
||||||
@property
|
|
||||||
def labels(self):
|
|
||||||
return dict(self._constellation.inverse).values()
|
|
||||||
|
|
||||||
|
|
||||||
@property
|
|
||||||
def constellation(self):
|
|
||||||
return dict(self._constellation)
|
|
||||||
|
|
||||||
|
|
||||||
@property
|
|
||||||
def lookup(self):
|
|
||||||
return dict(self._constellation.inverse)
|
|
||||||
|
|
||||||
|
|
||||||
def get_symbol(self, label):
|
|
||||||
return self._constellation[label]
|
|
||||||
|
|
||||||
|
|
||||||
def get_label(self, symbol):
|
|
||||||
return self._constellation.inverse[symbol]
|
|
||||||
|
|
||||||
|
|
||||||
def __setattr__(self, name: str, value) -> None:
|
|
||||||
if name == '_constellation':
|
|
||||||
self._generate_from_dict(value)
|
|
||||||
elif name == '_length':
|
|
||||||
self._generate_from_length(value)
|
|
||||||
elif name == '_radius':
|
|
||||||
self._set_radius(value)
|
|
||||||
else:
|
|
||||||
object.__setattr__(self, name, value)
|
|
||||||
|
|
||||||
|
|
||||||
def _set_radius(self, value):
|
|
||||||
if isinstance(value, (int, float)):
|
|
||||||
object.__setattr__(self, '_radius',value)
|
|
||||||
else:
|
|
||||||
raise TypeError('radius must be an integer or a float')
|
|
||||||
|
|
||||||
|
|
||||||
@cache
|
|
||||||
def gray_1d(k: int, label: int) -> int:
|
|
||||||
const_utils.gray_1d_input_validation(k, label)
|
|
||||||
|
|
||||||
if k == 1:
|
|
||||||
return 1 if label==0 else -1
|
|
||||||
|
|
||||||
b0, new_symbol = const_utils.next_symbol(k, label)
|
|
||||||
return (1-2*b0)*(2**(k-1)+gray_1d(k-1, new_symbol))
|
|
||||||
|
|
||||||
|
|
||||||
def gray_2d(n: int, m: int, label: int) -> tuple:
|
|
||||||
const_utils.gray_2d_input_validation(n, m, label)
|
|
||||||
|
|
||||||
if m == 0:
|
|
||||||
return (gray_1d(n, label), 0) # it's a 1d case in disguise!
|
|
||||||
|
|
||||||
symbol_i, symbol_q = const_utils.split_symbol(n, m, label)
|
|
||||||
return (gray_1d(n, symbol_i), gray_1d(m, symbol_q))
|
|
||||||
|
|
||||||
|
|
||||||
def hamming_dist(a, b):
|
|
||||||
if not isinstance(a, int):
|
|
||||||
raise ValueError('a must be an integer')
|
|
||||||
if not isinstance(b, int):
|
|
||||||
raise ValueError('b must be an integer')
|
|
||||||
|
|
||||||
|
|
||||||
def euclidean_distance(coord1, coord2):
|
|
||||||
if isinstance(coord1, int):
|
|
||||||
return abs(coord1 - coord2)
|
|
||||||
return np.sqrt((coord1[0] - coord2[0])**2 + (coord1[1] - coord2[1])**2)
|
|
||||||
|
|
||||||
|
|
||||||
def find_nearest(coord, coords):
|
|
||||||
min_distance = float('inf')
|
|
||||||
nearest_symbols = []
|
|
||||||
|
|
||||||
for c in coords:
|
|
||||||
dist = euclidean_distance(coord, c)
|
|
||||||
if dist == 0:
|
|
||||||
continue
|
|
||||||
elif dist < min_distance:
|
|
||||||
min_distance = dist
|
|
||||||
nearest_symbols = [c]
|
|
||||||
elif dist == min_distance:
|
|
||||||
nearest_symbols.append(c)
|
|
||||||
|
|
||||||
return nearest_symbols
|
|
||||||
|
|
||||||
|
|
||||||
def gray_penalty(constellation):
|
|
||||||
raise NotImplementedError
|
|
||||||
# constellation: {label_0:coordinate_0, label_1:coordinate_1, .., label_2^n-1:coordinate_2^n-1}
|
|
||||||
|
|
||||||
# 2^n-QAM -> 2^n symbols S_i, where i=0,1,..2^n-1, ex. S_0 = (-3,-3) or S_0 = -2
|
|
||||||
# N(S_i): set of (euclidean) nearest symbols S_j -> N((-3,-3)) = {(-3,-2), (-3,-4), (-2,-3), (-4,-3)}
|
|
||||||
# |N(S_i)|: size of set N(S_i)
|
|
||||||
# l(S): label given by mapping -> inverse of gray_Qd -> generate all symbols/labels for given constellation
|
|
||||||
# wt(l_1, l_2), hamming distance btw. two labels
|
|
||||||
|
|
||||||
t = constellation['meta']['len']
|
|
||||||
inverted_constellation = {tuple(symbol):label for label,symbol in constellation.items() if label != 'meta'} # -> invert constellation dict
|
|
||||||
syms = [symbol for _, symbol in constellation.items()]
|
|
||||||
|
|
||||||
const_utils.validate_intpow2()
|
|
||||||
|
|
||||||
G = 0
|
|
||||||
for li, si in constellation.items():
|
|
||||||
N = find_nearest(si, syms)
|
|
||||||
size_N = len(N)
|
|
||||||
wt = sum(hamming_dist(inverted_constellation[tuple(sj)], li) for sj in N)
|
|
||||||
G += wt/size_N
|
|
||||||
G /= t
|
|
||||||
|
|
||||||
return G
|
|
||||||
|
|
||||||
|
|
||||||
def generate_rectangular_constellation(length: int):
|
|
||||||
# const_utils.validate_int(length, 'length')
|
|
||||||
const_utils.validate_intpow2(length, 'length')
|
|
||||||
lengthexp = int(np.log2(length))
|
|
||||||
n = int(lengthexp/2+0.5) # ceil
|
|
||||||
m = lengthexp // 2 # floor
|
|
||||||
return {label:gray_2d(n, m, label) for label in range(length)}
|
|
||||||
|
|
||||||
|
|
||||||
def transform_rectangular_mapping(constellation: ConstellationPoints):
|
|
||||||
n, m = constellation['meta'] # TODO def generate_rectangular_constellation(n)
|
|
||||||
# r, c = find_rows_columns(constellation)
|
|
||||||
|
|
||||||
# # example: 32-qam -> 2^(2n+1) -> n = 2
|
|
||||||
|
|
||||||
# two_n1 = np.log2(len(constellation))
|
|
||||||
# if int(two_n1) != two_n1:
|
|
||||||
# raise ValueError('only constellations with 2^m points allowed')
|
|
||||||
|
|
||||||
# if r == 1: # 1D-constellation
|
|
||||||
# return constellation
|
|
||||||
|
|
||||||
# # get n and m for one quadrant
|
|
||||||
# n = c/2
|
|
||||||
# m = r/2
|
|
||||||
|
|
||||||
# const_utils._validate_integer(n, 'n')
|
|
||||||
# const_utils._validate_integer(m, 'm')
|
|
||||||
|
|
||||||
# [ ] set transformed flag in constellation?
|
|
||||||
|
|
||||||
if n == m: # square 2^(2n)-QAM
|
|
||||||
return constellation
|
|
||||||
|
|
||||||
if n == 2 and m == 1: # rectangular 8-QAM (4*2)
|
|
||||||
return transform_8QAM(constellation)
|
|
||||||
|
|
||||||
elif n == m+2: # REVIEW m+2 correct? p. 7
|
|
||||||
new_const = {}
|
|
||||||
s = 2**(n-1)
|
|
||||||
for label, symbol in constellation.items():
|
|
||||||
# STUB transfrom_rectangular_mapping(constellation) -> generalized non-square-QAM
|
|
||||||
raise NotImplementedError()
|
|
||||||
|
|
||||||
else:
|
|
||||||
# TODO define what should happen here
|
|
||||||
return constellation
|
|
||||||
|
|
||||||
# 2^(2n+1)
|
|
||||||
# for 32-QAM: 2^5 -> n = 2
|
|
||||||
# rectangular grid of 4*8 -> 2*4 per quadrant
|
|
||||||
|
|
||||||
# for 128-QAM: 2^7 -> n = 3
|
|
||||||
# rectangular grid of
|
|
||||||
|
|
||||||
|
|
||||||
def transform_8QAM(constellation):
|
|
||||||
new_const = {}
|
|
||||||
for label, symbol in constellation.items():
|
|
||||||
if symbol[0] < 3:
|
|
||||||
new_const[label] = symbol
|
|
||||||
else:
|
|
||||||
i_rct, q_rct = symbol
|
|
||||||
i_cr = -np.sign(i_rct)*(4-np.abs(i_rct))
|
|
||||||
q_cr = np.sign(q_rct)*(np.abs(q_rct)+2)
|
|
||||||
new_const[label] = [i_cr, q_cr]
|
|
||||||
return new_const# rectangular 2^(m+n)-QAM
|
|
||||||
|
|
||||||
|
|
||||||
if __name__ == '__main__':
|
|
||||||
const0 = ConstellationPoints()
|
|
||||||
const128 = ConstellationPoints(length=128)
|
|
||||||
const_ext = ConstellationPoints(constellation_dict=generate_rectangular_constellation(64))
|
|
||||||
|
|
||||||
print(vars(const0))
|
|
||||||
print(vars(const128))
|
|
||||||
print(vars(const_ext))
|
|
||||||
|
|
||||||
# constellation_128 = {label:gray_2d(3, 4, label) for label in range(128)}
|
|
||||||
@@ -1,67 +0,0 @@
|
|||||||
import pytest
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
from const import gray_1d
|
|
||||||
|
|
||||||
# Happy path tests with various realistic test values
|
|
||||||
@pytest.mark.parametrize("k, symbol, expected", [
|
|
||||||
(0, 0, 0), # k=0 case
|
|
||||||
(1, 0, 1), # k=1, symbol=0 case
|
|
||||||
(1, 1, -1), # k=1, symbol=1 case
|
|
||||||
(2, 0, 3), # k=2, symbol=0 case
|
|
||||||
(2, 1, 1), # k=2, symbol=1 case
|
|
||||||
(2, 2, -3), # k=2, symbol=2 case
|
|
||||||
(2, 3, -1), # k=2, symbol=3 case
|
|
||||||
], ids=[
|
|
||||||
"k=0_symbol=0",
|
|
||||||
"k=1_symbol=0",
|
|
||||||
"k=1_symbol=1",
|
|
||||||
"k=2_symbol=0",
|
|
||||||
"k=2_symbol=1",
|
|
||||||
"k=2_symbol=2",
|
|
||||||
"k=2_symbol=3"
|
|
||||||
])
|
|
||||||
def test_gray_1d_happy_path(k, symbol, expected):
|
|
||||||
# Act
|
|
||||||
result = gray_1d(k, symbol)
|
|
||||||
|
|
||||||
# Assert
|
|
||||||
assert result == expected
|
|
||||||
|
|
||||||
# Various edge cases
|
|
||||||
@pytest.mark.parametrize("k, symbol, expected", [
|
|
||||||
(3, 0, 7), # k=3, symbol=0 case
|
|
||||||
(3, 7, -3), # k=3, symbol=7 case
|
|
||||||
(4, 0, 15), # k=4, symbol=0 case
|
|
||||||
(4, 15, -5), # k=4, symbol=15 case
|
|
||||||
], ids=[
|
|
||||||
"k=3_symbol=0",
|
|
||||||
"k=3_symbol=7",
|
|
||||||
"k=4_symbol=0",
|
|
||||||
"k=4_symbol=15"
|
|
||||||
])
|
|
||||||
def test_gray_1d_edge_cases(k, symbol, expected):
|
|
||||||
# Act
|
|
||||||
result = gray_1d(k, symbol)
|
|
||||||
|
|
||||||
# Assert
|
|
||||||
assert result == expected
|
|
||||||
|
|
||||||
# Various error cases
|
|
||||||
@pytest.mark.parametrize("k, symbol, exception", [
|
|
||||||
(-1, 0, ValueError), # k is negative
|
|
||||||
(1, -1, ValueError), # symbol is negative
|
|
||||||
(1, 2, ValueError), # symbol is out of range for k=1
|
|
||||||
(2, 4, ValueError), # symbol is out of range for k=2
|
|
||||||
], ids=[
|
|
||||||
"k_negative",
|
|
||||||
"symbol_negative",
|
|
||||||
"symbol_out_of_range_k1",
|
|
||||||
"symbol_out_of_range_k2"
|
|
||||||
])
|
|
||||||
def test_gray_1d_error_cases(k, symbol, exception):
|
|
||||||
# Act and Assert
|
|
||||||
with pytest.raises(exception):
|
|
||||||
gray_1d(k, symbol)
|
|
||||||
|
|
||||||
109
const_utils.py
109
const_utils.py
@@ -1,109 +0,0 @@
|
|||||||
import numpy as np
|
|
||||||
|
|
||||||
def validate_int(value, name):
|
|
||||||
if not isinstance(value, int):
|
|
||||||
raise ValueError(f'{name} must be an integer')
|
|
||||||
|
|
||||||
|
|
||||||
def validate_range(value, name, min_val=None, max_val=None):
|
|
||||||
if max_val is not None and value > max_val:
|
|
||||||
raise ValueError(f'{name} must be \u2265 {max_val}')
|
|
||||||
if min_val is not None and value < min_val:
|
|
||||||
raise ValueError(f'{name} must be \u2264 {min_val}')
|
|
||||||
|
|
||||||
|
|
||||||
def gray_1d_input_validation(k, symbol):
|
|
||||||
validate_int(symbol, 'symbol')
|
|
||||||
validate_int(k, 'k')
|
|
||||||
validate_range(k, 'k', min_val=1)
|
|
||||||
validate_range(symbol, 'symbol', min_val=0, max_val=2**k-1)
|
|
||||||
|
|
||||||
|
|
||||||
def next_symbol(k, symbol):
|
|
||||||
bits = format(symbol, 'b').zfill(k)
|
|
||||||
b0 = int(bits[0])
|
|
||||||
new_symbol = int(bits[1:], 2)
|
|
||||||
return b0,new_symbol
|
|
||||||
|
|
||||||
|
|
||||||
def gray_2d_input_validation(n, m, symbol):
|
|
||||||
validate_int(n, 'n')
|
|
||||||
validate_int(m, 'm')
|
|
||||||
validate_int(symbol, 'symbol')
|
|
||||||
validate_range(m, 'm', min_val=0)
|
|
||||||
validate_range(n, 'n', min_val=m)
|
|
||||||
validate_range(symbol, 'symbol', min_val=0, max_val=2**(m+n)-1)
|
|
||||||
|
|
||||||
|
|
||||||
def split_symbol(n, m, symbol):
|
|
||||||
bits = format(symbol, 'b').zfill(n+m)
|
|
||||||
symbol_i = int(bits[:n], 2)
|
|
||||||
symbol_q = int(bits[n:], 2)
|
|
||||||
return symbol_i,symbol_q
|
|
||||||
|
|
||||||
|
|
||||||
def validate_intpow2(value, name):
|
|
||||||
exponent = np.log2(value)
|
|
||||||
if exponent != int(exponent):
|
|
||||||
raise ValueError(f'{name} must be an integer power of 2')
|
|
||||||
|
|
||||||
|
|
||||||
def validate_coords(constellation_dict):
|
|
||||||
# TODO validate all coords, not only first
|
|
||||||
# bit of a hack, only looking at first element
|
|
||||||
if not isinstance(temp:=next(iter(constellation_dict.values())), tuple) or not isinstance(temp[0], int):
|
|
||||||
raise ValueError('coords must be tuples of integers')
|
|
||||||
|
|
||||||
|
|
||||||
# https://www.ieee802.org/3/bn/public/nov13/prodan_3bn_02_1113.pdf
|
|
||||||
|
|
||||||
class bidict(dict):
|
|
||||||
'''
|
|
||||||
### Summary:
|
|
||||||
A bidirectional dictionary (bidict) that allows bidirectional mapping between keys and values.
|
|
||||||
|
|
||||||
### Explanation:
|
|
||||||
This class extends the functionality of a standard dictionary to maintain a bidirectional mapping between keys and values. It provides methods to set, delete items, and retrieve the inverse mapping efficiently.
|
|
||||||
|
|
||||||
### Methods:
|
|
||||||
- `__init__(*args, **kwargs)`: Initializes the bidict with optional initial key-value pairs.
|
|
||||||
- `__setitem__(key, value)`: Sets a key-value pair in the bidict, updating the inverse mapping.
|
|
||||||
- `__delitem__(key)`: Deletes a key-value pair from the bidict and updates the inverse mapping.
|
|
||||||
- `__repr__()`: Returns a string representation of the bidict instance.
|
|
||||||
|
|
||||||
### Attributes:
|
|
||||||
- `inverse`: A dictionary that stores the inverse mapping of values to keys.
|
|
||||||
|
|
||||||
### Returns:
|
|
||||||
- No explicit return value for methods. The bidict instance is modified in place.
|
|
||||||
|
|
||||||
### Raises:
|
|
||||||
- No specific exceptions are raised by the methods in this class.
|
|
||||||
|
|
||||||
### Source:
|
|
||||||
bidict by user 'Basj' at https://stackoverflow.com/a/21894086 (CC BY-SA 4.0)
|
|
||||||
'''
|
|
||||||
|
|
||||||
def __init__(self, *args, **kwargs):
|
|
||||||
super(bidict, self).__init__(*args, **kwargs)
|
|
||||||
self.inverse = {}
|
|
||||||
for key, value in self.items():
|
|
||||||
self.inverse.setdefault(value, []).append(key)
|
|
||||||
|
|
||||||
|
|
||||||
def __setitem__(self, key, value):
|
|
||||||
if key in self:
|
|
||||||
self.inverse[self[key]].remove(key)
|
|
||||||
super(bidict, self).__setitem__(key, value)
|
|
||||||
self.inverse.setdefault(value, []).append(key)
|
|
||||||
|
|
||||||
|
|
||||||
def __delitem__(self, key):
|
|
||||||
self.inverse.setdefault(self[key], []).remove(key)
|
|
||||||
if self[key] in self.inverse and not self.inverse[self[key]]:
|
|
||||||
del self.inverse[self[key]]
|
|
||||||
super(bidict, self).__delitem__(key)
|
|
||||||
|
|
||||||
|
|
||||||
def __repr__(self):
|
|
||||||
return f'<bidict @[{id(self)}]'
|
|
||||||
375
gray.py
Normal file
375
gray.py
Normal file
@@ -0,0 +1,375 @@
|
|||||||
|
import numpy as np
|
||||||
|
import matplotlib.pyplot as plt
|
||||||
|
from typing import Self
|
||||||
|
from functools import cache
|
||||||
|
import sys
|
||||||
|
|
||||||
|
|
||||||
|
class constellation_point:
|
||||||
|
def __init__(self, num: int = 0, n: int = 0, m: int = 0):
|
||||||
|
self.number = num
|
||||||
|
self.n = n
|
||||||
|
self.m = m
|
||||||
|
self.symbol = self.encode(self.number)
|
||||||
|
self.constellation = self._gray(self.symbol, n, m)
|
||||||
|
self.transformed = self.transform()
|
||||||
|
self.neighbours = []
|
||||||
|
|
||||||
|
@staticmethod
|
||||||
|
@cache
|
||||||
|
def encode(number):
|
||||||
|
if number <= 1:
|
||||||
|
return number
|
||||||
|
symbol = constellation_point.encode(number - 1)
|
||||||
|
j = constellation_point.determine_bit_position(number)
|
||||||
|
symbol = symbol ^ (1 << j)
|
||||||
|
return symbol
|
||||||
|
|
||||||
|
@staticmethod
|
||||||
|
def determine_bit_position(i):
|
||||||
|
@cache
|
||||||
|
def find_lowest_bit(x):
|
||||||
|
ld = np.log2(x)
|
||||||
|
if ld == int(ld):
|
||||||
|
return ld
|
||||||
|
else:
|
||||||
|
return find_lowest_bit(x - 2 ** int(ld))
|
||||||
|
|
||||||
|
if i % 2:
|
||||||
|
j = 0
|
||||||
|
else:
|
||||||
|
j = find_lowest_bit(i)
|
||||||
|
return int(j)
|
||||||
|
|
||||||
|
@staticmethod
|
||||||
|
def _gray_1d(symbol: int, k: int):
|
||||||
|
if symbol >= 2**k or symbol < 0:
|
||||||
|
raise ValueError(f"symbol {symbol} out of range (0..{2**k})")
|
||||||
|
|
||||||
|
# symbol += 2**(k-1)
|
||||||
|
# symbol %= 2**k
|
||||||
|
|
||||||
|
if k < 0:
|
||||||
|
raise ValueError(f"k must >= 0, is {k}")
|
||||||
|
if k == 0:
|
||||||
|
return 0
|
||||||
|
else:
|
||||||
|
b0 = (symbol & (2 ** (k - 1))) >> (k - 1) # highest value bit, shifted to the right
|
||||||
|
b1k = symbol & ((2 ** (k - 1)) - 1) # remaining bits
|
||||||
|
return (1 - 2 * b0) * (2 ** (k - 1) + constellation_point._gray_1d(b1k, k=k - 1))
|
||||||
|
|
||||||
|
@staticmethod
|
||||||
|
def _gray_2d(symbol: int, n: int, m: int):
|
||||||
|
# a b
|
||||||
|
# symbol: [0 1 .. n-1][0 1 .. m-1]
|
||||||
|
# mask: 1 1 .. 1 0 0 0 0
|
||||||
|
maskb = 2 ** (m) - 1
|
||||||
|
b = symbol & maskb
|
||||||
|
maska = 2 ** (m + n) - 1 & ~maskb
|
||||||
|
a = (symbol & maska) >> m
|
||||||
|
return np.array((constellation_point._gray_1d(a, n), constellation_point._gray_1d(b, m)))
|
||||||
|
|
||||||
|
@staticmethod
|
||||||
|
def _gray(symbol: int, k: int = 0, m: int = 0):
|
||||||
|
if m == 0:
|
||||||
|
return np.array((constellation_point._gray_1d(symbol, k), 0)).flatten()
|
||||||
|
return constellation_point._gray_2d(symbol, k, m).flatten()
|
||||||
|
|
||||||
|
def set_neighbours(self, points: tuple[Self]):
|
||||||
|
self.neighbours = tuple(points)
|
||||||
|
self.hamming_distance = sum(
|
||||||
|
(self.symbol ^ neighbour.symbol).bit_count() for neighbour in self.neighbours
|
||||||
|
) / len(self.neighbours)
|
||||||
|
return self.hamming_distance
|
||||||
|
|
||||||
|
# @staticmethod
|
||||||
|
# def calc_hamming_distance(a:Self, b:Self):
|
||||||
|
# return (a.symbol ^ b.symbol).bit_count()
|
||||||
|
|
||||||
|
def transform(self):
|
||||||
|
if self.n == self.m:
|
||||||
|
return self.constellation
|
||||||
|
if self.n == 2:
|
||||||
|
return self._transform_8()
|
||||||
|
return self._transform_2n1()
|
||||||
|
|
||||||
|
def _transform_2n1(self):
|
||||||
|
ir, qr = self.constellation
|
||||||
|
# ic, qc = self.transformed
|
||||||
|
|
||||||
|
s = 2 ** (self.m - 1)
|
||||||
|
|
||||||
|
# 3, -4, 2 for n = 2, m = 1
|
||||||
|
if np.abs(ir) < 3 * s:
|
||||||
|
return self.constellation
|
||||||
|
|
||||||
|
if np.abs(qr) > s:
|
||||||
|
ic = np.sign(ir) * (np.abs(ir) - 2 * s)
|
||||||
|
qc = np.sign(qr) * (4 * s - np.abs(qr))
|
||||||
|
else:
|
||||||
|
ic = np.sign(ir) * (4 * s - np.abs(ir))
|
||||||
|
qc = np.sign(qr) * (np.abs(qr) + 2 * s)
|
||||||
|
|
||||||
|
return np.array((ic, qc))
|
||||||
|
|
||||||
|
def _transform_8(self):
|
||||||
|
ir, qr = self.constellation
|
||||||
|
|
||||||
|
# 3, -4, 2 for n = 2, m = 1
|
||||||
|
if ir < 3:
|
||||||
|
return self.constellation
|
||||||
|
|
||||||
|
ic = np.sign(ir) * (np.abs(ir) - 4)
|
||||||
|
qc = np.sign(qr) * (2 + np.abs(qr))
|
||||||
|
|
||||||
|
return np.array((ic, qc))
|
||||||
|
|
||||||
|
def __repr__(self):
|
||||||
|
_digits_num = int(np.log10(2 ** (self.n + self.m))) + 1
|
||||||
|
_digits_x = int(np.log10(2 ** (self.n) - 1)) + 1
|
||||||
|
_digits_y = int(np.log10(2 ** (self.m) - 1)) + 1
|
||||||
|
sym = f"{self.number:>{_digits_num}d}"
|
||||||
|
sym_b = f"{self.symbol:0{self.n + self.m}b}"
|
||||||
|
const = f"({self.constellation[0]:+0{_digits_x}d}, {self.constellation[1]:+0{_digits_y}d})"
|
||||||
|
return f"{sym} [{sym_b}] - {const} | n:{self.n:d}, m:{self.m:d}"
|
||||||
|
|
||||||
|
|
||||||
|
class constellation:
|
||||||
|
def __init__(self, n_symbols: int = 1):
|
||||||
|
self.k = int(np.ceil(np.log2(n_symbols)))
|
||||||
|
self.n_symbols = n_symbols
|
||||||
|
self.m = self.k // 2
|
||||||
|
self.n = self.k - self.m
|
||||||
|
|
||||||
|
# self.n, self.m = self.most_square_fact(n_symbols)
|
||||||
|
# self.n = int(np.ceil(np.log2(self.n)))
|
||||||
|
# self.m = int(np.ceil(np.log2(self.m)))
|
||||||
|
self.points = [constellation_point(i, self.n, self.m) for i in range(self.n_symbols)]
|
||||||
|
self.distances = self.calc_distances(transformed=False)
|
||||||
|
self.set_nearest_neighbours()
|
||||||
|
self.gray_penalty_original = self.calc_gray_penalty()
|
||||||
|
|
||||||
|
self.distances = self.calc_distances(transformed=True)
|
||||||
|
self.set_nearest_neighbours()
|
||||||
|
self.gray_penalty_transformed = self.calc_gray_penalty()
|
||||||
|
|
||||||
|
self.scale = 1
|
||||||
|
|
||||||
|
def calc_distances(self, transformed=True):
|
||||||
|
distances = np.zeros((len(self.points), len(self.points)))
|
||||||
|
for i in range(len(self.points)):
|
||||||
|
for j in range(len(self.points)):
|
||||||
|
distances[i, j] = self.calc_distance(self.points[i], self.points[j], transformed=transformed)
|
||||||
|
distances = np.ma.masked_array(distances, np.eye(*distances.shape))
|
||||||
|
return distances
|
||||||
|
|
||||||
|
def set_nearest_neighbours(self):
|
||||||
|
# neighbours = []
|
||||||
|
for i, point in enumerate(self.points):
|
||||||
|
dists = self.distances[i]
|
||||||
|
indices = np.where(dists == dists.min())[0]
|
||||||
|
point.set_neighbours([self.points[index] for index in indices])
|
||||||
|
# for j, point2 in enumerate(self.points):
|
||||||
|
# if i == j:
|
||||||
|
# continue
|
||||||
|
# neighbours.append(point2)
|
||||||
|
# point.set_neighbours(neighbours)
|
||||||
|
...
|
||||||
|
|
||||||
|
def calc_gray_penalty(self):
|
||||||
|
return sum(pt.hamming_distance for pt in self.points) / self.n_symbols
|
||||||
|
|
||||||
|
# def calc_gray_penalty(self):
|
||||||
|
# hamming = 0
|
||||||
|
# euclidean = 0
|
||||||
|
# for i, point in enumerate(self.points):
|
||||||
|
# for j, point2 in enumerate(self.points):
|
||||||
|
# if i == j:
|
||||||
|
# continue
|
||||||
|
# hamming += (point.symbol ^ point2.symbol).bit_count()
|
||||||
|
# euclidean += np.linalg.norm(point.transformed - point2.transformed)
|
||||||
|
# penalty = hamming/euclidean
|
||||||
|
# return penalty
|
||||||
|
|
||||||
|
def get_points_for_plot(self):
|
||||||
|
pts = [point.constellation for point in self.points]
|
||||||
|
xs_pre = [pt[0] * self.scale for pt in pts]
|
||||||
|
ys_pre = [pt[1] * self.scale for pt in pts]
|
||||||
|
labels = [point.symbol for point in self.points]
|
||||||
|
pts = [point.transformed for point in self.points]
|
||||||
|
xs = [pt[0] * self.scale for pt in pts]
|
||||||
|
ys = [pt[1] * self.scale for pt in pts]
|
||||||
|
return xs, ys, xs_pre, ys_pre, labels
|
||||||
|
|
||||||
|
def plot(self, show=True, annotate=False, old=False):
|
||||||
|
ax = plt.subplot()
|
||||||
|
xs, ys, xs_pre, ys_pre, labels = self.get_points_for_plot()
|
||||||
|
if old:
|
||||||
|
ax.scatter(xs_pre, ys_pre, marker="o", facecolors="none", edgecolors="b")
|
||||||
|
# circ = Circle.welzl([point.constellation for point in self.points])
|
||||||
|
# circ.plot(ax, show=False)
|
||||||
|
ax.scatter(xs, ys, color="b")
|
||||||
|
circ = Circle.welzl(list(zip(xs, ys)))
|
||||||
|
circ.plot(ax, show=False)
|
||||||
|
if annotate:
|
||||||
|
for i, label in enumerate(labels):
|
||||||
|
ax.annotate(
|
||||||
|
f"{label:0{self.n + self.m}b}",
|
||||||
|
xy=(xs[i], ys[i] + self.scale / 4),
|
||||||
|
horizontalalignment="center",
|
||||||
|
verticalalignment="bottom",
|
||||||
|
)
|
||||||
|
# if old:
|
||||||
|
# ax.annotate(
|
||||||
|
# f"{label:0{self.n+self.m}b}",
|
||||||
|
# xy=(xs_pre[i], ys_pre[i] + self.scale/4),
|
||||||
|
# horizontalalignment='center',
|
||||||
|
# verticalalignment='bottom'
|
||||||
|
# )
|
||||||
|
|
||||||
|
lims = (min(*xs, *ys, *xs_pre, *ys_pre) - self.scale, max(*xs, *ys, *xs_pre, *ys_pre) + self.scale)
|
||||||
|
# ticks = tuple(range(lims[0], lims[1] + self.scale, 2 * self.scale))
|
||||||
|
ticks = np.arange(lims[0], lims[1]+self.scale, 2*self.scale)
|
||||||
|
ax.set_xticks(ticks)
|
||||||
|
ax.set_yticks(ticks)
|
||||||
|
lims = (min(circ.center[0] - circ.radius, circ.center[1] - circ.radius) - self.scale, max(circ.center[0] + circ.radius, circ.center[1] + circ.radius) + self.scale)
|
||||||
|
ax.vlines(0, *lims, color="k", linewidth=1)
|
||||||
|
ax.hlines(0, *lims, color="k", linewidth=1)
|
||||||
|
ax.set_xlim(*lims)
|
||||||
|
ax.set_ylim(*lims)
|
||||||
|
ax.grid("major")
|
||||||
|
ax.set_aspect("equal")
|
||||||
|
ax.set_title(f"G_pt = {self.gray_penalty_transformed:.4f}\nG_po = {self.gray_penalty_original:.4f}")
|
||||||
|
if show:
|
||||||
|
plt.show()
|
||||||
|
|
||||||
|
@staticmethod
|
||||||
|
def most_square_fact(n):
|
||||||
|
best_a = 1
|
||||||
|
best_b = n
|
||||||
|
for a in range(2, n // 2 + 1):
|
||||||
|
b = n // a
|
||||||
|
if b * a != n:
|
||||||
|
continue
|
||||||
|
if (a + b) < (best_a + best_b):
|
||||||
|
best_a = a
|
||||||
|
best_b = b
|
||||||
|
elif (a + b) == (best_a + best_b) and (b - a) < (best_b - best_a):
|
||||||
|
best_a = a
|
||||||
|
best_b = b
|
||||||
|
return (best_a, best_b)
|
||||||
|
|
||||||
|
@staticmethod
|
||||||
|
def calc_distance(a: constellation_point, b: constellation_point, transformed=True):
|
||||||
|
if transformed:
|
||||||
|
return np.linalg.norm(a.transformed - b.transformed)
|
||||||
|
return np.linalg.norm(a.constellation - b.constellation)
|
||||||
|
|
||||||
|
|
||||||
|
class Circle:
|
||||||
|
def __init__(self, center: tuple[float] | None = None, radius: float | None = None):
|
||||||
|
self.center = center
|
||||||
|
self.radius = radius
|
||||||
|
self.points = None
|
||||||
|
|
||||||
|
@staticmethod
|
||||||
|
def welzl(P: list, R: list = []):
|
||||||
|
if len(P) == 0 or len(R) == 3:
|
||||||
|
return Circle.from_points(R)
|
||||||
|
p = P.pop(np.random.randint(0, len(P)))
|
||||||
|
D = Circle.welzl(P.copy(), R.copy())
|
||||||
|
if p in D:
|
||||||
|
return D
|
||||||
|
|
||||||
|
R.append(p)
|
||||||
|
return Circle.welzl(P.copy(), R.copy())
|
||||||
|
|
||||||
|
@staticmethod
|
||||||
|
def from_points(R: list | tuple):
|
||||||
|
# if len(R) != 3:
|
||||||
|
# raise ValueError("R must have length 3")
|
||||||
|
if len(R) == 0:
|
||||||
|
return Circle((0, 0), 1)
|
||||||
|
if len(R) == 1:
|
||||||
|
return Circle(R[0], 1)
|
||||||
|
if len(R) == 2:
|
||||||
|
x = (R[0][0] + R[1][0]) / 2
|
||||||
|
y = (R[0][1] + R[1][1]) / 2
|
||||||
|
r = np.sqrt((R[0][0] - R[1][0]) ** 2 + (R[0][1] - R[1][1]) ** 2) / 2
|
||||||
|
return Circle((x, y), r)
|
||||||
|
if len(R) > 3:
|
||||||
|
raise ValueError("more than 3 points given")
|
||||||
|
z1, z2, z3 = (r[0] + 1j * r[1] for r in R)
|
||||||
|
if z1 == z2 or z2 == z3 or z3 == z1:
|
||||||
|
raise ValueError(f"Duplicate points: {z1}, {z2}, {z3}")
|
||||||
|
w = (z3 - z1) / (z2 - z1)
|
||||||
|
if w == w.real:
|
||||||
|
raise ValueError(f"Points are collinear: {z1}, {z2}, {z3}")
|
||||||
|
|
||||||
|
c = (z2 - z1) * (w - abs(w) ** 2) / (2j * w.imag) + z1
|
||||||
|
r = abs(z1 - c)
|
||||||
|
|
||||||
|
circ = Circle((c.real, c.imag), r)
|
||||||
|
circ.points = R
|
||||||
|
return circ
|
||||||
|
|
||||||
|
def __contains__(self, pt):
|
||||||
|
p = pt[0] + 1j * pt[1]
|
||||||
|
c = self.center[0] + 1j * self.center[1]
|
||||||
|
if np.linalg.norm(p - c) > self.radius:
|
||||||
|
return False
|
||||||
|
return True
|
||||||
|
|
||||||
|
def __repr__(self):
|
||||||
|
return f"Circle [c: ({self.center[0]:.2f}, {self.center[1]:.2f}), r:{self.radius:.2f}]"
|
||||||
|
|
||||||
|
def xy_from_phase(self, phi):
|
||||||
|
if isinstance(phi, float) or isinstance(phi, int):
|
||||||
|
x = self.center[0] + np.cos(phi) * self.radius
|
||||||
|
y = self.center[1] + np.sin(phi) * self.radius
|
||||||
|
return (x, y)
|
||||||
|
num = len(phi)
|
||||||
|
# coords = np.zeros(num)
|
||||||
|
center = np.array(self.center)
|
||||||
|
coords = np.stack((np.cos(phi) * self.radius, np.sin(phi) * self.radius)) + np.expand_dims(
|
||||||
|
center, axis=1
|
||||||
|
).repeat(num, axis=1)
|
||||||
|
return coords
|
||||||
|
|
||||||
|
def plot(self, ax=None, show=True, show_points=False):
|
||||||
|
ax = plt.subplot() if ax is None else ax
|
||||||
|
|
||||||
|
phi = np.arange(0, 2 * np.pi + np.pi / 50, np.pi / 50)
|
||||||
|
xs, ys = self.xy_from_phase(phi)
|
||||||
|
ax.plot(xs, ys, "k:")
|
||||||
|
if self.points is not None and show_points:
|
||||||
|
xs = [p[0] for p in self.points]
|
||||||
|
ys = [p[1] for p in self.points]
|
||||||
|
ax.scatter(xs, ys)
|
||||||
|
ax.grid()
|
||||||
|
ax.set_aspect("equal")
|
||||||
|
|
||||||
|
if show:
|
||||||
|
plt.show()
|
||||||
|
|
||||||
|
|
||||||
|
if __name__ == "__main__":
|
||||||
|
# prevs = []
|
||||||
|
# for i in range(32):
|
||||||
|
# gr = constellation_point.encode(i)
|
||||||
|
# if i > 0:
|
||||||
|
# prev_gr = constellation_point.encode(i-1)
|
||||||
|
# num_mark = i ^ (i -1)
|
||||||
|
# gr_mark = gr ^ (prev_gr)
|
||||||
|
# num_mark = f"{num_mark:05b}".replace("0", " ").replace("1", "|")
|
||||||
|
# gr_mark = f"{gr_mark:05b}".replace("0", " ").replace("1", "|")
|
||||||
|
# print(f" {num_mark} {gr_mark}")
|
||||||
|
# print(f"{i:>2d} [{i:05b}] -> {gr:>2d} [{gr:05b}] {"!" if gr in prevs else ""}")
|
||||||
|
# prevs.append(gr)
|
||||||
|
# exit()
|
||||||
|
N = 32
|
||||||
|
if len(sys.argv) >= 2:
|
||||||
|
N = int(sys.argv[1])
|
||||||
|
|
||||||
|
cnst = constellation(N)
|
||||||
|
cnst.plot(annotate=True)
|
||||||
Reference in New Issue
Block a user