a3x/gr-digitalhf
1
0
Fork 0
mirror of https://github.com/hb9fxq/gr-digitalhf synced 2024-12-22 15:10:00 +00:00
Code Issues Releases Wiki Activity

use more of numpy; improved doppler estimation

Browse source
This commit is contained in:
cmayer 2018-11-06 17:34:48 +01:00
parent fba138b288
commit 39a8426822
9 changed files with 302 additions and 177 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
examples/test_188-110A.grc
View file

@ -191,7 +191,7 @@
</param> </param>
<param> <param>
<key>value</key> <key>value</key>
<value>0.001</value> <value>0.005</value>
</param> </param>
<param> <param>
<key>_enabled</key> <key>_enabled</key>

126
examples/test_188-110C.grc
View file

@ -467,6 +467,57 @@
<value>1</value> <value>1</value>
</param> </param>
</block> </block>
<block>
<key>blocks_file_sink</key>
<param>
<key>append</key>
<value>False</value>
</param>
<param>
<key>alias</key>
<value></value>
</param>
<param>
<key>comment</key>
<value></value>
</param>
<param>
<key>affinity</key>
<value></value>
</param>
<param>
<key>_enabled</key>
<value>True</value>
</param>
<param>
<key>file</key>
<value>/Users/chm/Software/gr-digitalhf/examples/soft_dec.bin</value>
</param>
<param>
<key>_coordinate</key>
<value>(693, 677)</value>
</param>
<param>
<key>_rotation</key>
<value>0</value>
</param>
<param>
<key>id</key>
<value>blocks_file_sink_0</value>
</param>
<param>
<key>type</key>
<value>float</value>
</param>
<param>
<key>unbuffered</key>
<value>True</value>
</param>
<param>
<key>vlen</key>
<value>1</value>
</param>
</block>
<block> <block>
<key>blocks_float_to_complex</key> <key>blocks_float_to_complex</key>
<param> <param>
@ -561,6 +612,53 @@
<value>1</value> <value>1</value>
</param> </param>
</block> </block>
<block>
<key>blocks_pdu_to_tagged_stream</key>
<param>
<key>alias</key>
<value></value>
</param>
<param>
<key>comment</key>
<value></value>
</param>
<param>
<key>affinity</key>
<value></value>
</param>
<param>
<key>_enabled</key>
<value>True</value>
</param>
<param>
<key>_coordinate</key>
<value>(458, 581)</value>
</param>
<param>
<key>_rotation</key>
<value>0</value>
</param>
<param>
<key>id</key>
<value>blocks_pdu_to_tagged_stream_0</value>
</param>
<param>
<key>type</key>
<value>float</value>
</param>
<param>
<key>tag</key>
<value>packet_len</value>
</param>
<param>
<key>maxoutbuf</key>
<value>0</value>
</param>
<param>
<key>minoutbuf</key>
<value>0</value>
</param>
</block>
<block> <block>
<key>blocks_tag_debug</key> <key>blocks_tag_debug</key>
<param> <param>
@ -766,7 +864,7 @@
</param> </param>
<param> <param>
<key>alpha</key> <key>alpha</key>
<value>0.0005</value> <value>0.0001</value>
</param> </param>
<param> <param>
<key>mode</key> <key>mode</key>
@ -2022,11 +2120,11 @@
</param> </param>
<param> <param>
<key>size</key> <key>size</key>
<value>1024</value> <value>768</value>
</param> </param>
<param> <param>
<key>srate</key> <key>srate</key>
<value>samp_rate/sps/2</value> <value>samp_rate/sps/3</value>
</param> </param>
<param> <param>
<key>stemplot</key> <key>stemplot</key>
@ -2046,7 +2144,7 @@
</param> </param>
<param> <param>
<key>tr_mode</key> <key>tr_mode</key>
<value>qtgui.TRIG_MODE_AUTO</value> <value>qtgui.TRIG_MODE_TAG</value>
</param> </param>
<param> <param>
<key>tr_slope</key> <key>tr_slope</key>
@ -2054,11 +2152,11 @@
</param> </param>
<param> <param>
<key>tr_tag</key> <key>tr_tag</key>
<value>"soft_dec"</value> <value>"packet_len"</value>
</param> </param>
<param> <param>
<key>type</key> <key>type</key>
<value>msg_float</value> <value>float</value>
</param> </param>
<param> <param>
<key>update_time</key> <key>update_time</key>
@ -2322,6 +2420,18 @@
<source_key>0</source_key> <source_key>0</source_key>
<sink_key>0</sink_key> <sink_key>0</sink_key>
</connection> </connection>
<connection>
<source_block_id>blocks_pdu_to_tagged_stream_0</source_block_id>
<sink_block_id>blocks_file_sink_0</sink_block_id>
<source_key>0</source_key>
<sink_key>0</sink_key>
</connection>
<connection>
<source_block_id>blocks_pdu_to_tagged_stream_0</source_block_id>
<sink_block_id>qtgui_time_sink_x_1</sink_block_id>
<source_key>0</source_key>
<sink_key>0</sink_key>
</connection>
<connection> <connection>
<source_block_id>blocks_throttle_0</source_block_id> <source_block_id>blocks_throttle_0</source_block_id>
<sink_block_id>fir_filter_xxx_0</sink_block_id> <sink_block_id>fir_filter_xxx_0</sink_block_id>
@ -2348,9 +2458,9 @@
</connection> </connection>
<connection> <connection>
<source_block_id>digitalhf_adaptive_dfe_0</source_block_id> <source_block_id>digitalhf_adaptive_dfe_0</source_block_id>
<sink_block_id>qtgui_time_sink_x_1</sink_block_id> <sink_block_id>blocks_pdu_to_tagged_stream_0</sink_block_id>
<source_key>soft_dec</source_key> <source_key>soft_dec</source_key>
<sink_key>in</sink_key> <sink_key>pdus</sink_key>
</connection> </connection>
<connection> <connection>
<source_block_id>fir_filter_xxx_0</source_block_id> <source_block_id>fir_filter_xxx_0</source_block_id>

