-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathber_util.py
182 lines (159 loc) · 5.89 KB
/
ber_util.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
import numpy as np
import matplotlib.pyplot as plt
# from commpy.modulation import PSKModem, QAMModem # noqa
from commpy.filters import rrcosfilter # scikit-commpy
# import pandas as pd
dict_BpS = dict({
'bpsk': 1, # noqa
'qpsk': 2, # noqa
'qam4': 2,
'qam16': 4,
'qam32': 5,
'qam64': 6,
})
# def calc_ber():
# raise NotImplementedError
def get_h(fs=10, g_delay=4, debug=False):
# raised cosine (RC) filter (FIR) impulse response
# TODO // verify and validate
fd = 1
alpha = 0.3
n = 2*g_delay*fs + 1
h = rrcosfilter(N=n, alpha=alpha, Ts=fd, Fs=fs)[1]
h_norm = np.divide(h, np.sqrt(np.square(h).sum()), dtype=np.float16)
if debug:
plt.figure()
plt.plot(h_norm)
plt.show()
# https://stackoverflow.com/a/25858023
return h_norm
def gen_data(n, mod, seed):
"""
:param n: number of symbol to generate
:param mod: modulation type, valid options : ['bpsk','qpsk']
:param seed: seed value for random data generation
:return:
syms : (n, ) 1D numpy array, modulated symbols; no differentiation for real(I) and imaginary(Q) parts
bits : (nb, ) 1D numpy array, 0/1 message bits
"""
np.random.seed(seed)
# n : 'number of symbols'
# nb : 'number of bits'
nb = dict_BpS[mod] * n
bits = np.random.randint(low=0, high=2, size=nb, dtype=np.int8)
if mod == 'bpsk': # noqa
syms = 2 * bits - 1
# syms = np.reshape(syms, (-1, 1))
# make 1D real valued BPSK symbols to 2D (real, 0*imag) format
# syms = np.concatenate((syms, np.zeros((n, 1))), axis=1)
elif mod == 'qpsk': # noqa
# 1/np.sqrt(2) = 0.7071067811865475
# s = (2 * bits - 1) /np.sqrt(2)
# s = (2 * bits - 1) * 0.7071067811865475
# syms = bits.reshape(-1, 2).astype(np.float16)
syms = bits.reshape(-1).astype(np.float16) # flatten the symbols I and Q data into 1D array
# s[s==[0, 0]] = [-0.7071 -0.7071i]
# s[s==[0, 1]] = [-0.7071 +0.7071i]
# s[s==[1, 0]] = [0.7071 -0.7071i]
# s[s==[1, 1]] = [0.7071 +0.7071i]
syms[syms == 0] = -0.7071
syms[syms == 1] = 0.7071
# INFO data format for QPSK
# s.shape() = (N, 2) >> real(I) (x-axis) 0th col, imag(Q) (y-axis) 1th col
# [I1, Q1]
# [I2, Q2]
# ... ...
# [IN, QN]
else:
raise NotImplementedError
return syms, bits
# reshape notes:
# a = np.array([1, 2, 3, 4, 5, 6, 7, 8])
# b = a.reshape(-1, 2)
# [[1, 2],
# [3, 4],
# [5, 6],
# [7, 8]]
# c = b.flatten()
# array([1, 2, 3, 4, 5, 6, 7, 8])
# def awgn(signal, desired_snr):
# """
# Add AWGN to the input signal to achieve the desired SNR level.
#
# source: https://saturncloud.io/blog/python-numpy-implementing-an-additive-white-gaussian-noise-function/
# """
# # Calculate the power of the signal
# signal_power = np.mean(signal ** 2)
#
# # Calculate the noise power based on the desired SNR and signal power
# noise_power = signal_power / (10 ** (desired_snr / 10))
#
# # Generate the noise with the calculated power
# noise = np.random.normal(0, np.sqrt(noise_power), len(signal))
#
# # Add the noise to the original signal
# noisy_signal = signal + noise
#
# return noisy_signal
# def add_awgn2(x_volts, target_snr_db):
# # https://www.rfwireless-world.com/source-code/Python/AWGN-python-script.html
#
# x_watts = x_volts ** 2
# x_db = 10 * np.log10(x_watts)
#
# # Calculate signal power and convert to dB
# sig_avg_watts = np.mean(x_watts)
# sig_avg_db = 10 * np.log10(sig_avg_watts)
# # print("SMR in dB = ", sig_avg_db)
#
# # Calculate noise and convert it to watts
# noise_avg_db = sig_avg_db - target_snr_db
# # print("Average Noise power in dB = ", noise_avg_db)
# noise_avg_watts = 10 ** (noise_avg_db / 10)
# # Generate samples of white noise
# mean_noise = 0
# noise_volts = np.random.normal(mean_noise, np.sqrt(noise_avg_watts), len(x_watts))
#
# # Add noise to original sine waveform signal
# y_volts = x_volts + noise_volts
#
# return y_volts
def add_awgn(inputs, snr=10, seed=1234): # noqa
if snr == 'NoNoise':
return inputs
assert len(inputs.shape) == 1 or len(inputs.shape) == 2, 'Only 1D and 2D data are supported!'
np.random.seed(seed)
n = len(inputs)
outputs = np.empty(inputs.shape).astype(np.float16) # TODO FIX ERROR RuntimeWarning: overflow encountered in cast
# SNR = 10*log10(Eb/No)
# Eb/No = 10 ^(SNR/10)
# EB = 1 for BPSK # noqa
n0 = 10 ** (-snr / 10)
# noise = np.sqrt(N0/2)*(np.randn(n, 1) + np.randn(n, 1)) # noqa
# noise_r = np.multiply(np.sqrt(n0 / 2), np.random.standard_normal(n)) # standard_normal: (mean=0, stdev=1) # noqa
noise_r = np.multiply(np.sqrt(n0 / 2), np.random.randn(n).astype(np.float16)) # standard_normal: (mean=0, stdev=1) # noqa
if len(inputs.shape) == 1:
outputs = np.add(inputs, noise_r)
else: # 2d case
outputs[:, 0] = np.add(inputs[:, 0], noise_r)
if inputs.shape[1] == 2:
# noise_i = np.multiply(np.sqrt(n0 / 2), np.random.standard_normal(n)) # standard_normal: (mean=0, stdev=1) # noqa
noise_i = np.multiply(np.sqrt(n0 / 2), np.random.randn(n)) # standard_normal: (mean=0, stdev=1) # noqa
outputs[:, 1] = np.add(inputs[:, 1], noise_i)
# TODO verify noise addition in complex case
# plt.hist(noise, bins=40)
# plt.show()
return outputs
def bit_checker(bit_ref, bit_tc):
"""
:param bit_ref: reference bit
:param bit_tc: bits to check
:return: noe and nob
# noe : number of error
# nob : number of bit
"""
assert len(bit_ref) == len(bit_tc), 'mismatch on the input lengths'
# noe = abs(np.subtract(bit_ref, bit_tc)).sum()
nob = len(bit_ref)
noe = nob - np.equal(bit_ref, bit_tc).sum()
return noe, nob