Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 1 | /* -*- c++ -*- */ |
| 2 | /* |
| 3 | * @file |
| 4 | * @author (C) 2009-2017 by Piotr Krysik <ptrkrysik@gmail.com> |
| 5 | * @author Contributions by sysmocom - s.f.m.c. GmbH / Eric Wild <ewild@sysmocom.de> |
| 6 | * @section LICENSE |
| 7 | * |
| 8 | * Gr-gsm is free software; you can redistribute it and/or modify |
| 9 | * it under the terms of the GNU General Public License as published by |
| 10 | * the Free Software Foundation; either version 3, or (at your option) |
| 11 | * any later version. |
| 12 | * |
| 13 | * Gr-gsm is distributed in the hope that it will be useful, |
| 14 | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 15 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| 16 | * GNU General Public License for more details. |
| 17 | * |
| 18 | * You should have received a copy of the GNU General Public License |
| 19 | * along with gr-gsm; see the file COPYING. If not, write to |
| 20 | * the Free Software Foundation, Inc., 51 Franklin Street, |
| 21 | * Boston, MA 02110-1301, USA. |
| 22 | */ |
| 23 | |
| 24 | #include "constants.h" |
| 25 | |
| 26 | #ifdef HAVE_CONFIG_H |
| 27 | #include "config.h" |
| 28 | #endif |
| 29 | #include <complex> |
| 30 | |
| 31 | |
| 32 | #include <algorithm> |
| 33 | #include <string.h> |
| 34 | #include <iostream> |
| 35 | #include <numeric> |
| 36 | #include <vector> |
| 37 | #include <fstream> |
| 38 | |
| 39 | #include "viterbi_detector.h" |
| 40 | #include "grgsm_vitac.h" |
| 41 | |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 42 | gr_complex d_acc_training_seq[N_ACCESS_BITS]; ///<encoded training sequence of a RACH burst |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 43 | gr_complex d_sch_training_seq[N_SYNC_BITS]; ///<encoded training sequence of a SCH burst |
| 44 | gr_complex d_norm_training_seq[TRAIN_SEQ_NUM][N_TRAIN_BITS]; ///<encoded training sequences of a normal and dummy burst |
| 45 | const int d_chan_imp_length = CHAN_IMP_RESP_LENGTH; |
| 46 | |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 47 | void initvita() |
| 48 | { |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 49 | /** |
| 50 | * Prepare SCH sequence bits |
| 51 | * |
| 52 | * (TS_BITS + 2 * GUARD_PERIOD) |
| 53 | * Burst and two guard periods |
| 54 | * (one guard period is an arbitrary overlap) |
| 55 | */ |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 56 | gmsk_mapper(SYNC_BITS, N_SYNC_BITS, d_sch_training_seq, gr_complex(0.0, -1.0)); |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 57 | for (auto &i : d_sch_training_seq) |
| 58 | i = conj(i); |
| 59 | |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 60 | /* ab */ |
| 61 | gmsk_mapper(ACCESS_BITS, N_ACCESS_BITS, d_acc_training_seq, gr_complex(0.0, -1.0)); |
| 62 | for (auto &i : d_acc_training_seq) |
| 63 | i = conj(i); |
| 64 | |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 65 | /* Prepare bits of training sequences */ |
| 66 | for (int i = 0; i < TRAIN_SEQ_NUM; i++) { |
| 67 | /** |
| 68 | * If first bit of the sequence is 0 |
| 69 | * => first symbol is 1, else -1 |
| 70 | */ |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 71 | gr_complex startpoint = train_seq[i][0] == 0 ? gr_complex(1.0, 0.0) : gr_complex(-1.0, 0.0); |
| 72 | gmsk_mapper(train_seq[i], N_TRAIN_BITS, d_norm_training_seq[i], startpoint); |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 73 | for (auto &i : d_norm_training_seq[i]) |
| 74 | i = conj(i); |
| 75 | } |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 76 | } |
| 77 | |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 78 | template <unsigned int burst_size> |
| 79 | NO_UBSAN static void detect_burst_generic(const gr_complex *input, gr_complex *chan_imp_resp, int burst_start, |
| 80 | char *output_binary, int ss) |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 81 | { |
| 82 | std::vector<gr_complex> rhh_temp(CHAN_IMP_RESP_LENGTH * d_OSR); |
| 83 | unsigned int stop_states[2] = { 4, 12 }; |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 84 | gr_complex filtered_burst[burst_size]; |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 85 | gr_complex rhh[CHAN_IMP_RESP_LENGTH]; |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 86 | float output[burst_size]; |
| 87 | int start_state = ss; |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 88 | |
| 89 | autocorrelation(chan_imp_resp, &rhh_temp[0], d_chan_imp_length * d_OSR); |
| 90 | for (int ii = 0; ii < d_chan_imp_length; ii++) |
| 91 | rhh[ii] = conj(rhh_temp[ii * d_OSR]); |
| 92 | |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 93 | mafi(&input[burst_start], burst_size, chan_imp_resp, d_chan_imp_length * d_OSR, filtered_burst); |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 94 | |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 95 | viterbi_detector(filtered_burst, burst_size, rhh, start_state, stop_states, 2, output); |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 96 | |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 97 | for (unsigned int i = 0; i < burst_size; i++) |
Eric Wild | 56c7b77 | 2024-02-21 18:58:42 +0100 | [diff] [blame^] | 98 | output_binary[i] = output[i] > 0 ? -127 : 127; // pre flip bits! |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 99 | } |
| 100 | |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 101 | NO_UBSAN void detect_burst_nb(const gr_complex *input, gr_complex *chan_imp_resp, int burst_start, char *output_binary, |
| 102 | int ss) |
| 103 | { |
| 104 | return detect_burst_generic<BURST_SIZE>(input, chan_imp_resp, burst_start, output_binary, ss); |
| 105 | } |
| 106 | NO_UBSAN void detect_burst_ab(const gr_complex *input, gr_complex *chan_imp_resp, int burst_start, char *output_binary, |
| 107 | int ss) |
| 108 | { |
| 109 | return detect_burst_generic<8 + 41 + 36 + 3>(input, chan_imp_resp, burst_start, output_binary, ss); |
| 110 | } |
| 111 | |
| 112 | NO_UBSAN void detect_burst_nb(const gr_complex *input, gr_complex *chan_imp_resp, int burst_start, char *output_binary) |
| 113 | { |
| 114 | return detect_burst_nb(input, chan_imp_resp, burst_start, output_binary, 3); |
| 115 | } |
| 116 | NO_UBSAN void detect_burst_ab(const gr_complex *input, gr_complex *chan_imp_resp, int burst_start, char *output_binary) |
| 117 | { |
| 118 | return detect_burst_ab(input, chan_imp_resp, burst_start, output_binary, 3); |
| 119 | } |
| 120 | |
| 121 | void gmsk_mapper(const unsigned char *input, int nitems, gr_complex *gmsk_output, gr_complex start_point) |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 122 | { |
| 123 | gr_complex j = gr_complex(0.0, 1.0); |
| 124 | gmsk_output[0] = start_point; |
| 125 | |
| 126 | int previous_symbol = 2 * input[0] - 1; |
| 127 | int current_symbol; |
| 128 | int encoded_symbol; |
| 129 | |
| 130 | for (int i = 1; i < nitems; i++) { |
| 131 | /* Change bits representation to NRZ */ |
| 132 | current_symbol = 2 * input[i] - 1; |
| 133 | |
| 134 | /* Differentially encode */ |
| 135 | encoded_symbol = current_symbol * previous_symbol; |
| 136 | |
| 137 | /* And do GMSK mapping */ |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 138 | gmsk_output[i] = j * gr_complex(encoded_symbol, 0.0) * gmsk_output[i - 1]; |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 139 | |
| 140 | previous_symbol = current_symbol; |
| 141 | } |
| 142 | } |
| 143 | |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 144 | gr_complex correlate_sequence(const gr_complex *sequence, int length, const gr_complex *input) |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 145 | { |
| 146 | gr_complex result(0.0, 0.0); |
| 147 | |
| 148 | for (int ii = 0; ii < length; ii++) |
| 149 | result += sequence[ii] * input[ii * d_OSR]; |
| 150 | |
| 151 | return conj(result) / gr_complex(length, 0); |
| 152 | } |
| 153 | |
| 154 | /* Computes autocorrelation for positive arguments */ |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 155 | inline void autocorrelation(const gr_complex *input, gr_complex *out, int nitems) |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 156 | { |
| 157 | for (int k = nitems - 1; k >= 0; k--) { |
| 158 | out[k] = gr_complex(0, 0); |
| 159 | for (int i = k; i < nitems; i++) |
| 160 | out[k] += input[i] * conj(input[i - k]); |
| 161 | } |
| 162 | } |
| 163 | |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 164 | inline void mafi(const gr_complex *input, int nitems, gr_complex *filter, int filter_length, gr_complex *output) |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 165 | { |
| 166 | for (int n = 0; n < nitems; n++) { |
| 167 | int a = n * d_OSR; |
| 168 | output[n] = 0; |
| 169 | |
| 170 | for (int ii = 0; ii < filter_length; ii++) { |
| 171 | if ((a + ii) >= nitems * d_OSR) |
| 172 | break; |
| 173 | |
| 174 | output[n] += input[a + ii] * filter[ii]; |
| 175 | } |
| 176 | } |
| 177 | } |
| 178 | |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 179 | int get_chan_imp_resp(const gr_complex *input, gr_complex *chan_imp_resp, int search_start_pos, int search_stop_pos, |
| 180 | gr_complex *tseq, int tseqlen, float *corr_max) |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 181 | { |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 182 | const int num_search_windows = search_stop_pos - search_start_pos; |
| 183 | const int power_search_window_len = d_chan_imp_length * d_OSR; |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 184 | std::vector<float> window_energy_buffer; |
| 185 | std::vector<float> power_buffer; |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 186 | std::vector<gr_complex> correlation_buffer; |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 187 | |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 188 | power_buffer.reserve(num_search_windows); |
| 189 | correlation_buffer.reserve(num_search_windows); |
| 190 | window_energy_buffer.reserve(num_search_windows); |
| 191 | |
| 192 | for (int ii = 0; ii < num_search_windows; ii++) { |
| 193 | gr_complex correlation = correlate_sequence(tseq, tseqlen, &input[search_start_pos + ii]); |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 194 | correlation_buffer.push_back(correlation); |
| 195 | power_buffer.push_back(std::pow(abs(correlation), 2)); |
| 196 | } |
| 197 | |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 198 | /* Compute window energies */ |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 199 | float windowSum = 0; |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 200 | |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 201 | // first window |
| 202 | for (int i = 0; i < power_search_window_len; i++) { |
| 203 | windowSum += power_buffer[i]; |
| 204 | } |
| 205 | window_energy_buffer.push_back(windowSum); |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 206 | |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 207 | // slide windows |
| 208 | for (int i = power_search_window_len; i < num_search_windows; i++) { |
| 209 | windowSum += power_buffer[i] - power_buffer[i - power_search_window_len]; |
| 210 | window_energy_buffer.push_back(windowSum); |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 211 | } |
| 212 | |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 213 | int strongest_window_nr = std::max_element(window_energy_buffer.begin(), window_energy_buffer.end()) - |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 214 | window_energy_buffer.begin(); |
| 215 | |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 216 | float max_correlation = 0; |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 217 | for (int ii = 0; ii < power_search_window_len; ii++) { |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 218 | gr_complex correlation = correlation_buffer[strongest_window_nr + ii]; |
| 219 | if (abs(correlation) > max_correlation) |
| 220 | max_correlation = abs(correlation); |
| 221 | chan_imp_resp[ii] = correlation; |
| 222 | } |
| 223 | |
| 224 | *corr_max = max_correlation; |
| 225 | |
| 226 | /** |
| 227 | * Compute first sample position, which corresponds |
| 228 | * to the first sample of the impulse response |
| 229 | */ |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 230 | return search_start_pos + strongest_window_nr; |
| 231 | } |
| 232 | |
| 233 | /* |
| 234 | 8 ext tail bits |
| 235 | 41 sync seq |
| 236 | 36 encrypted bits |
| 237 | 3 tail bits |
| 238 | 68.25 extended tail bits (!) |
| 239 | |
| 240 | center at 8+5 (actually known tb -> known isi, start at 8?) FIXME |
| 241 | */ |
| 242 | int get_access_imp_resp(const gr_complex *input, gr_complex *chan_imp_resp, float *corr_max, int max_delay) |
| 243 | { |
| 244 | const int search_center = 8 + 5; |
| 245 | const int search_start_pos = (search_center - 5) * d_OSR + 1; |
| 246 | const int search_stop_pos = (search_center + 5 + d_chan_imp_length + max_delay) * d_OSR; |
| 247 | const auto tseq = &d_acc_training_seq[TRAIN_BEGINNING]; |
| 248 | const auto tseqlen = N_ACCESS_BITS - (2 * TRAIN_BEGINNING); |
| 249 | return get_chan_imp_resp(input, chan_imp_resp, search_start_pos, search_stop_pos, tseq, tseqlen, corr_max) - |
| 250 | search_center * d_OSR; |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 251 | } |
| 252 | |
| 253 | /* |
| 254 | |
| 255 | 3 + 57 + 1 + 26 + 1 + 57 + 3 + 8.25 |
| 256 | |
| 257 | search center = 3 + 57 + 1 + 5 (due to tsc 5+16+5 split) |
| 258 | this is +-5 samples around (+5 beginning) of truncated t16 tsc |
| 259 | |
| 260 | */ |
| 261 | int get_norm_chan_imp_resp(const gr_complex *input, gr_complex *chan_imp_resp, float *corr_max, int bcc) |
| 262 | { |
| 263 | const int search_center = TRAIN_POS; |
| 264 | const int search_start_pos = (search_center - 5) * d_OSR + 1; |
| 265 | const int search_stop_pos = (search_center + 5 + d_chan_imp_length) * d_OSR; |
| 266 | const auto tseq = &d_norm_training_seq[bcc][TRAIN_BEGINNING]; |
| 267 | const auto tseqlen = N_TRAIN_BITS - (2 * TRAIN_BEGINNING); |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 268 | return get_chan_imp_resp(input, chan_imp_resp, search_start_pos, search_stop_pos, tseq, tseqlen, corr_max) - |
| 269 | search_center * d_OSR; |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 270 | } |
| 271 | |
| 272 | /* |
| 273 | |
| 274 | 3 tail | 39 data | 64 tsc | 39 data | 3 tail | 8.25 guard |
| 275 | start 3+39 - 10 |
| 276 | end 3+39 + SYNC_SEARCH_RANGE |
| 277 | |
| 278 | */ |
| 279 | int get_sch_chan_imp_resp(const gr_complex *input, gr_complex *chan_imp_resp) |
| 280 | { |
| 281 | const int search_center = SYNC_POS + TRAIN_BEGINNING; |
| 282 | const int search_start_pos = (search_center - 10) * d_OSR; |
| 283 | const int search_stop_pos = (search_center + SYNC_SEARCH_RANGE) * d_OSR; |
| 284 | const auto tseq = &d_sch_training_seq[TRAIN_BEGINNING]; |
| 285 | const auto tseqlen = N_SYNC_BITS - (2 * TRAIN_BEGINNING); |
| 286 | |
| 287 | // strongest_window_nr + chan_imp_resp_center + SYNC_POS *d_OSR - 48 * d_OSR - 2 * d_OSR + 2 ; |
| 288 | float corr_max; |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 289 | return get_chan_imp_resp(input, chan_imp_resp, search_start_pos, search_stop_pos, tseq, tseqlen, &corr_max) - |
| 290 | search_center * d_OSR; |
| 291 | ; |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 292 | } |
| 293 | |
| 294 | int get_sch_buffer_chan_imp_resp(const gr_complex *input, gr_complex *chan_imp_resp, unsigned int len, float *corr_max) |
| 295 | { |
| 296 | const auto tseqlen = N_SYNC_BITS - (2 * TRAIN_BEGINNING); |
| 297 | const int search_center = SYNC_POS + TRAIN_BEGINNING; |
| 298 | const int search_start_pos = 0; |
| 299 | // FIXME: proper end offset |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 300 | const int search_stop_pos = len - (N_SYNC_BITS * 8); |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 301 | auto tseq = &d_sch_training_seq[TRAIN_BEGINNING]; |
| 302 | |
Eric | a5a2275 | 2023-01-09 22:46:41 +0100 | [diff] [blame] | 303 | return get_chan_imp_resp(input, chan_imp_resp, search_start_pos, search_stop_pos, tseq, tseqlen, corr_max) - |
| 304 | search_center * d_OSR; |
Eric | ac726b1 | 2022-11-28 18:57:58 +0100 | [diff] [blame] | 305 | } |