54
lib/adaptive_dfe_impl.cc
View file

@ -93,6 +93,7 @@ adaptive_dfe_impl::adaptive_dfe_impl(int sps, // samples per symbol
, _nW(nW) , _nW(nW)
, _mu(mu) , _mu(mu)
, _alpha(alpha) , _alpha(alpha)
, _use_symbol_taps(true)
, _py_module_name(python_module_name) , _py_module_name(python_module_name)
, _physicalLayer() , _physicalLayer()
, _taps_samples(nullptr) , _taps_samples(nullptr)
@ -183,6 +184,7 @@ adaptive_dfe_impl::general_work(int noutput_items,
std::fill_n(_hist_symbols, 2*_nW, gr_complex(0)); std::fill_n(_hist_symbols, 2*_nW, gr_complex(0));
std::fill_n(_taps_samples, _nB+_nF+1, gr_complex(0)); std::fill_n(_taps_samples, _nB+_nF+1, gr_complex(0));
std::fill_n(_taps_symbols, _nW, gr_complex(0)); std::fill_n(_taps_symbols, _nW, gr_complex(0));
_use_symbol_taps = true;
_samples.clear(); _samples.clear();
_phase = -phase_est; _phase = -phase_est;
_taps_samples[_nB+1] = 0.01; _taps_samples[_nB+1] = 0.01;
@ -279,8 +281,10 @@ adaptive_dfe_impl::general_work(int noutput_items,
} }
// publish soft decisions // publish soft decisions
if (!_vec_soft_decisions.empty()) { if (!_vec_soft_decisions.empty()) {
std::cout << "soft_dec " << _vec_soft_decisions.size() << "\n";
unsigned int const bits_per_symbol = _constellations[_constellation_index]->bits_per_symbol(); unsigned int const bits_per_symbol = _constellations[_constellation_index]->bits_per_symbol();
_msg_metadata = pmt::dict_add(_msg_metadata, pmt::mp("bits_per_symbol"), pmt::from_long(bits_per_symbol)); _msg_metadata = pmt::dict_add(_msg_metadata, pmt::mp("bits_per_symbol"), pmt::from_long(bits_per_symbol));
_msg_metadata = pmt::dict_add(_msg_metadata, pmt::mp("packet_len"), pmt::mp(_vec_soft_decisions.size()));
message_port_pub(_msg_port_name, message_port_pub(_msg_port_name,
pmt::cons(_msg_metadata, pmt::cons(_msg_metadata,
pmt::init_f32vector(_vec_soft_decisions.size(), _vec_soft_decisions))); pmt::init_f32vector(_vec_soft_decisions.size(), _vec_soft_decisions)));
@ -382,16 +386,21 @@ gr_complex adaptive_dfe_impl::filter() {
_taps_samples, _taps_samples,
_nB+_nF+1); _nB+_nF+1);
gr_complex dot_symbols=0; gr_complex dot_symbols=0;
for (int l=0; l<_nW; ++l) { gr::digital::constellation_sptr constell = _constellations[_constellation_index];
assert(_hist_symbol_index+l < 2*_nW); bool const update_taps = constell->bits_per_symbol() <= 3;
dot_symbols += _hist_symbols[_hist_symbol_index+l]*_taps_symbols[l]; if (constell->bits_per_symbol() > 3)
_use_symbol_taps = false;
if (_use_symbol_taps) {
for (int l=0; l<_nW; ++l) {
assert(_hist_symbol_index+l < 2*_nW);
dot_symbols += _hist_symbols[_hist_symbol_index+l]*_taps_symbols[l];
}
filter_output += dot_symbols;
} }
filter_output += dot_symbols;
gr_complex known_symbol = _symbols[_symbol_counter]; gr_complex known_symbol = _symbols[_symbol_counter];
bool const is_known = std::abs(known_symbol) > 1e-5; bool const is_known = std::abs(known_symbol) > 1e-5;
if (not is_known) { // not known if (not is_known) { // not known
gr_complex const descrambled_filter_output = std::conj(_scramble[_symbol_counter]) * filter_output; gr_complex const descrambled_filter_output = std::conj(_scramble[_symbol_counter]) * filter_output;
gr::digital::constellation_sptr constell = _constellations[_constellation_index];
unsigned int jc = constell->decision_maker(&descrambled_filter_output); unsigned int jc = constell->decision_maker(&descrambled_filter_output);
gr_complex descrambled_symbol = 0; gr_complex descrambled_symbol = 0;
constell->map_to_points(jc, &descrambled_symbol); constell->map_to_points(jc, &descrambled_symbol);
@ -403,29 +412,24 @@ gr_complex adaptive_dfe_impl::filter() {
_npwr[_constellation_index] = (1-alpha)*_npwr[_constellation_index] + alpha*err; _npwr[_constellation_index] = (1-alpha)*_npwr[_constellation_index] + alpha*err;
std::vector<float> const soft_dec = constell->calc_soft_dec(descrambled_filter_output, _npwr[_constellation_index]); std::vector<float> const soft_dec = constell->calc_soft_dec(descrambled_filter_output, _npwr[_constellation_index]);
std::copy(soft_dec.begin(), soft_dec.end(), std::back_inserter<std::vector<float> >(_vec_soft_decisions)); std::copy(soft_dec.begin(), soft_dec.end(), std::back_inserter<std::vector<float> >(_vec_soft_decisions));
// std::cout << "soft_dec " << _npwr[_constellation_index] << " : ";
// for (int k=0; k<soft_dec.size(); ++k) {
// std::cout << soft_dec[k] << " ";
// }
// std::cout << "\n";
} }
known_symbol = _scramble[_symbol_counter] * descrambled_symbol; known_symbol = _scramble[_symbol_counter] * descrambled_symbol;
} }
gr_complex err = filter_output - known_symbol; if (is_known || update_taps) {
for (int j=0; j<_nB+_nF+1; ++j) { gr_complex const err = filter_output - known_symbol;
assert(_hist_sample_index+j < 2*(_nB+_nF+1)); for (int j=0; j<_nB+_nF+1; ++j) {
_taps_samples[j] -= _mu*err*std::conj(_hist_samples[_hist_sample_index+j]); assert(_hist_sample_index+j < 2*(_nB+_nF+1));
} _taps_samples[j] -= _mu*err*std::conj(_hist_samples[_hist_sample_index+j]);
for (int j=0; j<_nW; ++j) { }
assert(_hist_symbol_index+j < 2*_nW); if (_use_symbol_taps) {
_taps_symbols[j] -= _mu*err*std::conj(_hist_symbols[_hist_symbol_index+j]) + _alpha*_taps_symbols[j]; for (int j=0; j<_nW; ++j) {
} assert(_hist_symbol_index+j < 2*_nW);
// if (_sample_counter < 80*5) _taps_symbols[j] -= _mu*err*std::conj(_hist_symbols[_hist_symbol_index+j]) + _alpha*_taps_symbols[j];
// std::cout << "filter: " << _symbol_counter << " " << _sample_counter << " " << filter_output << " " << known_symbol << " " << std::abs(err) << std::endl; }
if (is_known || true) { _hist_symbols[_hist_symbol_index] = _hist_symbols[_hist_symbol_index + _nW] = known_symbol;
_hist_symbols[_hist_symbol_index] = _hist_symbols[_hist_symbol_index + _nW] = known_symbol; if (++_hist_symbol_index == _nW)
if (++_hist_symbol_index == _nW) _hist_symbol_index = 0;
_hist_symbol_index = 0; }
} }
_descrambled_symbols[_symbol_counter] = filter_output*std::conj(_scramble[_symbol_counter]); _descrambled_symbols[_symbol_counter] = filter_output*std::conj(_scramble[_symbol_counter]);
return filter_output*std::conj(_scramble[_symbol_counter++]); return filter_output*std::conj(_scramble[_symbol_counter++]);

2
lib/adaptive_dfe_impl.h
View file

@ -37,6 +37,8 @@ private:
float _mu; float _mu;
float _alpha; float _alpha;
bool _use_symbol_taps;
// module name w.r.t. digitalhf.physical_layer containing a PhysicalLayer class // module name w.r.t. digitalhf.physical_layer containing a PhysicalLayer class
std::string _py_module_name; std::string _py_module_name;
boost::python::object _physicalLayer; // class instance of physical layer description boost::python::object _physicalLayer; // class instance of physical layer description

1
python/physical_layer/CMakeLists.txt
View file

@ -32,6 +32,7 @@ endif()
GR_PYTHON_INSTALL( GR_PYTHON_INSTALL(
FILES FILES
__init__.py __init__.py
common.py
STANAG_4285.py STANAG_4285.py
MIL_STD_188_110A.py MIL_STD_188_110A.py
MIL_STD_188_110C.py MIL_STD_188_110C.py

110
python/physical_layer/MIL_STD_188_110A.py
View file

@ -2,6 +2,7 @@
from __future__ import print_function from __future__ import print_function
import numpy as np import numpy as np
from common import *
## ---- Walsh-4 codes ----------------------------------------------------------- ## ---- Walsh-4 codes -----------------------------------------------------------
WALSH = np.array([[0,0,0,0, 0,0,0,0], WALSH = np.array([[0,0,0,0, 0,0,0,0],
@ -14,12 +15,9 @@ WALSH = np.array([[0,0,0,0, 0,0,0,0],
[0,1,1,0, 1,0,0,1]], [0,1,1,0, 1,0,0,1]],
dtype=np.uint8) dtype=np.uint8)
def walsh_to_num(w):
return sum(w*(1<<np.arange(8)[::-1]))
FROM_WALSH = -np.ones(256, dtype=np.int8) FROM_WALSH = -np.ones(256, dtype=np.int8)
for i in range(8): for i in range(8):
FROM_WALSH[walsh_to_num(WALSH[i][:])] = i FROM_WALSH[np.packbits(WALSH[i][:])[0]] = i
## ---- tri-bit codes ----------------------------------------------------------- ## ---- tri-bit codes -----------------------------------------------------------
TRIBIT = np.zeros((8,32), dtype=np.uint8) TRIBIT = np.zeros((8,32), dtype=np.uint8)
@ -31,9 +29,6 @@ TRIBIT_SCRAMBLE = np.array(
[7,4,3,0,5,1,5,0,2,2,1,1,5,7,4,3,5,0,2,6,2,1,6,2,0,0,5,0,5,2,6,6], [7,4,3,0,5,1,5,0,2,2,1,1,5,7,4,3,5,0,2,6,2,1,6,2,0,0,5,0,5,2,6,6],
dtype=np.uint8) dtype=np.uint8)
def n_psk(n,x):
return np.complex64(np.exp(2j*np.pi*x/n))
## ---- preamble symbols --------------------------------------------------------- ## ---- preamble symbols ---------------------------------------------------------
D1=D2=C1=C2=C3=0 ## not known D1=D2=C1=C2=C3=0 ## not known
PRE_SYMBOLS = n_psk(2, np.concatenate( PRE_SYMBOLS = n_psk(2, np.concatenate(
@ -41,7 +36,7 @@ PRE_SYMBOLS = n_psk(2, np.concatenate(
PRE_SYMBOLS[9*32:14*32] = 0 PRE_SYMBOLS[9*32:14*32] = 0
## ---- preamble scramble symbols ------------------------------------------------ ## ---- preamble scramble symbols ------------------------------------------------
PRE_SCRAMBLE = n_psk(8, np.concatenate([TRIBIT_SCRAMBLE for i in range(15)])) PRE_SCRAMBLE = n_psk(8, np.concatenate([TRIBIT_SCRAMBLE for _ in range(15)]))
## ---- data scrambler ----------------------------------------------------------- ## ---- data scrambler -----------------------------------------------------------
class ScrambleData(object): class ScrambleData(object):
@ -50,13 +45,13 @@ class ScrambleData(object):
self.reset() self.reset()
def reset(self): def reset(self):
self._state = 0xBAD self._state = 0xBAD
self._counter = 0 self._counter = 0
def next(self): def next(self):
if self._counter == 160: if self._counter == 160:
self.reset() self.reset()
for j in range(8): for _ in range(8):
self._advance() self._advance()
self._counter += 1 self._counter += 1
return self._state&7 return self._state&7
@ -68,13 +63,18 @@ class ScrambleData(object):
self._state ^= 0x053 self._state ^= 0x053
return self._state return self._state
## ---- constellatios -----------------------------------------------------------
BPSK=np.array(zip(np.exp(2j*np.pi*np.arange(2)/2), [0,1]), CONST_DTYPE)
QPSK=np.array(zip(np.exp(2j*np.pi*np.arange(4)/4), [0,1,3,2]), CONST_DTYPE)
PSK8=np.array(zip(np.exp(2j*np.pi*np.arange(8)/8), [0,1,3,2,7,6,4,5]), CONST_DTYPE)
## ---- constellation indices --------------------------------------------------- ## ---- constellation indices ---------------------------------------------------
MODE_BPSK=0 MODE_BPSK=0
MODE_QPSK=1 MODE_QPSK=1
MODE_8PSK=2 MODE_8PSK=2
## ---- mode definitions -------------------------------------------------------- ## ---- mode definitions --------------------------------------------------------
MODE = [[{} for x in range(8)] for y in range(8)] MODE = [[{} for _ in range(8)] for _ in range(8)]
MODE[7][6] = {'bit_rate':4800, 'ci':MODE_8PSK, 'interleaver':['N', 1, 1], 'unknown':32,'known':16, 'nsymb': 1, 'coding_rate': -1 } MODE[7][6] = {'bit_rate':4800, 'ci':MODE_8PSK, 'interleaver':['N', 1, 1], 'unknown':32,'known':16, 'nsymb': 1, 'coding_rate': -1 }
MODE[7][7] = {'bit_rate':2400, 'ci':MODE_8PSK, 'interleaver':['N', 1, 1], 'unknown':32,'known':16, 'nsymb': 1, 'coding_rate':1./2} MODE[7][7] = {'bit_rate':2400, 'ci':MODE_8PSK, 'interleaver':['N', 1, 1], 'unknown':32,'known':16, 'nsymb': 1, 'coding_rate':1./2}
@ -104,15 +104,12 @@ class PhysicalLayer(object):
"""intialization""" """intialization"""
self._sps = sps self._sps = sps
self._frame_counter = -1 self._frame_counter = -1
self._constellations = [self.make_psk(2, [0,1]), self._constellations = [BPSK, QPSK, PSK8]
self.make_psk(4, [0,1,3,2]),
self.make_psk(8, [0,1,3,2,7,6,4,5])] ## TODO: check 8PSK gray code
self._preamble = self.get_preamble() self._preamble = self.get_preamble()
self._pre_counter = -1 self._pre_counter = -1
self._d1d2 = [-1,-1] ## D1,D2 self._d1d2 = [-1,-1] ## D1,D2
self._mode = {} self._mode = {}
self._scr_data = ScrambleData() self._scr_data = ScrambleData()
##self._data = self.get_data()
def get_constellations(self): def get_constellations(self):
return self._constellations return self._constellations
@ -133,18 +130,19 @@ class PhysicalLayer(object):
## ----- data frame ------ ## ----- data frame ------
if self._frame_counter == self._num_frames_per_block: if self._frame_counter == self._num_frames_per_block:
self._frame_counter = 0 self._frame_counter = 0
a = np.zeros(self._frame_len, dtype=[('symb', np.complex64), scramble_for_frame = n_psk(8, np.array([self._scr_data.next()
('scramble', np.complex64)]) for _ in range(self._frame_len)]))
a = np.array(zip(scramble_for_frame,
scramble_for_frame),
dtype=[('symb', np.complex64),
('scramble', np.complex64)])
n_unknown = self._mode['unknown'] n_unknown = self._mode['unknown']
a['symb'] = 1;
a['symb'][0:n_unknown] = 0 a['symb'][0:n_unknown] = 0
if self._frame_counter >= self._num_frames_per_block-2: if self._frame_counter >= self._num_frames_per_block-2:
idx_d1d2 = self._frame_counter - self._num_frames_per_block + 2; idx_d1d2 = self._frame_counter - self._num_frames_per_block + 2;
a['symb'][n_unknown :n_unknown+ 8] *= n_psk(2, WALSH[self._d1d2[idx_d1d2]][:]) a['symb'][n_unknown :n_unknown+ 8] *= n_psk(2, WALSH[self._d1d2[idx_d1d2]][:])
a['symb'][n_unknown+8:n_unknown+16] *= n_psk(2, WALSH[self._d1d2[idx_d1d2]][:]) a['symb'][n_unknown+8:n_unknown+16] *= n_psk(2, WALSH[self._d1d2[idx_d1d2]][:])
a['scramble'] = n_psk(8, np.array([self._scr_data.next() for _ in range(self._frame_len)]))
a['symb'] *= a['scramble']
self._frame_counter += 1 self._frame_counter += 1
return [a, self._mode['ci'],False,True] return [a, self._mode['ci'],False,True]
@ -154,8 +152,7 @@ class PhysicalLayer(object):
[1] ... doppler estimate (rad/symbol) if available""" [1] ... doppler estimate (rad/symbol) if available"""
print('-------------------- get_doppler --------------------', print('-------------------- get_doppler --------------------',
self._frame_counter,len(symbols),len(iq_samples)) self._frame_counter,len(symbols),len(iq_samples))
success = False success,doppler = False,0
doppler = 0
if self._frame_counter == -1: ## -- preamble ---- if self._frame_counter == -1: ## -- preamble ----
success,doppler = self.get_doppler_from_preamble(symbols, iq_samples) success,doppler = self.get_doppler_from_preamble(symbols, iq_samples)
if len(symbols) != 0: if len(symbols) != 0:
@ -174,35 +171,37 @@ class PhysicalLayer(object):
def get_doppler_from_preamble(self, symbols, iq_samples): def get_doppler_from_preamble(self, symbols, iq_samples):
"""quality check and doppler estimation for preamble""" """quality check and doppler estimation for preamble"""
success = True success,doppler = True,0
doppler = 0
if len(iq_samples) != 0: if len(iq_samples) != 0:
zp = np.conj(self.get_preamble_z(self._sps))[9*self._sps:] sps = self._sps
cc = np.array([np.sum(iq_samples[i*self._sps:(3*32+i-9)*self._sps]*zp) zp = np.array([z for z in PhysicalLayer.get_preamble()['symb']
for i in range(4*32)]) for _ in range(sps)], dtype=np.complex64)
acc = np.abs(cc) ## find starting point
for i in range(0,len(cc),32): cc = np.correlate(iq_samples, zp[0:3*32*sps])
print('i=%3d: '%i,end='') imax = np.argmax(np.abs(cc[0:2*32*sps]))
for j in range(32): pks = cc[(imax, imax+3*32*sps),]
print('%3.0f ' % acc[i+j], end='') tpks = cc[imax+3*16*sps:imax+5*16*sps]
print() print('imax=', imax, 'apks=',np.abs(pks),
imax = np.argmax(np.abs(cc[0:2*32])) np.mean(np.abs(pks)), np.mean(np.abs(tpks)), np.abs(tpks))
pks = cc[(imax,imax+3*16,imax+3*16+1,imax+3*32),] success = np.mean(np.abs(pks)) > 2*np.mean(np.abs(tpks))
apks = np.abs(pks) doppler = np.diff(np.unwrap(np.angle(pks)))[0]/(3*32) if success else 0
print('imax=', imax, 'apks=',apks) if success:
success = np.mean(apks[(0,3),]) > 2*np.mean(apks[(1,2),]) idx = np.arange(32*sps)
doppler = np.diff(np.unwrap(np.angle(pks[(0,3),])))[0]/(3*32) if success else 0 pks = [np.correlate(iq_samples[imax+i*32*sps+idx],
zp[i*32*sps+idx])[0]
for i in range(9)]
doppler = freq_est(pks)/32
print('success=', success, 'doppler=', doppler) print('success=', success, 'doppler=', doppler)
return success,doppler return success,doppler
def decode_preamble(self, symbols): def decode_preamble(self, symbols):
data = [FROM_WALSH[walsh_to_num data = [FROM_WALSH[np.packbits
(np.real (np.real
(np.sum (np.sum
(symbols[i:i+32].reshape((4,8)),0))<0)] (symbols[i:i+32].reshape((4,8)),0))<0)[0]]
for i in range(0,15*32,32)] for i in range(0,15*32,32)]
print('data=',data) print('data=',data)
self._pre_counter = sum((np.array(data[11:14])&3)*(1<<2*np.arange(3)[::-1])) self._pre_counter = sum([(x&3)*(1<<2*y) for (x,y) in zip(data[11:14][::-1], range(3))])
self._d1d2 = data[9:11] self._d1d2 = data[9:11]
self._mode = MODE[data[9]][data[10]] self._mode = MODE[data[9]][data[10]]
self._block_len = 11520 if self._mode['interleaver'][0] == 'L' else 1440 self._block_len = 11520 if self._mode['interleaver'][0] == 'L' else 1440
@ -213,27 +212,17 @@ class PhysicalLayer(object):
@staticmethod @staticmethod
def get_preamble(): def get_preamble():
"""preamble symbols + scrambler""" """preamble symbols + scrambler"""
a=np.zeros(15*32, dtype=[('symb', np.complex64), return np.array(zip(PRE_SCRAMBLE*PRE_SYMBOLS,
('scramble', np.complex64)]) PRE_SCRAMBLE),
a['symb'] = PRE_SCRAMBLE*PRE_SYMBOLS dtype=[('symb', np.complex64),
a['scramble'] = PRE_SCRAMBLE ('scramble', np.complex64)])
return a
@staticmethod @staticmethod
def get_preamble_z(sps): def get_preamble_z(sps):
"""preamble symbols for preamble correlation""" """preamble symbols for preamble correlation"""
a = PhysicalLayer.get_preamble() a = PhysicalLayer.get_preamble()
return np.array([z for z in a['symb'][0:32*3] return np.array([z for z in a['symb'][0:32*3]
for i in range(sps)]) for _ in range(sps)])
@staticmethod
def make_psk(n, gray_code):
"""generates n-PSK constellation data"""
c = np.zeros(n, dtype=[('points', np.complex64),
('symbols', np.uint8)])
c['points'] = n_psk(n,np.arange(n))
c['symbols'] = gray_code
return c
if __name__ == '__main__': if __name__ == '__main__':
def gen_data_scramble(): def gen_data_scramble():
@ -245,13 +234,13 @@ if __name__ == '__main__':
a = np.zeros(160, dtype=np.uint8) a = np.zeros(160, dtype=np.uint8)
s = 0xBAD s = 0xBAD
for i in range(160): for i in range(160):
for j in range(8): s = advance(s) for _ in range(8): s = advance(s)
a[i] = s&7; a[i] = s&7;
return a return a
p=PhysicalLayer(5) p=PhysicalLayer(5)
z1=np.array([x for x in PRE_SYMBOLS for i in range(5)]) z1=np.array([x for x in PRE_SYMBOLS for _ in range(5)])
z2=np.array([x for x in PRE_SCRAMBLE for i in range(5)]) z2=np.array([x for x in PRE_SCRAMBLE for _ in range(5)])
z=z1*z2 z=z1*z2
for i in range(3): for i in range(3):
@ -264,6 +253,5 @@ if __name__ == '__main__':
print(gen_data_scramble()) print(gen_data_scramble())
s=ScrambleData() s=ScrambleData()
print(type(s))
print([s.next() for _ in range(160)]) print([s.next() for _ in range(160)])
print([s.next() for _ in range(160)]) print([s.next() for _ in range(160)])

117
python/physical_layer/MIL_STD_188_110C.py
View file

@ -2,13 +2,9 @@
from __future__ import print_function from __future__ import print_function
import numpy as np import numpy as np
from common import *
def n_psk(n,x):
return np.complex64(np.exp(2j*np.pi*x/n))
## ---- constellations ----------------------------------------------------------- ## ---- constellations -----------------------------------------------------------
CONST_DTYPE=np.dtype([('points', np.complex64),
('symbols', np.uint8)])
BPSK=np.array(zip(np.exp(2j*np.pi*np.arange(2)/2), [0,1]), CONST_DTYPE) BPSK=np.array(zip(np.exp(2j*np.pi*np.arange(2)/2), [0,1]), CONST_DTYPE)
QPSK=np.array(zip(np.exp(2j*np.pi*np.arange(4)/4), [0,1,3,2]), CONST_DTYPE) QPSK=np.array(zip(np.exp(2j*np.pi*np.arange(4)/4), [0,1,3,2]), CONST_DTYPE)
PSK8=np.array(zip(np.exp(2j*np.pi*np.arange(8)/8), [0,1,3,2,7,6,4,5]), CONST_DTYPE) PSK8=np.array(zip(np.exp(2j*np.pi*np.arange(8)/8), [0,1,3,2,7,6,4,5]), CONST_DTYPE)
@ -47,6 +43,10 @@ QAM64=np.array(
-0.360142-0.932897j, -0.353057-0.588429j, -0.353057-0.117686j, -0.353057-0.353057j], -0.360142-0.932897j, -0.353057-0.588429j, -0.353057-0.117686j, -0.353057-0.353057j],
range(64)), CONST_DTYPE) range(64)), CONST_DTYPE)
## for test
QAM64 = QAM64[(7,3,24,56,35,39,60,28),]
QAM64['symbols'] = [1, 0, 2, 6, 4, 5, 7, 3]
## ---- constellation indices --------------------------------------------------- ## ---- constellation indices ---------------------------------------------------
MODE_BPSK = 0 MODE_BPSK = 0
MODE_QPSK = 1 MODE_QPSK = 1
@ -62,20 +62,18 @@ class ScrambleData(object):
self.reset() self.reset()
def reset(self): def reset(self):
self._state = 1 self._state = np.array([0,0,0,0,0,0,0,0,1], dtype=np.bool)
self._taps = np.array([0,0,0,0,1,0,0,0,1], dtype=np.bool)
def next(self, num_bits): def next(self, num_bits):
r = self._state & ((1<<num_bits)-1) r = np.packbits(self._state[1:])[0]&((1<<num_bits)-1)
for i in range(num_bits): for _ in range(num_bits):
self._advance() self._advance()
return r return r
def _advance(self): def _advance(self):
lsb = self._state&1 self._state = np.concatenate(([np.sum(self._state&self._taps)&1],
self._state = (self._state>>1)&511 self._state[0:-1]))
if lsb:
self._state ^= 0x10B
return self._state
## ---- preamble definitions --------------------------------------------------- ## ---- preamble definitions ---------------------------------------------------
## 184 = 8*23 ## 184 = 8*23
@ -132,10 +130,13 @@ class PhysicalLayer(object):
## --- preamble frame ---- ## --- preamble frame ----
if self._frame_counter == -1: if self._frame_counter == -1:
return [self._preamble,MODE_BPSK,True,False] return [self._preamble,MODE_BPSK,True,False]
## ----- data frame ------ ## ----- (re)inserted preamble ------
if self._frame_counter == 0: if self._frame_counter == 0:
self.a = self.make_reinserted_preamble() self.a = self.make_reinserted_preamble()
return [self.a, MODE_QPSK,False,True] return [self.a, MODE_QPSK,False,False]
if self._frame_counter >= 1:
self.a = self.make_data_frame()
return [self.a, MODE_64QAM,False,True]
def get_doppler(self, symbols, iq_samples): def get_doppler(self, symbols, iq_samples):
"""returns a tuple """returns a tuple
@ -145,50 +146,45 @@ class PhysicalLayer(object):
self._frame_counter,len(symbols),len(iq_samples)) self._frame_counter,len(symbols),len(iq_samples))
#if len(symbols)!=0: #if len(symbols)!=0:
# print('symb=', symbols) # print('symb=', symbols)
success = False success,doppler = False,0
doppler = 0
if self._frame_counter == -1: ## -- preamble ---- if self._frame_counter == -1: ## -- preamble ----
success,doppler = self.get_doppler_from_preamble(symbols, iq_samples) success,doppler = self.get_doppler_from_preamble(symbols, iq_samples)
if len(symbols) != 0: if len(symbols) != 0:
for s in symbols: for s in symbols:
print(s) print(s)
self._frame_counter = 0 self._frame_counter = 0
elif self._frame_counter >= 0 and self._frame_counter <5: ## -- reinserted preamble ----
for s in symbols:
print(s)
self._frame_counter += 1
success = True
else: ## ------------------------ data frame ---- else: ## ------------------------ data frame ----
if len(symbols) != 0: for s in symbols:
for s in symbols: print(s)
print(s)
success = False success = False
self._frame_counter = -1 self._frame_counter = -1
return success,doppler return success,doppler
def get_doppler_from_preamble(self, symbols, iq_samples): def get_doppler_from_preamble(self, symbols, iq_samples):
"""quality check and doppler estimation for preamble""" """quality check and doppler estimation for preamble"""
success = True success,doppler = True,0
doppler = 0
shift=9
if len(iq_samples) != 0: if len(iq_samples) != 0:
zp = np.conj(self.get_preamble_z(self._sps)[shift*self._sps:]) sps = self._sps
cc = np.array([np.sum(iq_samples[i:i+23*self._sps] * idx = np.arange(23*sps)
zp[0:23*self._sps]) zp = self.get_preamble_z(self._sps)
for i in range(23*3*self._sps)]) cc = np.correlate(iq_samples, zp[idx])
acc = np.abs(cc) imax = np.argmax(np.abs(cc[0:23*sps]))
for i in range(0,len(cc),23*self._sps): pks = [np.correlate(iq_samples[imax+i*23*sps+idx],
print('i=%3d: '%i,end='') zp[i*23*sps+idx])[0]
for j in range(23*self._sps): for i in range(7)]
print('%3.0f ' % acc[i+j], end='') success = np.mean(np.abs(pks)) > 2*np.mean(np.abs(cc[imax+11*sps+range(-sps,sps)]))
print() print('test:',imax, np.mean(np.abs(pks)), np.mean(np.abs(cc[imax+11*sps+range(-sps,sps)])))
if success:
imax = np.argmax(np.abs(cc[0:3*23*self._sps])) print('doppler apks', np.abs(pks))
print(imax) print('doppler ppks', np.angle(pks),
pks = np.array([np.sum(iq_samples[(imax+23*i*self._sps): np.diff(np.unwrap(np.angle(pks)))/23,
(imax+23*i*self._sps+23*self._sps)] * np.mean(np.diff(np.unwrap(np.angle(pks)))/23))
zp[(23*i*self._sps): doppler = freq_est(pks[1:])/23;
(23*i*self._sps+23*self._sps)])
for i in range(1,5)])
print('doppler apks', np.abs(pks))
print('doppler ppks', np.angle(pks), np.diff(np.unwrap(np.angle(pks)))/23)
doppler = np.mean(np.diff(np.unwrap(np.angle(pks))))/23
success = True
print('success=', success, 'doppler=', doppler) print('success=', success, 'doppler=', doppler)
return success,doppler return success,doppler
@ -199,6 +195,13 @@ class PhysicalLayer(object):
('scramble', np.complex64)]) ('scramble', np.complex64)])
a['symb'][32:32+3*13] = 0 ## D0,D1,D2 a['symb'][32:32+3*13] = 0 ## D0,D1,D2
return a return a
def make_data_frame(self):
a=np.zeros(256+31, dtype=[('symb', np.complex64),
('scramble', np.complex64)])
a['scramble'] = 1
a['symb'][256:] = MINI_PROBE_MINUS
a['scramble'][256:] = MINI_PROBE_MINUS
return a
@staticmethod @staticmethod
def get_preamble(): def get_preamble():
@ -210,7 +213,7 @@ class PhysicalLayer(object):
@staticmethod @staticmethod
def get_preamble_z(sps): def get_preamble_z(sps):
"""preamble symbols for preamble correlation""" """preamble symbols for preamble correlation"""
return np.array([z for z in PREAMBLE for i in range(sps)]) return np.array([z for z in PREAMBLE for _ in range(sps)])
if __name__ == '__main__': if __name__ == '__main__':
print(PREAMBLE) print(PREAMBLE)
@ -224,12 +227,16 @@ if __name__ == '__main__':
print([s.next(1) for _ in range(511)]) print([s.next(1) for _ in range(511)])
print([s.next(1) for _ in range(511)] == print([s.next(1) for _ in range(511)] ==
[s.next(1) for _ in range(511)]) [s.next(1) for _ in range(511)])
print(QAM64) #print(QAM64)
print(QAM32) #print(QAM32)
print(QAM16) #print(QAM16)
print(PSK8) #print(PSK8)
print(QPSK) #print(QPSK)
print(BPSK) #print(BPSK)
print(MINI_PROBE_PLUS) #print(MINI_PROBE_PLUS)
print(MINI_PROBE_MINUS) #print(MINI_PROBE_MINUS)
print(MINI_PROBE_PLUS*MINI_PROBE_MINUS) #print(MINI_PROBE_PLUS*MINI_PROBE_MINUS)
#for i in range(len(QAM64)):
# print(QAM64['points'][i])
print([s.next(6) for _ in range(256)])

38
python/physical_layer/STANAG_4285.py
View file

@ -57,13 +57,15 @@ class PhysicalLayer(object):
success = True success = True
doppler = 0 doppler = 0
if len(iq_samples) != 0: if len(iq_samples) != 0:
zp = [x for x in self._preamble['symb'][9:40] for i in range(self._sps)] sps = self._sps
cc = np.array([np.sum(iq_samples[ i*5:(31+i)*5]*zp) for i in range(49)]) zp = np.array([x for x in self._preamble['symb'][9:40]
imax = np.argmax(np.abs(cc[0:18])) for _ in range(sps)], dtype=np.complex64)
pks = cc[(imax,imax+15,imax+16,imax+31),] cc = np.correlate(iq_samples, zp)
apks = np.abs(pks) imax = np.argmax(np.abs(cc[0:18*sps]))
success = np.mean(apks[(0,3),]) > 2*np.mean(apks[(1,2),]) pks = cc[(imax,imax+31*sps),]
doppler = np.diff(np.unwrap(np.angle(pks[(0,3),])))[0]/31 if success else 0 tpks = cc[imax+15*sps:imax+16*sps]
success = np.mean(np.abs(pks)) > 2*np.mean(np.abs(tpks))
doppler = np.diff(np.unwrap(np.angle(pks)))[0]/31 if success else 0
if len(symbols) != 0: if len(symbols) != 0:
idx = range(30,80) if self._is_first_frame else range(80) idx = range(30,80) if self._is_first_frame else range(80)
z = symbols[idx]*np.conj(self._preamble['symb'][idx]) z = symbols[idx]*np.conj(self._preamble['symb'][idx])
@ -102,7 +104,7 @@ class PhysicalLayer(object):
p = np.zeros(176, dtype=np.uint8) p = np.zeros(176, dtype=np.uint8)
for i in range(176): for i in range(176):
p[i] = np.sum(state[-3:]*[4,2,1]) p[i] = np.sum(state[-3:]*[4,2,1])
for j in range(3): for _ in range(3):
state = np.concatenate(([np.sum(state&taps)&1], state[0:-1])) state = np.concatenate(([np.sum(state&taps)&1], state[0:-1]))
a=np.zeros(176, dtype=[('symb',np.complex64), ('scramble', np.complex64)]) a=np.zeros(176, dtype=[('symb',np.complex64), ('scramble', np.complex64)])
## 8PSK modulation ## 8PSK modulation
@ -116,24 +118,6 @@ class PhysicalLayer(object):
def make_psk(n, gray_code): def make_psk(n, gray_code):
"""generates n-PSK constellation data""" """generates n-PSK constellation data"""
c = np.zeros(n, dtype=[('points', np.complex64), ('symbols', np.uint8)]) c = np.zeros(n, dtype=[('points', np.complex64), ('symbols', np.uint8)])
c['points'] = np.exp(2*np.pi*1j*np.array(range(n))/n) c['points'] = np.exp(2*np.pi*1j*np.arange(n)/n)
c['symbols'] = gray_code c['symbols'] = gray_code
return c return c
## for now not used (doppler estimation after adaptive filtering does not work)
@staticmethod
def data_aided_frequency_estimation(x,c):
"""Data-Aided Frequency Estimation for Burst Digital Transmission,
Umberto Mengali and M. Morelli, IEEE TRANSACTIONS ON COMMUNICATIONS,
VOL. 45, NO. 1, JANUARY 1997"""
z = x*np.conj(c) ## eq (2)
L0 = len(z)
N = L0//2
R = np.zeros(N, dtype=np.complex64)
for i in range(N):
R[i] = 1.0/(L0-i)*np.sum(z[i:]*np.conj(z[0:L0-i])) ## eq (3)
m = np.array(range(N), dtype=np.float)
w = 3*((L0-m)*(L0-m+1)-N*(L0-N))/(N*(4*N*N - 6*N*L0 + 3*L0*L0-1)) ## eq (9)
mod_2pi = lambda x : np.mod(x-np.pi, 2*np.pi) - np.pi
fd = np.sum(w[1:] * mod_2pi(np.diff(np.angle(R)))) ## eq (8)
return fd

29
python/physical_layer/common.py Normal file
View file

@ -0,0 +1,29 @@
## -*- python -*-
import numpy as np
CONST_DTYPE=np.dtype([('points', np.complex64),
('symbols', np.uint8)])
def n_psk(n,x):
"""n-ary PSK constellation"""
return np.complex64(np.exp(2j*np.pi*x/n))
def freq_est(z):
"""Data-Aided Frequency Estimation for Burst Digital Transmission,
Umberto Mengali and M. Morelli, IEEE TRANSACTIONS ON COMMUNICATIONS,
VOL. 45, NO. 1, JANUARY 1997"""
L0 = len(z)
N = L0//2
R = np.zeros(N+1, dtype=np.complex64)
for i in range(N+1):
R[i] = 1.0/(L0-i)*np.sum(z[i:]*np.conj(z[0:L0-i])) ## eq (3)
m = np.arange(N+1, dtype=np.float32)
w = 3*((L0-m)*(L0-m+1)-N*(L0-N))/(N*(4*N*N - 6*N*L0 + 3*L0*L0-1)) ## eq (9)
mod_2pi = lambda x : np.mod(x-np.pi, 2*np.pi) - np.pi
return np.sum(w[1:] * mod_2pi(np.diff(np.angle(R)))) ## eq (8)
if __name__ == '__main__':
idx=np.arange(3)
z=np.exp(1j*idx*0.056+1j)
print(freq_est(z)/0.056)