blob: 2182550a64e0cdb1663d7352a35a7698e75e1a34 [file] [log] [blame]
dburgessb3a0ca42011-10-12 07:44:40 +00001/*
kurtis.heimerla198d452011-11-26 03:19:28 +00002* Copyright 2008, 2011 Free Software Foundation, Inc.
dburgessb3a0ca42011-10-12 07:44:40 +00003*
4* This software is distributed under the terms of the GNU Affero Public License.
5* See the COPYING file in the main directory for details.
6*
7* This use of this software may be subject to additional restrictions.
8* See the LEGAL file in the main directory for details.
9
10 This program is free software: you can redistribute it and/or modify
11 it under the terms of the GNU Affero General Public License as published by
12 the Free Software Foundation, either version 3 of the License, or
13 (at your option) any later version.
14
15 This program is distributed in the hope that it will be useful,
16 but WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 GNU Affero General Public License for more details.
19
20 You should have received a copy of the GNU Affero General Public License
21 along with this program. If not, see <http://www.gnu.org/licenses/>.
22
23*/
24
Thomas Tsou7e4e5362013-10-30 21:18:55 -040025#ifdef HAVE_CONFIG_H
26#include "config.h"
27#endif
28
dburgessb3a0ca42011-10-12 07:44:40 +000029#include "sigProcLib.h"
30#include "GSMCommon.h"
Alexander Chemeris954b1182015-06-04 15:39:41 -040031#include "Logger.h"
Tom Tsoud3253432016-03-06 03:08:01 -080032#include "Resampler.h"
dburgessb3a0ca42011-10-12 07:44:40 +000033
Thomas Tsou3eaae802013-08-20 19:31:14 -040034extern "C" {
35#include "convolve.h"
Thomas Tsou7e4e5362013-10-30 21:18:55 -040036#include "scale.h"
Thomas Tsou0a3dc4c2013-11-09 02:29:55 -050037#include "mult.h"
Thomas Tsou3eaae802013-08-20 19:31:14 -040038}
39
Thomas Tsou7e4e5362013-10-30 21:18:55 -040040using namespace GSM;
41
Thomas Tsouf79c4d02013-11-09 15:51:56 -060042#define TABLESIZE 1024
43#define DELAYFILTS 64
dburgessb3a0ca42011-10-12 07:44:40 +000044
Tom Tsou577cd022015-05-18 13:57:54 -070045/* Clipping detection threshold */
46#define CLIP_THRESH 30000.0f
47
dburgessb3a0ca42011-10-12 07:44:40 +000048/** Lookup tables for trigonometric approximation */
49float cosTable[TABLESIZE+1]; // add 1 element for wrap around
50float sinTable[TABLESIZE+1];
Thomas Tsou0e0e1f42013-11-09 22:08:51 -050051float sincTable[TABLESIZE+1];
dburgessb3a0ca42011-10-12 07:44:40 +000052
53/** Constants */
54static const float M_PI_F = (float)M_PI;
55static const float M_2PI_F = (float)(2.0*M_PI);
56static const float M_1_2PI_F = 1/M_2PI_F;
57
Thomas Tsouc1f7c422013-10-11 13:49:55 -040058/* Precomputed rotation vectors */
Tom Tsou2079a3c2016-03-06 00:58:56 -080059static signalVector *GMSKRotation4 = NULL;
60static signalVector *GMSKReverseRotation4 = NULL;
Thomas Tsouc1f7c422013-10-11 13:49:55 -040061static signalVector *GMSKRotation1 = NULL;
62static signalVector *GMSKReverseRotation1 = NULL;
dburgessb3a0ca42011-10-12 07:44:40 +000063
Thomas Tsouf79c4d02013-11-09 15:51:56 -060064/* Precomputed fractional delay filters */
65static signalVector *delayFilters[DELAYFILTS];
66
Tom Tsoud3253432016-03-06 03:08:01 -080067static Complex<float> psk8_table[8] = {
68 Complex<float>(-0.70710678, 0.70710678),
69 Complex<float>( 0.0, -1.0),
70 Complex<float>( 0.0, 1.0),
71 Complex<float>( 0.70710678, -0.70710678),
72 Complex<float>(-1.0, 0.0),
73 Complex<float>(-0.70710678, -0.70710678),
74 Complex<float>( 0.70710678, 0.70710678),
75 Complex<float>( 1.0, 0.0),
76};
77
78/* Downsampling filterbank - 4 SPS to 1 SPS */
79#define DOWNSAMPLE_IN_LEN 624
80#define DOWNSAMPLE_OUT_LEN 156
81
82static Resampler *dnsampler = NULL;
83static signalVector *dnsampler_in = NULL;
84
Thomas Tsoue5dcfc42013-08-20 16:27:12 -040085/*
Thomas Tsou3eaae802013-08-20 19:31:14 -040086 * RACH and midamble correlation waveforms. Store the buffer separately
87 * because we need to allocate it explicitly outside of the signal vector
88 * constructor. This is because C++ (prior to C++11) is unable to natively
89 * perform 16-byte memory alignment required by many SSE instructions.
Thomas Tsoue5dcfc42013-08-20 16:27:12 -040090 */
91struct CorrelationSequence {
Thomas Tsou3eaae802013-08-20 19:31:14 -040092 CorrelationSequence() : sequence(NULL), buffer(NULL)
Thomas Tsoue5dcfc42013-08-20 16:27:12 -040093 {
94 }
95
96 ~CorrelationSequence()
97 {
98 delete sequence;
Thomas Tsou3eaae802013-08-20 19:31:14 -040099 free(buffer);
Thomas Tsoue5dcfc42013-08-20 16:27:12 -0400100 }
101
dburgessb3a0ca42011-10-12 07:44:40 +0000102 signalVector *sequence;
Thomas Tsou3eaae802013-08-20 19:31:14 -0400103 void *buffer;
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400104 float toa;
dburgessb3a0ca42011-10-12 07:44:40 +0000105 complex gain;
Thomas Tsoue5dcfc42013-08-20 16:27:12 -0400106};
dburgessb3a0ca42011-10-12 07:44:40 +0000107
Thomas Tsou83e06892013-08-20 16:10:01 -0400108/*
Thomas Tsou3eaae802013-08-20 19:31:14 -0400109 * Gaussian and empty modulation pulses. Like the correlation sequences,
110 * store the runtime (Gaussian) buffer separately because of needed alignment
111 * for SSE instructions.
Thomas Tsou83e06892013-08-20 16:10:01 -0400112 */
113struct PulseSequence {
Tom Tsoud3253432016-03-06 03:08:01 -0800114 PulseSequence() : c0(NULL), c1(NULL), c0_inv(NULL), empty(NULL),
115 c0_buffer(NULL), c1_buffer(NULL), c0_inv_buffer(NULL)
Thomas Tsou83e06892013-08-20 16:10:01 -0400116 {
117 }
118
119 ~PulseSequence()
120 {
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800121 delete c0;
122 delete c1;
Tom Tsoud3253432016-03-06 03:08:01 -0800123 delete c0_inv;
Thomas Tsou83e06892013-08-20 16:10:01 -0400124 delete empty;
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800125 free(c0_buffer);
126 free(c1_buffer);
Thomas Tsou83e06892013-08-20 16:10:01 -0400127 }
128
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800129 signalVector *c0;
130 signalVector *c1;
Tom Tsoud3253432016-03-06 03:08:01 -0800131 signalVector *c0_inv;
Thomas Tsou83e06892013-08-20 16:10:01 -0400132 signalVector *empty;
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800133 void *c0_buffer;
134 void *c1_buffer;
Tom Tsoud3253432016-03-06 03:08:01 -0800135 void *c0_inv_buffer;
Thomas Tsou83e06892013-08-20 16:10:01 -0400136};
137
Tom Tsoud3253432016-03-06 03:08:01 -0800138static CorrelationSequence *gMidambles[] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
139static CorrelationSequence *gEdgeMidambles[] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
140static CorrelationSequence *gRACHSequence = NULL;
141static PulseSequence *GSMPulse1 = NULL;
142static PulseSequence *GSMPulse4 = NULL;
dburgessb3a0ca42011-10-12 07:44:40 +0000143
Thomas Tsoue5dcfc42013-08-20 16:27:12 -0400144void sigProcLibDestroy()
145{
dburgessb3a0ca42011-10-12 07:44:40 +0000146 for (int i = 0; i < 8; i++) {
Thomas Tsoue5dcfc42013-08-20 16:27:12 -0400147 delete gMidambles[i];
Tom Tsoud3253432016-03-06 03:08:01 -0800148 delete gEdgeMidambles[i];
Thomas Tsoue5dcfc42013-08-20 16:27:12 -0400149 gMidambles[i] = NULL;
Tom Tsoud3253432016-03-06 03:08:01 -0800150 gEdgeMidambles[i] = NULL;
dburgessb3a0ca42011-10-12 07:44:40 +0000151 }
Thomas Tsoue5dcfc42013-08-20 16:27:12 -0400152
Thomas Tsouf79c4d02013-11-09 15:51:56 -0600153 for (int i = 0; i < DELAYFILTS; i++) {
154 delete delayFilters[i];
155 delayFilters[i] = NULL;
156 }
157
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400158 delete GMSKRotation1;
159 delete GMSKReverseRotation1;
Tom Tsou2079a3c2016-03-06 00:58:56 -0800160 delete GMSKRotation4;
161 delete GMSKReverseRotation4;
Thomas Tsoue5dcfc42013-08-20 16:27:12 -0400162 delete gRACHSequence;
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400163 delete GSMPulse1;
Tom Tsou2079a3c2016-03-06 00:58:56 -0800164 delete GSMPulse4;
Tom Tsoud3253432016-03-06 03:08:01 -0800165 delete dnsampler;
166 delete dnsampler_in;
Thomas Tsoue5dcfc42013-08-20 16:27:12 -0400167
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400168 GMSKRotation1 = NULL;
Tom Tsou2079a3c2016-03-06 00:58:56 -0800169 GMSKRotation4 = NULL;
170 GMSKReverseRotation4 = NULL;
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400171 GMSKReverseRotation1 = NULL;
Thomas Tsoue5dcfc42013-08-20 16:27:12 -0400172 gRACHSequence = NULL;
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400173 GSMPulse1 = NULL;
Tom Tsou2079a3c2016-03-06 00:58:56 -0800174 GSMPulse4 = NULL;
dburgessb3a0ca42011-10-12 07:44:40 +0000175}
176
dburgessb3a0ca42011-10-12 07:44:40 +0000177// dB relative to 1.0.
178// if > 1.0, then return 0 dB
179float dB(float x) {
180
181 float arg = 1.0F;
182 float dB = 0.0F;
183
184 if (x >= 1.0F) return 0.0F;
185 if (x <= 0.0F) return -200.0F;
186
187 float prevArg = arg;
188 float prevdB = dB;
189 float stepSize = 16.0F;
190 float dBstepSize = 12.0F;
191 while (stepSize > 1.0F) {
192 do {
193 prevArg = arg;
194 prevdB = dB;
195 arg /= stepSize;
196 dB -= dBstepSize;
197 } while (arg > x);
198 arg = prevArg;
199 dB = prevdB;
200 stepSize *= 0.5F;
201 dBstepSize -= 3.0F;
202 }
203 return ((arg-x)*(dB-3.0F) + (x-arg*0.5F)*dB)/(arg - arg*0.5F);
204
205}
206
207// 10^(-dB/10), inverse of dB func.
208float dBinv(float x) {
209
210 float arg = 1.0F;
211 float dB = 0.0F;
212
213 if (x >= 0.0F) return 1.0F;
214 if (x <= -200.0F) return 0.0F;
215
216 float prevArg = arg;
217 float prevdB = dB;
218 float stepSize = 16.0F;
219 float dBstepSize = 12.0F;
220 while (stepSize > 1.0F) {
221 do {
222 prevArg = arg;
223 prevdB = dB;
224 arg /= stepSize;
225 dB -= dBstepSize;
226 } while (dB > x);
227 arg = prevArg;
228 dB = prevdB;
229 stepSize *= 0.5F;
230 dBstepSize -= 3.0F;
231 }
232
233 return ((dB-x)*(arg*0.5F)+(x-(dB-3.0F))*(arg))/3.0F;
234
235}
236
237float vectorNorm2(const signalVector &x)
238{
239 signalVector::const_iterator xPtr = x.begin();
240 float Energy = 0.0;
241 for (;xPtr != x.end();xPtr++) {
242 Energy += xPtr->norm2();
243 }
244 return Energy;
245}
246
247
248float vectorPower(const signalVector &x)
249{
250 return vectorNorm2(x)/x.size();
251}
252
253/** compute cosine via lookup table */
254float cosLookup(const float x)
255{
256 float arg = x*M_1_2PI_F;
257 while (arg > 1.0F) arg -= 1.0F;
258 while (arg < 0.0F) arg += 1.0F;
259
260 const float argT = arg*((float)TABLESIZE);
261 const int argI = (int)argT;
262 const float delta = argT-argI;
263 const float iDelta = 1.0F-delta;
264 return iDelta*cosTable[argI] + delta*cosTable[argI+1];
265}
266
267/** compute sine via lookup table */
268float sinLookup(const float x)
269{
270 float arg = x*M_1_2PI_F;
271 while (arg > 1.0F) arg -= 1.0F;
272 while (arg < 0.0F) arg += 1.0F;
273
274 const float argT = arg*((float)TABLESIZE);
275 const int argI = (int)argT;
276 const float delta = argT-argI;
277 const float iDelta = 1.0F-delta;
278 return iDelta*sinTable[argI] + delta*sinTable[argI+1];
279}
280
281
282/** compute e^(-jx) via lookup table. */
Tom Tsou2079a3c2016-03-06 00:58:56 -0800283static complex expjLookup(float x)
dburgessb3a0ca42011-10-12 07:44:40 +0000284{
285 float arg = x*M_1_2PI_F;
286 while (arg > 1.0F) arg -= 1.0F;
287 while (arg < 0.0F) arg += 1.0F;
288
289 const float argT = arg*((float)TABLESIZE);
290 const int argI = (int)argT;
291 const float delta = argT-argI;
292 const float iDelta = 1.0F-delta;
293 return complex(iDelta*cosTable[argI] + delta*cosTable[argI+1],
294 iDelta*sinTable[argI] + delta*sinTable[argI+1]);
295}
296
297/** Library setup functions */
Tom Tsou2079a3c2016-03-06 00:58:56 -0800298static void initTrigTables() {
dburgessb3a0ca42011-10-12 07:44:40 +0000299 for (int i = 0; i < TABLESIZE+1; i++) {
300 cosTable[i] = cos(2.0*M_PI*i/TABLESIZE);
301 sinTable[i] = sin(2.0*M_PI*i/TABLESIZE);
302 }
303}
304
Tom Tsou2079a3c2016-03-06 00:58:56 -0800305/*
306 * Initialize 4 sps and 1 sps rotation tables
307 */
308static void initGMSKRotationTables()
Thomas Tsoud24cc2c2013-08-20 15:41:45 -0400309{
Tom Tsou2079a3c2016-03-06 00:58:56 -0800310 size_t len1 = 157, len4 = 625;
311
312 GMSKRotation4 = new signalVector(len4);
313 GMSKReverseRotation4 = new signalVector(len4);
314 signalVector::iterator rotPtr = GMSKRotation4->begin();
315 signalVector::iterator revPtr = GMSKReverseRotation4->begin();
dburgessb3a0ca42011-10-12 07:44:40 +0000316 float phase = 0.0;
Tom Tsou2079a3c2016-03-06 00:58:56 -0800317 while (rotPtr != GMSKRotation4->end()) {
dburgessb3a0ca42011-10-12 07:44:40 +0000318 *rotPtr++ = expjLookup(phase);
319 *revPtr++ = expjLookup(-phase);
Tom Tsou2079a3c2016-03-06 00:58:56 -0800320 phase += M_PI_F / 2.0F / 4.0;
dburgessb3a0ca42011-10-12 07:44:40 +0000321 }
dburgessb3a0ca42011-10-12 07:44:40 +0000322
Tom Tsou2079a3c2016-03-06 00:58:56 -0800323 GMSKRotation1 = new signalVector(len1);
324 GMSKReverseRotation1 = new signalVector(len1);
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400325 rotPtr = GMSKRotation1->begin();
326 revPtr = GMSKReverseRotation1->begin();
327 phase = 0.0;
328 while (rotPtr != GMSKRotation1->end()) {
329 *rotPtr++ = expjLookup(phase);
330 *revPtr++ = expjLookup(-phase);
331 phase += M_PI_F / 2.0F;
Thomas Tsoue57004d2013-08-20 18:55:33 -0400332 }
dburgessb3a0ca42011-10-12 07:44:40 +0000333}
334
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400335static void GMSKRotate(signalVector &x, int sps)
336{
Thomas Tsou0a3dc4c2013-11-09 02:29:55 -0500337#if HAVE_NEON
338 size_t len;
339 signalVector *a, *b, *out;
340
341 a = &x;
342 out = &x;
343 len = out->size();
344
345 if (len == 157)
346 len--;
347
348 if (sps == 1)
349 b = GMSKRotation1;
350 else
Tom Tsou2079a3c2016-03-06 00:58:56 -0800351 b = GMSKRotation4;
Thomas Tsou0a3dc4c2013-11-09 02:29:55 -0500352
353 mul_complex((float *) out->begin(),
354 (float *) a->begin(),
355 (float *) b->begin(), len);
356#else
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400357 signalVector::iterator rotPtr, xPtr = x.begin();
358
359 if (sps == 1)
360 rotPtr = GMSKRotation1->begin();
361 else
Tom Tsou2079a3c2016-03-06 00:58:56 -0800362 rotPtr = GMSKRotation4->begin();
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400363
Thomas Tsou20eb6d62013-11-09 14:30:41 -0500364 if (x.isReal()) {
dburgessb3a0ca42011-10-12 07:44:40 +0000365 while (xPtr < x.end()) {
366 *xPtr = *rotPtr++ * (xPtr->real());
367 xPtr++;
368 }
369 }
370 else {
371 while (xPtr < x.end()) {
372 *xPtr = *rotPtr++ * (*xPtr);
373 xPtr++;
374 }
375 }
Thomas Tsou0a3dc4c2013-11-09 02:29:55 -0500376#endif
dburgessb3a0ca42011-10-12 07:44:40 +0000377}
378
Tom Tsou2079a3c2016-03-06 00:58:56 -0800379static bool GMSKReverseRotate(signalVector &x, int sps)
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400380{
381 signalVector::iterator rotPtr, xPtr= x.begin();
382
383 if (sps == 1)
384 rotPtr = GMSKReverseRotation1->begin();
Tom Tsou2079a3c2016-03-06 00:58:56 -0800385 else if (sps == 4)
386 rotPtr = GMSKReverseRotation4->begin();
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400387 else
Tom Tsou2079a3c2016-03-06 00:58:56 -0800388 return false;
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400389
Thomas Tsou20eb6d62013-11-09 14:30:41 -0500390 if (x.isReal()) {
dburgessb3a0ca42011-10-12 07:44:40 +0000391 while (xPtr < x.end()) {
392 *xPtr = *rotPtr++ * (xPtr->real());
393 xPtr++;
394 }
395 }
396 else {
397 while (xPtr < x.end()) {
398 *xPtr = *rotPtr++ * (*xPtr);
399 xPtr++;
400 }
401 }
Tom Tsou2079a3c2016-03-06 00:58:56 -0800402
403 return true;
dburgessb3a0ca42011-10-12 07:44:40 +0000404}
405
Thomas Tsou3eaae802013-08-20 19:31:14 -0400406signalVector *convolve(const signalVector *x,
407 const signalVector *h,
408 signalVector *y,
Thomas Tsou3f32ab52013-11-15 16:32:54 -0500409 ConvType spanType, size_t start,
410 size_t len, size_t step, int offset)
dburgessb3a0ca42011-10-12 07:44:40 +0000411{
Thomas Tsou3f32ab52013-11-15 16:32:54 -0500412 int rc;
413 size_t head = 0, tail = 0;
Thomas Tsou3eaae802013-08-20 19:31:14 -0400414 bool alloc = false, append = false;
415 const signalVector *_x = NULL;
dburgessb3a0ca42011-10-12 07:44:40 +0000416
Thomas Tsou3eaae802013-08-20 19:31:14 -0400417 if (!x || !h)
dburgessb3a0ca42011-10-12 07:44:40 +0000418 return NULL;
419
Thomas Tsou3eaae802013-08-20 19:31:14 -0400420 switch (spanType) {
421 case START_ONLY:
422 start = 0;
Thomas Tsou6f4906e2013-11-09 02:40:18 -0500423 head = h->size() - 1;
Thomas Tsou3eaae802013-08-20 19:31:14 -0400424 len = x->size();
Thomas Tsou6f4906e2013-11-09 02:40:18 -0500425
Thomas Tsou20eb6d62013-11-09 14:30:41 -0500426 if (x->getStart() < head)
Thomas Tsou6f4906e2013-11-09 02:40:18 -0500427 append = true;
dburgessb3a0ca42011-10-12 07:44:40 +0000428 break;
Thomas Tsou3eaae802013-08-20 19:31:14 -0400429 case NO_DELAY:
430 start = h->size() / 2;
431 head = start;
432 tail = start;
433 len = x->size();
434 append = true;
435 break;
436 case CUSTOM:
437 if (start < h->size() - 1) {
438 head = h->size() - start;
439 append = true;
440 }
441 if (start + len > x->size()) {
442 tail = start + len - x->size();
443 append = true;
dburgessb3a0ca42011-10-12 07:44:40 +0000444 }
445 break;
446 default:
447 return NULL;
dburgessb3a0ca42011-10-12 07:44:40 +0000448 }
dburgessb3a0ca42011-10-12 07:44:40 +0000449
Thomas Tsou3eaae802013-08-20 19:31:14 -0400450 /*
451 * Error if the output vector is too small. Create the output vector
452 * if the pointer is NULL.
453 */
454 if (y && (len > y->size()))
455 return NULL;
456 if (!y) {
457 y = new signalVector(len);
458 alloc = true;
459 }
460
461 /* Prepend or post-pend the input vector if the parameters require it */
462 if (append)
463 _x = new signalVector(*x, head, tail);
464 else
465 _x = x;
466
467 /*
468 * Four convovle types:
469 * 1. Complex-Real (aligned)
470 * 2. Complex-Complex (aligned)
471 * 3. Complex-Real (!aligned)
472 * 4. Complex-Complex (!aligned)
473 */
Thomas Tsou20eb6d62013-11-09 14:30:41 -0500474 if (h->isReal() && h->isAligned()) {
Thomas Tsou3eaae802013-08-20 19:31:14 -0400475 rc = convolve_real((float *) _x->begin(), _x->size(),
476 (float *) h->begin(), h->size(),
477 (float *) y->begin(), y->size(),
478 start, len, step, offset);
Thomas Tsou20eb6d62013-11-09 14:30:41 -0500479 } else if (!h->isReal() && h->isAligned()) {
Thomas Tsou3eaae802013-08-20 19:31:14 -0400480 rc = convolve_complex((float *) _x->begin(), _x->size(),
481 (float *) h->begin(), h->size(),
482 (float *) y->begin(), y->size(),
483 start, len, step, offset);
Thomas Tsou20eb6d62013-11-09 14:30:41 -0500484 } else if (h->isReal() && !h->isAligned()) {
Thomas Tsou3eaae802013-08-20 19:31:14 -0400485 rc = base_convolve_real((float *) _x->begin(), _x->size(),
486 (float *) h->begin(), h->size(),
487 (float *) y->begin(), y->size(),
488 start, len, step, offset);
Thomas Tsou20eb6d62013-11-09 14:30:41 -0500489 } else if (!h->isReal() && !h->isAligned()) {
Thomas Tsou3eaae802013-08-20 19:31:14 -0400490 rc = base_convolve_complex((float *) _x->begin(), _x->size(),
491 (float *) h->begin(), h->size(),
492 (float *) y->begin(), y->size(),
493 start, len, step, offset);
494 } else {
495 rc = -1;
496 }
497
498 if (append)
499 delete _x;
500
501 if (rc < 0) {
502 if (alloc)
503 delete y;
504 return NULL;
505 }
506
507 return y;
508}
dburgessb3a0ca42011-10-12 07:44:40 +0000509
Tom Tsoud3253432016-03-06 03:08:01 -0800510/*
511 * Generate static EDGE linear equalizer. This equalizer is not adaptive.
512 * Filter taps are generated from the inverted 1 SPS impulse response of
513 * the EDGE pulse shape captured after the downsampling filter.
514 */
515static bool generateInvertC0Pulse(PulseSequence *pulse)
516{
517 if (!pulse)
518 return false;
519
520 pulse->c0_inv_buffer = convolve_h_alloc(5);
521 pulse->c0_inv = new signalVector((complex *) pulse->c0_inv_buffer, 0, 5);
522 pulse->c0_inv->isReal(true);
523 pulse->c0_inv->setAligned(false);
524
525 signalVector::iterator xP = pulse->c0_inv->begin();
526 *xP++ = 0.15884;
527 *xP++ = -0.43176;
528 *xP++ = 1.00000;
529 *xP++ = -0.42608;
530 *xP++ = 0.14882;
531
532 return true;
533}
534
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400535static bool generateC1Pulse(int sps, PulseSequence *pulse)
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800536{
537 int len;
538
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400539 if (!pulse)
540 return false;
541
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800542 switch (sps) {
543 case 4:
544 len = 8;
545 break;
546 default:
547 return false;
548 }
549
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400550 pulse->c1_buffer = convolve_h_alloc(len);
551 pulse->c1 = new signalVector((complex *)
552 pulse->c1_buffer, 0, len);
Thomas Tsou20eb6d62013-11-09 14:30:41 -0500553 pulse->c1->isReal(true);
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800554
555 /* Enable alignment for SSE usage */
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400556 pulse->c1->setAligned(true);
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800557
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400558 signalVector::iterator xP = pulse->c1->begin();
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800559
560 switch (sps) {
561 case 4:
562 /* BT = 0.30 */
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400563 *xP++ = 0.0;
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800564 *xP++ = 8.16373112e-03;
565 *xP++ = 2.84385729e-02;
566 *xP++ = 5.64158904e-02;
567 *xP++ = 7.05463553e-02;
568 *xP++ = 5.64158904e-02;
569 *xP++ = 2.84385729e-02;
570 *xP++ = 8.16373112e-03;
571 }
572
573 return true;
574}
575
Tom Tsou2079a3c2016-03-06 00:58:56 -0800576static PulseSequence *generateGSMPulse(int sps)
dburgessb3a0ca42011-10-12 07:44:40 +0000577{
Thomas Tsou83e06892013-08-20 16:10:01 -0400578 int len;
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800579 float arg, avg, center;
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400580 PulseSequence *pulse;
Thomas Tsou83e06892013-08-20 16:10:01 -0400581
Tom Tsoud3253432016-03-06 03:08:01 -0800582 if ((sps != 1) && (sps != 4))
583 return NULL;
584
Thomas Tsou83e06892013-08-20 16:10:01 -0400585 /* Store a single tap filter used for correlation sequence generation */
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400586 pulse = new PulseSequence();
587 pulse->empty = new signalVector(1);
Thomas Tsou20eb6d62013-11-09 14:30:41 -0500588 pulse->empty->isReal(true);
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400589 *(pulse->empty->begin()) = 1.0f;
Thomas Tsou83e06892013-08-20 16:10:01 -0400590
Thomas Tsou9ccd9f22013-08-21 13:59:52 -0400591 /*
592 * For 4 samples-per-symbol use a precomputed single pulse Laurent
593 * approximation. This should yields below 2 degrees of phase error at
594 * the modulator output. Use the existing pulse approximation for all
595 * other oversampling factors.
596 */
597 switch (sps) {
598 case 4:
599 len = 16;
600 break;
Tom Tsoud3253432016-03-06 03:08:01 -0800601 case 1:
Thomas Tsou9ccd9f22013-08-21 13:59:52 -0400602 default:
Tom Tsou2079a3c2016-03-06 00:58:56 -0800603 len = 4;
Thomas Tsou9ccd9f22013-08-21 13:59:52 -0400604 }
Thomas Tsou3eaae802013-08-20 19:31:14 -0400605
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400606 pulse->c0_buffer = convolve_h_alloc(len);
607 pulse->c0 = new signalVector((complex *) pulse->c0_buffer, 0, len);
Thomas Tsou20eb6d62013-11-09 14:30:41 -0500608 pulse->c0->isReal(true);
Thomas Tsou3eaae802013-08-20 19:31:14 -0400609
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800610 /* Enable alingnment for SSE usage */
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400611 pulse->c0->setAligned(true);
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800612
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400613 signalVector::iterator xP = pulse->c0->begin();
Thomas Tsou83e06892013-08-20 16:10:01 -0400614
Thomas Tsou9ccd9f22013-08-21 13:59:52 -0400615 if (sps == 4) {
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800616 *xP++ = 0.0;
Thomas Tsou9ccd9f22013-08-21 13:59:52 -0400617 *xP++ = 4.46348606e-03;
618 *xP++ = 2.84385729e-02;
619 *xP++ = 1.03184855e-01;
620 *xP++ = 2.56065552e-01;
621 *xP++ = 4.76375085e-01;
622 *xP++ = 7.05961177e-01;
623 *xP++ = 8.71291644e-01;
624 *xP++ = 9.29453645e-01;
625 *xP++ = 8.71291644e-01;
626 *xP++ = 7.05961177e-01;
627 *xP++ = 4.76375085e-01;
628 *xP++ = 2.56065552e-01;
629 *xP++ = 1.03184855e-01;
630 *xP++ = 2.84385729e-02;
631 *xP++ = 4.46348606e-03;
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400632 generateC1Pulse(sps, pulse);
Thomas Tsou9ccd9f22013-08-21 13:59:52 -0400633 } else {
634 center = (float) (len - 1.0) / 2.0;
Thomas Tsou83e06892013-08-20 16:10:01 -0400635
Thomas Tsou9ccd9f22013-08-21 13:59:52 -0400636 /* GSM pulse approximation */
637 for (int i = 0; i < len; i++) {
638 arg = ((float) i - center) / (float) sps;
639 *xP++ = 0.96 * exp(-1.1380 * arg * arg -
640 0.527 * arg * arg * arg * arg);
641 }
dburgessb3a0ca42011-10-12 07:44:40 +0000642
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400643 avg = sqrtf(vectorNorm2(*pulse->c0) / sps);
644 xP = pulse->c0->begin();
645 for (int i = 0; i < len; i++)
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800646 *xP++ /= avg;
647 }
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400648
Tom Tsoud3253432016-03-06 03:08:01 -0800649 /*
650 * Current form of the EDGE equalization filter non-realizable at 4 SPS.
651 * Load the onto both 1 SPS and 4 SPS objects for convenience. Note that
652 * the EDGE demodulator downsamples to 1 SPS prior to equalization.
653 */
654 generateInvertC0Pulse(pulse);
655
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400656 return pulse;
dburgessb3a0ca42011-10-12 07:44:40 +0000657}
658
659signalVector* frequencyShift(signalVector *y,
660 signalVector *x,
661 float freq,
662 float startPhase,
663 float *finalPhase)
664{
665
666 if (!x) return NULL;
667
668 if (y==NULL) {
669 y = new signalVector(x->size());
Thomas Tsou20eb6d62013-11-09 14:30:41 -0500670 y->isReal(x->isReal());
dburgessb3a0ca42011-10-12 07:44:40 +0000671 if (y==NULL) return NULL;
672 }
673
674 if (y->size() < x->size()) return NULL;
675
676 float phase = startPhase;
677 signalVector::iterator yP = y->begin();
678 signalVector::iterator xPEnd = x->end();
679 signalVector::iterator xP = x->begin();
680
Thomas Tsou20eb6d62013-11-09 14:30:41 -0500681 if (x->isReal()) {
dburgessb3a0ca42011-10-12 07:44:40 +0000682 while (xP < xPEnd) {
683 (*yP++) = expjLookup(phase)*( (xP++)->real() );
684 phase += freq;
685 }
686 }
687 else {
688 while (xP < xPEnd) {
689 (*yP++) = (*xP++)*expjLookup(phase);
690 phase += freq;
Thomas Tsou34bbef72013-11-13 22:58:15 -0500691 if (phase > 2 * M_PI)
692 phase -= 2 * M_PI;
693 else if (phase < -2 * M_PI)
694 phase += 2 * M_PI;
dburgessb3a0ca42011-10-12 07:44:40 +0000695 }
696 }
697
698
699 if (finalPhase) *finalPhase = phase;
700
701 return y;
702}
703
704signalVector* reverseConjugate(signalVector *b)
705{
706 signalVector *tmp = new signalVector(b->size());
Thomas Tsou20eb6d62013-11-09 14:30:41 -0500707 tmp->isReal(b->isReal());
dburgessb3a0ca42011-10-12 07:44:40 +0000708 signalVector::iterator bP = b->begin();
709 signalVector::iterator bPEnd = b->end();
710 signalVector::iterator tmpP = tmp->end()-1;
Thomas Tsou20eb6d62013-11-09 14:30:41 -0500711 if (!b->isReal()) {
dburgessb3a0ca42011-10-12 07:44:40 +0000712 while (bP < bPEnd) {
713 *tmpP-- = bP->conj();
714 bP++;
715 }
716 }
717 else {
718 while (bP < bPEnd) {
719 *tmpP-- = bP->real();
720 bP++;
721 }
722 }
723
724 return tmp;
725}
726
Tom Tsoud3253432016-03-06 03:08:01 -0800727bool vectorSlicer(SoftVector *x)
728{
729 SoftVector::iterator xP = x->begin();
730 SoftVector::iterator xPEnd = x->end();
731 while (xP < xPEnd) {
732 *xP = 0.5 * (*xP + 1.0f);
733 if (*xP > 1.0)
734 *xP = 1.0;
735 if (*xP < 0.0)
736 *xP = 0.0;
737 xP++;
738 }
739 return true;
740}
741
742bool vectorSlicer(signalVector *x)
dburgessb3a0ca42011-10-12 07:44:40 +0000743{
744
745 signalVector::iterator xP = x->begin();
746 signalVector::iterator xPEnd = x->end();
747 while (xP < xPEnd) {
748 *xP = (complex) (0.5*(xP->real()+1.0F));
749 if (xP->real() > 1.0) *xP = 1.0;
750 if (xP->real() < 0.0) *xP = 0.0;
751 xP++;
752 }
753 return true;
754}
Thomas Tsou3eaae802013-08-20 19:31:14 -0400755
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800756static signalVector *rotateBurst(const BitVector &wBurst,
757 int guardPeriodLength, int sps)
758{
759 int burst_len;
760 signalVector *pulse, rotated, *shaped;
761 signalVector::iterator itr;
762
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400763 pulse = GSMPulse1->empty;
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800764 burst_len = sps * (wBurst.size() + guardPeriodLength);
765 rotated = signalVector(burst_len);
766 itr = rotated.begin();
767
768 for (unsigned i = 0; i < wBurst.size(); i++) {
769 *itr = 2.0 * (wBurst[i] & 0x01) - 1.0;
770 itr += sps;
771 }
772
Thomas Tsouc1f7c422013-10-11 13:49:55 -0400773 GMSKRotate(rotated, sps);
Thomas Tsou20eb6d62013-11-09 14:30:41 -0500774 rotated.isReal(false);
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800775
776 /* Dummy filter operation */
777 shaped = convolve(&rotated, pulse, NULL, START_ONLY);
778 if (!shaped)
779 return NULL;
780
781 return shaped;
782}
783
Tom Tsoud3253432016-03-06 03:08:01 -0800784static void rotateBurst2(signalVector &burst, double phase)
785{
786 Complex<float> rot = Complex<float>(cos(phase), sin(phase));
787
788 for (size_t i = 0; i < burst.size(); i++)
789 burst[i] = burst[i] * rot;
790}
791
Tom Tsou4dfd64a2016-03-06 20:31:51 -0800792/*
793 * Ignore the guard length argument in the GMSK modulator interface
794 * because it results in 624/628 sized bursts instead of the preferred
795 * burst length of 625. Only 4 SPS is supported.
796 */
797static signalVector *modulateBurstLaurent(const BitVector &bits)
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800798{
Tom Tsou4dfd64a2016-03-06 20:31:51 -0800799 int burst_len, sps = 4;
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800800 float phase;
Thomas Tsou6f4906e2013-11-09 02:40:18 -0500801 signalVector *c0_pulse, *c1_pulse, *c0_burst;
802 signalVector *c1_burst, *c0_shaped, *c1_shaped;
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800803 signalVector::iterator c0_itr, c1_itr;
804
Tom Tsou2079a3c2016-03-06 00:58:56 -0800805 c0_pulse = GSMPulse4->c0;
806 c1_pulse = GSMPulse4->c1;
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800807
Tom Tsou4dfd64a2016-03-06 20:31:51 -0800808 if (bits.size() > 156)
809 return NULL;
810
811 burst_len = 625;
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800812
Thomas Tsou6f4906e2013-11-09 02:40:18 -0500813 c0_burst = new signalVector(burst_len, c0_pulse->size());
Thomas Tsou20eb6d62013-11-09 14:30:41 -0500814 c0_burst->isReal(true);
Thomas Tsou6f4906e2013-11-09 02:40:18 -0500815 c0_itr = c0_burst->begin();
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800816
Thomas Tsou6f4906e2013-11-09 02:40:18 -0500817 c1_burst = new signalVector(burst_len, c1_pulse->size());
Thomas Tsou20eb6d62013-11-09 14:30:41 -0500818 c1_burst->isReal(true);
Thomas Tsou6f4906e2013-11-09 02:40:18 -0500819 c1_itr = c1_burst->begin();
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800820
821 /* Padded differential start bits */
Alexander Chemeris351fd762015-05-24 20:16:51 -0400822 *c0_itr = 2.0 * (0x01 & 0x01) - 1.0;
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800823 c0_itr += sps;
824
825 /* Main burst bits */
826 for (unsigned i = 0; i < bits.size(); i++) {
827 *c0_itr = 2.0 * (bits[i] & 0x01) - 1.0;
828 c0_itr += sps;
829 }
830
831 /* Padded differential end bits */
832 *c0_itr = 2.0 * (0x01 & 0x01) - 1.0;
833
834 /* Generate C0 phase coefficients */
Thomas Tsou6f4906e2013-11-09 02:40:18 -0500835 GMSKRotate(*c0_burst, sps);
Thomas Tsou20eb6d62013-11-09 14:30:41 -0500836 c0_burst->isReal(false);
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800837
Thomas Tsou6f4906e2013-11-09 02:40:18 -0500838 c0_itr = c0_burst->begin();
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800839 c0_itr += sps * 2;
840 c1_itr += sps * 2;
841
842 /* Start magic */
843 phase = 2.0 * ((0x01 & 0x01) ^ (0x01 & 0x01)) - 1.0;
844 *c1_itr = *c0_itr * Complex<float>(0, phase);
845 c0_itr += sps;
846 c1_itr += sps;
847
848 /* Generate C1 phase coefficients */
849 for (unsigned i = 2; i < bits.size(); i++) {
850 phase = 2.0 * ((bits[i - 1] & 0x01) ^ (bits[i - 2] & 0x01)) - 1.0;
851 *c1_itr = *c0_itr * Complex<float>(0, phase);
852
853 c0_itr += sps;
854 c1_itr += sps;
855 }
856
857 /* End magic */
858 int i = bits.size();
859 phase = 2.0 * ((bits[i-1] & 0x01) ^ (bits[i-2] & 0x01)) - 1.0;
860 *c1_itr = *c0_itr * Complex<float>(0, phase);
861
862 /* Primary (C0) and secondary (C1) pulse shaping */
Thomas Tsou6f4906e2013-11-09 02:40:18 -0500863 c0_shaped = convolve(c0_burst, c0_pulse, NULL, START_ONLY);
864 c1_shaped = convolve(c1_burst, c1_pulse, NULL, START_ONLY);
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800865
866 /* Sum shaped outputs into C0 */
867 c0_itr = c0_shaped->begin();
868 c1_itr = c1_shaped->begin();
869 for (unsigned i = 0; i < c0_shaped->size(); i++ )
870 *c0_itr++ += *c1_itr++;
871
Thomas Tsou6f4906e2013-11-09 02:40:18 -0500872 delete c0_burst;
873 delete c1_burst;
Thomas Tsoua57bc8a2013-09-05 08:16:47 +0800874 delete c1_shaped;
875
876 return c0_shaped;
877}
878
Tom Tsoud3253432016-03-06 03:08:01 -0800879static signalVector *rotateEdgeBurst(const signalVector &symbols, int sps)
880{
881 signalVector *burst;
882 signalVector::iterator burst_itr;
883
884 burst = new signalVector(symbols.size() * sps);
885 burst_itr = burst->begin();
886
887 for (size_t i = 0; i < symbols.size(); i++) {
888 float phase = i * 3.0f * M_PI / 8.0f;
889 Complex<float> rot = Complex<float>(cos(phase), sin(phase));
890
891 *burst_itr = symbols[i] * rot;
892 burst_itr += sps;
893 }
894
895 return burst;
896}
897
898static signalVector *derotateEdgeBurst(const signalVector &symbols, int sps)
899{
900 signalVector *burst;
901 signalVector::iterator burst_itr;
902
903 if (symbols.size() % sps)
904 return NULL;
905
906 burst = new signalVector(symbols.size() / sps);
907 burst_itr = burst->begin();
908
909 for (size_t i = 0; i < burst->size(); i++) {
910 float phase = (float) (i % 16) * 3.0f * M_PI / 8.0f;
911 Complex<float> rot = Complex<float>(cosf(phase), -sinf(phase));
912
913 *burst_itr = symbols[sps * i] * rot;
914 burst_itr++;
915 }
916
917 return burst;
918}
919
920static signalVector *mapEdgeSymbols(const BitVector &bits)
921{
922 if (bits.size() % 3)
923 return NULL;
924
925 signalVector *symbols = new signalVector(bits.size() / 3);
926
927 for (size_t i = 0; i < symbols->size(); i++) {
928 unsigned index = (((unsigned) bits[3 * i + 0] & 0x01) << 0) |
929 (((unsigned) bits[3 * i + 1] & 0x01) << 1) |
930 (((unsigned) bits[3 * i + 2] & 0x01) << 2);
931
932 (*symbols)[i] = psk8_table[index];
933 }
934
935 return symbols;
936}
937
938static signalVector *shapeEdgeBurst(const signalVector &symbols)
939{
940 size_t nsyms, nsamps = 625;
941 signalVector *burst, *shape;
942 signalVector::iterator burst_itr;
943
944 nsyms = symbols.size();
945
946 if (nsyms * 4 > nsamps)
947 nsyms = 156;
948
949 burst = new signalVector(nsamps, GSMPulse4->c0->size());
950 burst_itr = burst->begin();
951
952 for (size_t i = 0; i < nsyms; i++) {
953 float phase = i * 3.0f * M_PI / 8.0f;
954 Complex<float> rot = Complex<float>(cos(phase), sin(phase));
955
956 *burst_itr = symbols[i] * rot;
957 burst_itr += 4;
958 }
959
960 /* Single Gaussian pulse approximation shaping */
961 shape = convolve(burst, GSMPulse4->c0, NULL, START_ONLY);
962 delete burst;
963
964 return shape;
965}
966
967/*
Tom Tsou8ee2f382016-03-06 20:57:34 -0800968 * Generate a random GSM normal burst.
969 */
970signalVector *genRandNormalBurst(int tsc, int sps, int tn)
971{
972 if ((tsc < 0) || (tsc > 7) || (tn < 0) || (tn > 7))
973 return NULL;
974 if ((sps != 1) && (sps != 4))
975 return NULL;
976
977 int i = 0;
978 BitVector *bits = new BitVector(148);
979 signalVector *burst;
980
981 /* Tail bits */
982 for (; i < 4; i++)
983 (*bits)[i] = 0;
984
985 /* Random bits */
986 for (; i < 61; i++)
987 (*bits)[i] = rand() % 2;
988
989 /* Training sequence */
990 for (int n = 0; i < 87; i++, n++)
991 (*bits)[i] = gTrainingSequence[tsc][n];
992
993 /* Random bits */
994 for (; i < 144; i++)
995 (*bits)[i] = rand() % 2;
996
997 /* Tail bits */
998 for (; i < 148; i++)
999 (*bits)[i] = 0;
1000
1001 int guard = 8 + !(tn % 4);
1002 burst = modulateBurst(*bits, guard, sps);
1003 delete bits;
1004
1005 return burst;
1006}
1007
Alexander Chemeris5efe0502016-03-23 17:06:32 +03001008/*
1009 * Generate a random GSM access burst.
1010 */
1011signalVector *genRandAccessBurst(int sps, int tn)
1012{
1013 if ((tn < 0) || (tn > 7))
1014 return NULL;
1015 if ((sps != 1) && (sps != 4))
1016 return NULL;
1017
1018 int i = 0;
1019 BitVector *bits = new BitVector(88);
1020 signalVector *burst;
1021
1022 /* head and synch bits */
1023 for (int n = 0; i < 49; i++, n++)
1024 (*bits)[i] = gRACHBurst[n];
1025
1026 /* Random bits */
1027 for (; i < 85; i++)
1028 (*bits)[i] = rand() % 2;
1029
1030 /* Tail bits */
1031 for (; i < 88; i++)
1032 (*bits)[i] = 0;
1033
1034 int guard = 68 + !(tn % 4);
1035 burst = modulateBurst(*bits, guard, sps);
1036 delete bits;
1037
1038 return burst;
1039}
1040
Tom Tsou8ee2f382016-03-06 20:57:34 -08001041signalVector *generateEmptyBurst(int sps, int tn)
1042{
1043 if ((tn < 0) || (tn > 7))
1044 return NULL;
1045
1046 if (sps == 4)
1047 return new signalVector(625);
1048 else if (sps == 1)
1049 return new signalVector(148 + 8 + !(tn % 4));
1050 else
1051 return NULL;
1052}
1053
1054signalVector *generateDummyBurst(int sps, int tn)
1055{
1056 if (((sps != 1) && (sps != 4)) || (tn < 0) || (tn > 7))
1057 return NULL;
1058
1059 return modulateBurst(gDummyBurst, 8 + !(tn % 4), sps);
1060}
1061
1062/*
Tom Tsoud3253432016-03-06 03:08:01 -08001063 * Generate a random 8-PSK EDGE burst. Only 4 SPS is supported with
1064 * the returned burst being 625 samples in length.
1065 */
1066signalVector *generateEdgeBurst(int tsc)
1067{
1068 int tail = 9 / 3;
1069 int data = 174 / 3;
1070 int train = 78 / 3;
1071
1072 if ((tsc < 0) || (tsc > 7))
1073 return NULL;
1074
1075 signalVector *shape, *burst = new signalVector(148);
1076 const BitVector *midamble = &gEdgeTrainingSequence[tsc];
1077
1078 /* Tail */
1079 int n, i = 0;
1080 for (; i < tail; i++)
1081 (*burst)[i] = psk8_table[7];
1082
1083 /* Body */
1084 for (; i < tail + data; i++)
1085 (*burst)[i] = psk8_table[rand() % 8];
1086
1087 /* TSC */
1088 for (n = 0; i < tail + data + train; i++, n++) {
1089 unsigned index = (((unsigned) (*midamble)[3 * n + 0] & 0x01) << 0) |
1090 (((unsigned) (*midamble)[3 * n + 1] & 0x01) << 1) |
1091 (((unsigned) (*midamble)[3 * n + 2] & 0x01) << 2);
1092
1093 (*burst)[i] = psk8_table[index];
1094 }
1095
1096 /* Body */
1097 for (; i < tail + data + train + data; i++)
1098 (*burst)[i] = psk8_table[rand() % 8];
1099
1100 /* Tail */
1101 for (; i < tail + data + train + data + tail; i++)
1102 (*burst)[i] = psk8_table[7];
1103
1104 shape = shapeEdgeBurst(*burst);
1105 delete burst;
1106
1107 return shape;
1108}
1109
1110/*
1111 * Modulate 8-PSK burst. When empty pulse shaping (rotation only)
1112 * is enabled, the output vector length will be bit sequence length
1113 * times the SPS value. When pulse shaping is enabled, the output
1114 * vector length is fixed at 625 samples (156.25 sybols at 4 SPS).
1115 * Pulse shaped bit sequences that go beyond one burst are truncated.
1116 * Pulse shaping at anything but 4 SPS is not supported.
1117 */
1118signalVector *modulateEdgeBurst(const BitVector &bits,
1119 int sps, bool empty)
1120{
1121 signalVector *shape, *burst;
1122
1123 if ((sps != 4) && !empty)
1124 return NULL;
1125
1126 burst = mapEdgeSymbols(bits);
1127 if (!burst)
1128 return NULL;
1129
1130 if (empty)
1131 shape = rotateEdgeBurst(*burst, sps);
1132 else
1133 shape = shapeEdgeBurst(*burst);
1134
1135 delete burst;
1136 return shape;
1137}
1138
Thomas Tsoua57bc8a2013-09-05 08:16:47 +08001139static signalVector *modulateBurstBasic(const BitVector &bits,
1140 int guard_len, int sps)
1141{
1142 int burst_len;
Thomas Tsou6f4906e2013-11-09 02:40:18 -05001143 signalVector *pulse, *burst, *shaped;
Thomas Tsoua57bc8a2013-09-05 08:16:47 +08001144 signalVector::iterator burst_itr;
1145
Thomas Tsouc1f7c422013-10-11 13:49:55 -04001146 if (sps == 1)
1147 pulse = GSMPulse1->c0;
1148 else
Tom Tsou2079a3c2016-03-06 00:58:56 -08001149 pulse = GSMPulse4->c0;
Thomas Tsouc1f7c422013-10-11 13:49:55 -04001150
Thomas Tsoua57bc8a2013-09-05 08:16:47 +08001151 burst_len = sps * (bits.size() + guard_len);
1152
Thomas Tsou6f4906e2013-11-09 02:40:18 -05001153 burst = new signalVector(burst_len, pulse->size());
Thomas Tsou20eb6d62013-11-09 14:30:41 -05001154 burst->isReal(true);
Thomas Tsou6f4906e2013-11-09 02:40:18 -05001155 burst_itr = burst->begin();
Thomas Tsoua57bc8a2013-09-05 08:16:47 +08001156
1157 /* Raw bits are not differentially encoded */
1158 for (unsigned i = 0; i < bits.size(); i++) {
1159 *burst_itr = 2.0 * (bits[i] & 0x01) - 1.0;
1160 burst_itr += sps;
1161 }
1162
Thomas Tsou6f4906e2013-11-09 02:40:18 -05001163 GMSKRotate(*burst, sps);
Thomas Tsou20eb6d62013-11-09 14:30:41 -05001164 burst->isReal(false);
Thomas Tsoua57bc8a2013-09-05 08:16:47 +08001165
1166 /* Single Gaussian pulse approximation shaping */
Thomas Tsou6f4906e2013-11-09 02:40:18 -05001167 shaped = convolve(burst, pulse, NULL, START_ONLY);
1168
1169 delete burst;
Thomas Tsoua57bc8a2013-09-05 08:16:47 +08001170
1171 return shaped;
1172}
1173
Thomas Tsou3eaae802013-08-20 19:31:14 -04001174/* Assume input bits are not differentially encoded */
Thomas Tsou83e06892013-08-20 16:10:01 -04001175signalVector *modulateBurst(const BitVector &wBurst, int guardPeriodLength,
1176 int sps, bool emptyPulse)
dburgessb3a0ca42011-10-12 07:44:40 +00001177{
Thomas Tsou83e06892013-08-20 16:10:01 -04001178 if (emptyPulse)
Thomas Tsoua57bc8a2013-09-05 08:16:47 +08001179 return rotateBurst(wBurst, guardPeriodLength, sps);
1180 else if (sps == 4)
Tom Tsou4dfd64a2016-03-06 20:31:51 -08001181 return modulateBurstLaurent(wBurst);
Thomas Tsou83e06892013-08-20 16:10:01 -04001182 else
Thomas Tsoua57bc8a2013-09-05 08:16:47 +08001183 return modulateBurstBasic(wBurst, guardPeriodLength, sps);
dburgessb3a0ca42011-10-12 07:44:40 +00001184}
1185
Tom Tsou2079a3c2016-03-06 00:58:56 -08001186static void generateSincTable()
Thomas Tsou0e0e1f42013-11-09 22:08:51 -05001187{
1188 float x;
1189
1190 for (int i = 0; i < TABLESIZE; i++) {
1191 x = (float) i / TABLESIZE * 8 * M_PI;
1192 if (fabs(x) < 0.01) {
1193 sincTable[i] = 1.0f;
1194 continue;
1195 }
1196
1197 sincTable[i] = sinf(x) / x;
1198 }
1199}
1200
Thomas Tsoua57bc8a2013-09-05 08:16:47 +08001201float sinc(float x)
dburgessb3a0ca42011-10-12 07:44:40 +00001202{
Thomas Tsou0e0e1f42013-11-09 22:08:51 -05001203 if (fabs(x) >= 8 * M_PI)
1204 return 0.0;
1205
1206 int index = (int) floorf(fabs(x) / (8 * M_PI) * TABLESIZE);
1207
1208 return sincTable[index];
dburgessb3a0ca42011-10-12 07:44:40 +00001209}
1210
Thomas Tsouf79c4d02013-11-09 15:51:56 -06001211/*
1212 * Create fractional delay filterbank with Blackman-harris windowed
1213 * sinc function generator. The number of filters generated is specified
1214 * by the DELAYFILTS value.
1215 */
1216void generateDelayFilters()
1217{
1218 int h_len = 20;
1219 complex *data;
1220 signalVector *h;
1221 signalVector::iterator itr;
1222
1223 float k, sum;
1224 float a0 = 0.35875;
1225 float a1 = 0.48829;
1226 float a2 = 0.14128;
1227 float a3 = 0.01168;
1228
1229 for (int i = 0; i < DELAYFILTS; i++) {
1230 data = (complex *) convolve_h_alloc(h_len);
1231 h = new signalVector(data, 0, h_len);
1232 h->setAligned(true);
1233 h->isReal(true);
1234
1235 sum = 0.0;
1236 itr = h->end();
1237 for (int n = 0; n < h_len; n++) {
1238 k = (float) n;
1239 *--itr = (complex) sinc(M_PI_F *
1240 (k - (float) h_len / 2.0 - (float) i / DELAYFILTS));
1241 *itr *= a0 -
1242 a1 * cos(2 * M_PI * n / (h_len - 1)) +
1243 a2 * cos(4 * M_PI * n / (h_len - 1)) -
1244 a3 * cos(6 * M_PI * n / (h_len - 1));
1245
1246 sum += itr->real();
1247 }
1248
1249 itr = h->begin();
1250 for (int n = 0; n < h_len; n++)
1251 *itr++ /= sum;
1252
1253 delayFilters[i] = h;
1254 }
1255}
1256
Thomas Tsou94edaae2013-11-09 22:19:19 -05001257signalVector *delayVector(signalVector *in, signalVector *out, float delay)
dburgessb3a0ca42011-10-12 07:44:40 +00001258{
Thomas Tsouf79c4d02013-11-09 15:51:56 -06001259 int whole, index;
Thomas Tsou2c282f52013-10-08 21:34:35 -04001260 float frac;
Thomas Tsou94edaae2013-11-09 22:19:19 -05001261 signalVector *h, *shift, *fshift = NULL;
Thomas Tsou3eaae802013-08-20 19:31:14 -04001262
Thomas Tsou2c282f52013-10-08 21:34:35 -04001263 whole = floor(delay);
1264 frac = delay - whole;
1265
1266 /* Sinc interpolated fractional shift (if allowable) */
1267 if (fabs(frac) > 1e-2) {
Thomas Tsouf79c4d02013-11-09 15:51:56 -06001268 index = floorf(frac * (float) DELAYFILTS);
1269 h = delayFilters[index];
Thomas Tsou2c282f52013-10-08 21:34:35 -04001270
Thomas Tsou94edaae2013-11-09 22:19:19 -05001271 fshift = convolve(in, h, NULL, NO_DELAY);
1272 if (!fshift)
1273 return NULL;
dburgessb3a0ca42011-10-12 07:44:40 +00001274 }
1275
Thomas Tsou94edaae2013-11-09 22:19:19 -05001276 if (!fshift)
1277 shift = new signalVector(*in);
1278 else
1279 shift = fshift;
1280
Thomas Tsou2c282f52013-10-08 21:34:35 -04001281 /* Integer sample shift */
1282 if (whole < 0) {
1283 whole = -whole;
Thomas Tsou94edaae2013-11-09 22:19:19 -05001284 signalVector::iterator wBurstItr = shift->begin();
1285 signalVector::iterator shiftedItr = shift->begin() + whole;
Thomas Tsou2c282f52013-10-08 21:34:35 -04001286
Thomas Tsou94edaae2013-11-09 22:19:19 -05001287 while (shiftedItr < shift->end())
dburgessb3a0ca42011-10-12 07:44:40 +00001288 *wBurstItr++ = *shiftedItr++;
Thomas Tsou2c282f52013-10-08 21:34:35 -04001289
Thomas Tsou94edaae2013-11-09 22:19:19 -05001290 while (wBurstItr < shift->end())
1291 *wBurstItr++ = 0.0;
1292 } else if (whole >= 0) {
1293 signalVector::iterator wBurstItr = shift->end() - 1;
1294 signalVector::iterator shiftedItr = shift->end() - 1 - whole;
1295
1296 while (shiftedItr >= shift->begin())
dburgessb3a0ca42011-10-12 07:44:40 +00001297 *wBurstItr-- = *shiftedItr--;
Thomas Tsou94edaae2013-11-09 22:19:19 -05001298
1299 while (wBurstItr >= shift->begin())
dburgessb3a0ca42011-10-12 07:44:40 +00001300 *wBurstItr-- = 0.0;
1301 }
Thomas Tsou2c282f52013-10-08 21:34:35 -04001302
Thomas Tsou94edaae2013-11-09 22:19:19 -05001303 if (!out)
1304 return shift;
1305
1306 out->clone(*shift);
1307 delete shift;
1308 return out;
dburgessb3a0ca42011-10-12 07:44:40 +00001309}
Thomas Tsou2c282f52013-10-08 21:34:35 -04001310
dburgessb3a0ca42011-10-12 07:44:40 +00001311signalVector *gaussianNoise(int length,
1312 float variance,
1313 complex mean)
1314{
1315
1316 signalVector *noise = new signalVector(length);
1317 signalVector::iterator nPtr = noise->begin();
1318 float stddev = sqrtf(variance);
1319 while (nPtr < noise->end()) {
1320 float u1 = (float) rand()/ (float) RAND_MAX;
1321 while (u1==0.0)
1322 u1 = (float) rand()/ (float) RAND_MAX;
1323 float u2 = (float) rand()/ (float) RAND_MAX;
1324 float arg = 2.0*M_PI*u2;
1325 *nPtr = mean + stddev*complex(cos(arg),sin(arg))*sqrtf(-2.0*log(u1));
1326 nPtr++;
1327 }
1328
1329 return noise;
1330}
1331
1332complex interpolatePoint(const signalVector &inSig,
1333 float ix)
1334{
1335
1336 int start = (int) (floor(ix) - 10);
1337 if (start < 0) start = 0;
1338 int end = (int) (floor(ix) + 11);
1339 if ((unsigned) end > inSig.size()-1) end = inSig.size()-1;
1340
1341 complex pVal = 0.0;
Thomas Tsou20eb6d62013-11-09 14:30:41 -05001342 if (!inSig.isReal()) {
dburgessb3a0ca42011-10-12 07:44:40 +00001343 for (int i = start; i < end; i++)
1344 pVal += inSig[i] * sinc(M_PI_F*(i-ix));
1345 }
1346 else {
1347 for (int i = start; i < end; i++)
1348 pVal += inSig[i].real() * sinc(M_PI_F*(i-ix));
1349 }
1350
1351 return pVal;
1352}
1353
Thomas Tsou8181b012013-08-20 21:17:19 -04001354static complex fastPeakDetect(const signalVector &rxBurst, float *index)
1355{
1356 float val, max = 0.0f;
1357 complex amp;
1358 int _index = -1;
1359
Thomas Tsou3f32ab52013-11-15 16:32:54 -05001360 for (size_t i = 0; i < rxBurst.size(); i++) {
Thomas Tsou8181b012013-08-20 21:17:19 -04001361 val = rxBurst[i].norm2();
1362 if (val > max) {
1363 max = val;
1364 _index = i;
1365 amp = rxBurst[i];
1366 }
1367 }
1368
1369 if (index)
1370 *index = (float) _index;
1371
1372 return amp;
1373}
1374
dburgessb3a0ca42011-10-12 07:44:40 +00001375complex peakDetect(const signalVector &rxBurst,
1376 float *peakIndex,
1377 float *avgPwr)
1378{
1379
1380
1381 complex maxVal = 0.0;
1382 float maxIndex = -1;
1383 float sumPower = 0.0;
1384
1385 for (unsigned int i = 0; i < rxBurst.size(); i++) {
1386 float samplePower = rxBurst[i].norm2();
1387 if (samplePower > maxVal.real()) {
1388 maxVal = samplePower;
1389 maxIndex = i;
1390 }
1391 sumPower += samplePower;
1392 }
1393
1394 // interpolate around the peak
1395 // to save computation, we'll use early-late balancing
1396 float earlyIndex = maxIndex-1;
1397 float lateIndex = maxIndex+1;
1398
1399 float incr = 0.5;
1400 while (incr > 1.0/1024.0) {
1401 complex earlyP = interpolatePoint(rxBurst,earlyIndex);
1402 complex lateP = interpolatePoint(rxBurst,lateIndex);
1403 if (earlyP < lateP)
1404 earlyIndex += incr;
1405 else if (earlyP > lateP)
1406 earlyIndex -= incr;
1407 else break;
1408 incr /= 2.0;
1409 lateIndex = earlyIndex + 2.0;
1410 }
1411
1412 maxIndex = earlyIndex + 1.0;
1413 maxVal = interpolatePoint(rxBurst,maxIndex);
1414
1415 if (peakIndex!=NULL)
1416 *peakIndex = maxIndex;
1417
1418 if (avgPwr!=NULL)
1419 *avgPwr = (sumPower-maxVal.norm2()) / (rxBurst.size()-1);
1420
1421 return maxVal;
1422
1423}
1424
1425void scaleVector(signalVector &x,
1426 complex scale)
1427{
Thomas Tsou7e4e5362013-10-30 21:18:55 -04001428#ifdef HAVE_NEON
1429 int len = x.size();
1430
1431 scale_complex((float *) x.begin(),
1432 (float *) x.begin(),
1433 (float *) &scale, len);
1434#else
dburgessb3a0ca42011-10-12 07:44:40 +00001435 signalVector::iterator xP = x.begin();
1436 signalVector::iterator xPEnd = x.end();
Thomas Tsou20eb6d62013-11-09 14:30:41 -05001437 if (!x.isReal()) {
dburgessb3a0ca42011-10-12 07:44:40 +00001438 while (xP < xPEnd) {
1439 *xP = *xP * scale;
1440 xP++;
1441 }
1442 }
1443 else {
1444 while (xP < xPEnd) {
1445 *xP = xP->real() * scale;
1446 xP++;
1447 }
1448 }
Thomas Tsou7e4e5362013-10-30 21:18:55 -04001449#endif
dburgessb3a0ca42011-10-12 07:44:40 +00001450}
1451
1452/** in-place conjugation */
1453void conjugateVector(signalVector &x)
1454{
Thomas Tsou20eb6d62013-11-09 14:30:41 -05001455 if (x.isReal()) return;
dburgessb3a0ca42011-10-12 07:44:40 +00001456 signalVector::iterator xP = x.begin();
1457 signalVector::iterator xPEnd = x.end();
1458 while (xP < xPEnd) {
1459 *xP = xP->conj();
1460 xP++;
1461 }
1462}
1463
1464
1465// in-place addition!!
1466bool addVector(signalVector &x,
1467 signalVector &y)
1468{
1469 signalVector::iterator xP = x.begin();
1470 signalVector::iterator yP = y.begin();
1471 signalVector::iterator xPEnd = x.end();
1472 signalVector::iterator yPEnd = y.end();
1473 while ((xP < xPEnd) && (yP < yPEnd)) {
1474 *xP = *xP + *yP;
1475 xP++; yP++;
1476 }
1477 return true;
1478}
1479
1480// in-place multiplication!!
1481bool multVector(signalVector &x,
1482 signalVector &y)
1483{
1484 signalVector::iterator xP = x.begin();
1485 signalVector::iterator yP = y.begin();
1486 signalVector::iterator xPEnd = x.end();
1487 signalVector::iterator yPEnd = y.end();
1488 while ((xP < xPEnd) && (yP < yPEnd)) {
1489 *xP = (*xP) * (*yP);
1490 xP++; yP++;
1491 }
1492 return true;
1493}
1494
Tom Tsou2079a3c2016-03-06 00:58:56 -08001495static bool generateMidamble(int sps, int tsc)
dburgessb3a0ca42011-10-12 07:44:40 +00001496{
Thomas Tsoue5dcfc42013-08-20 16:27:12 -04001497 bool status = true;
Thomas Tsouc1f7c422013-10-11 13:49:55 -04001498 float toa;
Thomas Tsoue5dcfc42013-08-20 16:27:12 -04001499 complex *data = NULL;
1500 signalVector *autocorr = NULL, *midamble = NULL;
Thomas Tsou3eaae802013-08-20 19:31:14 -04001501 signalVector *midMidamble = NULL, *_midMidamble = NULL;
Thomas Tsoue5dcfc42013-08-20 16:27:12 -04001502
Thomas Tsou3eaae802013-08-20 19:31:14 -04001503 if ((tsc < 0) || (tsc > 7))
dburgessb3a0ca42011-10-12 07:44:40 +00001504 return false;
1505
Thomas Tsoue5dcfc42013-08-20 16:27:12 -04001506 delete gMidambles[tsc];
Thomas Tsou3eaae802013-08-20 19:31:14 -04001507
Thomas Tsoue5dcfc42013-08-20 16:27:12 -04001508 /* Use middle 16 bits of each TSC. Correlation sequence is not pulse shaped */
1509 midMidamble = modulateBurst(gTrainingSequence[tsc].segment(5,16), 0, sps, true);
1510 if (!midMidamble)
1511 return false;
1512
Thomas Tsou3eaae802013-08-20 19:31:14 -04001513 /* Simulated receive sequence is pulse shaped */
Thomas Tsoue5dcfc42013-08-20 16:27:12 -04001514 midamble = modulateBurst(gTrainingSequence[tsc], 0, sps, false);
1515 if (!midamble) {
1516 status = false;
1517 goto release;
dburgessb3a0ca42011-10-12 07:44:40 +00001518 }
Thomas Tsou3eaae802013-08-20 19:31:14 -04001519
dburgessb3a0ca42011-10-12 07:44:40 +00001520 // NOTE: Because ideal TSC 16-bit midamble is 66 symbols into burst,
1521 // the ideal TSC has an + 180 degree phase shift,
1522 // due to the pi/2 frequency shift, that
1523 // needs to be accounted for.
1524 // 26-midamble is 61 symbols into burst, has +90 degree phase shift.
Thomas Tsoue5dcfc42013-08-20 16:27:12 -04001525 scaleVector(*midMidamble, complex(-1.0, 0.0));
1526 scaleVector(*midamble, complex(0.0, 1.0));
dburgessb3a0ca42011-10-12 07:44:40 +00001527
Thomas Tsoue5dcfc42013-08-20 16:27:12 -04001528 conjugateVector(*midMidamble);
dburgessb3a0ca42011-10-12 07:44:40 +00001529
Thomas Tsou3eaae802013-08-20 19:31:14 -04001530 /* For SSE alignment, reallocate the midamble sequence on 16-byte boundary */
1531 data = (complex *) convolve_h_alloc(midMidamble->size());
1532 _midMidamble = new signalVector(data, 0, midMidamble->size());
1533 _midMidamble->setAligned(true);
1534 memcpy(_midMidamble->begin(), midMidamble->begin(),
1535 midMidamble->size() * sizeof(complex));
1536
1537 autocorr = convolve(midamble, _midMidamble, NULL, NO_DELAY);
Thomas Tsoue5dcfc42013-08-20 16:27:12 -04001538 if (!autocorr) {
1539 status = false;
1540 goto release;
1541 }
dburgessb3a0ca42011-10-12 07:44:40 +00001542
Thomas Tsoue5dcfc42013-08-20 16:27:12 -04001543 gMidambles[tsc] = new CorrelationSequence;
Thomas Tsou3eaae802013-08-20 19:31:14 -04001544 gMidambles[tsc]->buffer = data;
1545 gMidambles[tsc]->sequence = _midMidamble;
Thomas Tsouc1f7c422013-10-11 13:49:55 -04001546 gMidambles[tsc]->gain = peakDetect(*autocorr, &toa, NULL);
1547
1548 /* For 1 sps only
1549 * (Half of correlation length - 1) + midpoint of pulse shape + remainder
1550 * 13.5 = (16 / 2 - 1) + 1.5 + (26 - 10) / 2
1551 */
1552 if (sps == 1)
1553 gMidambles[tsc]->toa = toa - 13.5;
1554 else
1555 gMidambles[tsc]->toa = 0;
dburgessb3a0ca42011-10-12 07:44:40 +00001556
Thomas Tsoue5dcfc42013-08-20 16:27:12 -04001557release:
dburgessb3a0ca42011-10-12 07:44:40 +00001558 delete autocorr;
1559 delete midamble;
Thomas Tsou3eaae802013-08-20 19:31:14 -04001560 delete midMidamble;
dburgessb3a0ca42011-10-12 07:44:40 +00001561
Thomas Tsoue5dcfc42013-08-20 16:27:12 -04001562 if (!status) {
Thomas Tsou3eaae802013-08-20 19:31:14 -04001563 delete _midMidamble;
1564 free(data);
Thomas Tsoue5dcfc42013-08-20 16:27:12 -04001565 gMidambles[tsc] = NULL;
1566 }
1567
1568 return status;
dburgessb3a0ca42011-10-12 07:44:40 +00001569}
1570
Tom Tsoud3253432016-03-06 03:08:01 -08001571CorrelationSequence *generateEdgeMidamble(int tsc)
1572{
1573 complex *data = NULL;
1574 signalVector *midamble = NULL, *_midamble = NULL;
1575 CorrelationSequence *seq;
1576
1577 if ((tsc < 0) || (tsc > 7))
1578 return NULL;
1579
1580 /* Use middle 48 bits of each TSC. Correlation sequence is not pulse shaped */
1581 const BitVector *bits = &gEdgeTrainingSequence[tsc];
1582 midamble = modulateEdgeBurst(bits->segment(15, 48), 1, true);
1583 if (!midamble)
1584 return NULL;
1585
1586 conjugateVector(*midamble);
1587
1588 data = (complex *) convolve_h_alloc(midamble->size());
1589 _midamble = new signalVector(data, 0, midamble->size());
1590 _midamble->setAligned(true);
1591 memcpy(_midamble->begin(), midamble->begin(),
1592 midamble->size() * sizeof(complex));
1593
1594 /* Channel gain is an empirically measured value */
1595 seq = new CorrelationSequence;
1596 seq->buffer = data;
1597 seq->sequence = _midamble;
1598 seq->gain = Complex<float>(-19.6432, 19.5006) / 1.18;
1599 seq->toa = 0;
1600
1601 delete midamble;
1602
1603 return seq;
1604}
1605
Tom Tsou2079a3c2016-03-06 00:58:56 -08001606static bool generateRACHSequence(int sps)
dburgessb3a0ca42011-10-12 07:44:40 +00001607{
Thomas Tsoue5dcfc42013-08-20 16:27:12 -04001608 bool status = true;
Thomas Tsouc1f7c422013-10-11 13:49:55 -04001609 float toa;
Thomas Tsoue5dcfc42013-08-20 16:27:12 -04001610 complex *data = NULL;
1611 signalVector *autocorr = NULL;
Thomas Tsou3eaae802013-08-20 19:31:14 -04001612 signalVector *seq0 = NULL, *seq1 = NULL, *_seq1 = NULL;
Thomas Tsoue5dcfc42013-08-20 16:27:12 -04001613
1614 delete gRACHSequence;
1615
1616 seq0 = modulateBurst(gRACHSynchSequence, 0, sps, false);
1617 if (!seq0)
1618 return false;
1619
1620 seq1 = modulateBurst(gRACHSynchSequence.segment(0, 40), 0, sps, true);
1621 if (!seq1) {
1622 status = false;
1623 goto release;
dburgessb3a0ca42011-10-12 07:44:40 +00001624 }
1625
Thomas Tsoue5dcfc42013-08-20 16:27:12 -04001626 conjugateVector(*seq1);
dburgessb3a0ca42011-10-12 07:44:40 +00001627
Thomas Tsou3eaae802013-08-20 19:31:14 -04001628 /* For SSE alignment, reallocate the midamble sequence on 16-byte boundary */
1629 data = (complex *) convolve_h_alloc(seq1->size());
1630 _seq1 = new signalVector(data, 0, seq1->size());
1631 _seq1->setAligned(true);
1632 memcpy(_seq1->begin(), seq1->begin(), seq1->size() * sizeof(complex));
1633
1634 autocorr = convolve(seq0, _seq1, autocorr, NO_DELAY);
1635 if (!autocorr) {
Thomas Tsoue5dcfc42013-08-20 16:27:12 -04001636 status = false;
1637 goto release;
1638 }
dburgessb3a0ca42011-10-12 07:44:40 +00001639
1640 gRACHSequence = new CorrelationSequence;
Thomas Tsou3eaae802013-08-20 19:31:14 -04001641 gRACHSequence->sequence = _seq1;
1642 gRACHSequence->buffer = data;
Thomas Tsouc1f7c422013-10-11 13:49:55 -04001643 gRACHSequence->gain = peakDetect(*autocorr, &toa, NULL);
1644
1645 /* For 1 sps only
1646 * (Half of correlation length - 1) + midpoint of pulse shaping filer
1647 * 20.5 = (40 / 2 - 1) + 1.5
1648 */
1649 if (sps == 1)
1650 gRACHSequence->toa = toa - 20.5;
1651 else
1652 gRACHSequence->toa = 0.0;
Thomas Tsoue5dcfc42013-08-20 16:27:12 -04001653
1654release:
dburgessb3a0ca42011-10-12 07:44:40 +00001655 delete autocorr;
Thomas Tsoue5dcfc42013-08-20 16:27:12 -04001656 delete seq0;
Thomas Tsou3eaae802013-08-20 19:31:14 -04001657 delete seq1;
dburgessb3a0ca42011-10-12 07:44:40 +00001658
Thomas Tsoue5dcfc42013-08-20 16:27:12 -04001659 if (!status) {
Thomas Tsou3eaae802013-08-20 19:31:14 -04001660 delete _seq1;
1661 free(data);
Thomas Tsoue5dcfc42013-08-20 16:27:12 -04001662 gRACHSequence = NULL;
1663 }
dburgessb3a0ca42011-10-12 07:44:40 +00001664
Thomas Tsoue5dcfc42013-08-20 16:27:12 -04001665 return status;
dburgessb3a0ca42011-10-12 07:44:40 +00001666}
Thomas Tsou3eaae802013-08-20 19:31:14 -04001667
Thomas Tsou865bca42013-08-21 20:58:00 -04001668static float computePeakRatio(signalVector *corr,
1669 int sps, float toa, complex amp)
dburgessb3a0ca42011-10-12 07:44:40 +00001670{
Thomas Tsou865bca42013-08-21 20:58:00 -04001671 int num = 0;
1672 complex *peak;
1673 float rms, avg = 0.0;
dburgessb3a0ca42011-10-12 07:44:40 +00001674
Thomas Tsou865bca42013-08-21 20:58:00 -04001675 /* Check for bogus results */
1676 if ((toa < 0.0) || (toa > corr->size()))
1677 return 0.0;
dburgessb3a0ca42011-10-12 07:44:40 +00001678
Alexander Chemeris1e9b4d52015-06-04 19:05:28 -04001679 peak = corr->begin() + (int) rint(toa);
1680
Thomas Tsou3eaae802013-08-20 19:31:14 -04001681 for (int i = 2 * sps; i <= 5 * sps; i++) {
Thomas Tsou865bca42013-08-21 20:58:00 -04001682 if (peak - i >= corr->begin()) {
Thomas Tsou3eaae802013-08-20 19:31:14 -04001683 avg += (peak - i)->norm2();
1684 num++;
1685 }
Thomas Tsou865bca42013-08-21 20:58:00 -04001686 if (peak + i < corr->end()) {
Thomas Tsou3eaae802013-08-20 19:31:14 -04001687 avg += (peak + i)->norm2();
1688 num++;
1689 }
dburgessb3a0ca42011-10-12 07:44:40 +00001690 }
1691
Thomas Tsou3eaae802013-08-20 19:31:14 -04001692 if (num < 2)
Thomas Tsou865bca42013-08-21 20:58:00 -04001693 return 0.0;
dburgessb3a0ca42011-10-12 07:44:40 +00001694
Thomas Tsou3eaae802013-08-20 19:31:14 -04001695 rms = sqrtf(avg / (float) num) + 0.00001;
dburgessb3a0ca42011-10-12 07:44:40 +00001696
Thomas Tsou865bca42013-08-21 20:58:00 -04001697 return (amp.abs()) / rms;
dburgessb3a0ca42011-10-12 07:44:40 +00001698}
1699
1700bool energyDetect(signalVector &rxBurst,
1701 unsigned windowLength,
1702 float detectThreshold,
1703 float *avgPwr)
1704{
1705
1706 signalVector::const_iterator windowItr = rxBurst.begin(); //+rxBurst.size()/2 - 5*windowLength/2;
1707 float energy = 0.0;
1708 if (windowLength < 0) windowLength = 20;
1709 if (windowLength > rxBurst.size()) windowLength = rxBurst.size();
1710 for (unsigned i = 0; i < windowLength; i++) {
1711 energy += windowItr->norm2();
1712 windowItr+=4;
1713 }
1714 if (avgPwr) *avgPwr = energy/windowLength;
dburgessb3a0ca42011-10-12 07:44:40 +00001715 return (energy/windowLength > detectThreshold*detectThreshold);
1716}
dburgessb3a0ca42011-10-12 07:44:40 +00001717
Thomas Tsou865bca42013-08-21 20:58:00 -04001718/*
1719 * Detect a burst based on correlation and peak-to-average ratio
1720 *
1721 * For one sampler-per-symbol, perform fast peak detection (no interpolation)
1722 * for initial gating. We do this because energy detection should be disabled.
1723 * For higher oversampling values, we assume the energy detector is in place
1724 * and we run full interpolating peak detection.
1725 */
1726static int detectBurst(signalVector &burst,
1727 signalVector &corr, CorrelationSequence *sync,
1728 float thresh, int sps, complex *amp, float *toa,
1729 int start, int len)
dburgessb3a0ca42011-10-12 07:44:40 +00001730{
Tom Tsoud3253432016-03-06 03:08:01 -08001731 signalVector *corr_in, *dec = NULL;
1732
1733 if (sps == 4) {
1734 dec = downsampleBurst(burst);
1735 corr_in = dec;
1736 sps = 1;
1737 } else {
1738 corr_in = &burst;
1739 }
1740
Thomas Tsou865bca42013-08-21 20:58:00 -04001741 /* Correlate */
Tom Tsoud3253432016-03-06 03:08:01 -08001742 if (!convolve(corr_in, sync->sequence, &corr,
1743 CUSTOM, start, len, 1, 0)) {
1744 delete dec;
Thomas Tsou3eaae802013-08-20 19:31:14 -04001745 return -1;
dburgessb3a0ca42011-10-12 07:44:40 +00001746 }
1747
Tom Tsoud3253432016-03-06 03:08:01 -08001748 delete dec;
1749
1750 /* Running at the downsampled rate at this point */
1751 sps = 1;
1752
Thomas Tsou865bca42013-08-21 20:58:00 -04001753 /* Peak detection - place restrictions at correlation edges */
1754 *amp = fastPeakDetect(corr, toa);
Thomas Tsou3eaae802013-08-20 19:31:14 -04001755
Thomas Tsou865bca42013-08-21 20:58:00 -04001756 if ((*toa < 3 * sps) || (*toa > len - 3 * sps))
1757 return 0;
Thomas Tsou3eaae802013-08-20 19:31:14 -04001758
Thomas Tsou865bca42013-08-21 20:58:00 -04001759 /* Peak -to-average ratio */
1760 if (computePeakRatio(&corr, sps, *toa, *amp) < thresh)
1761 return 0;
1762
1763 /* Compute peak-to-average ratio. Reject if we don't have enough values */
1764 *amp = peakDetect(corr, toa, NULL);
1765
1766 /* Normalize our channel gain */
1767 *amp = *amp / sync->gain;
1768
Thomas Tsoua57bc8a2013-09-05 08:16:47 +08001769 /* Compenate for residual rotation with dual Laurent pulse */
1770 if (sps == 4)
1771 *amp = *amp * complex(0.0, 1.0);
1772
Thomas Tsouc1f7c422013-10-11 13:49:55 -04001773 /* Compensate for residuate time lag */
1774 *toa = *toa - sync->toa;
1775
Thomas Tsou865bca42013-08-21 20:58:00 -04001776 return 1;
1777}
1778
Alexander Chemeris954b1182015-06-04 15:39:41 -04001779static float maxAmplitude(signalVector &burst)
Tom Tsou577cd022015-05-18 13:57:54 -07001780{
Alexander Chemeris954b1182015-06-04 15:39:41 -04001781 float max = 0.0;
1782 for (size_t i = 0; i < burst.size(); i++) {
1783 if (fabs(burst[i].real()) > max)
1784 max = fabs(burst[i].real());
1785 if (fabs(burst[i].imag()) > max)
1786 max = fabs(burst[i].imag());
1787 }
Tom Tsou577cd022015-05-18 13:57:54 -07001788
Alexander Chemeris954b1182015-06-04 15:39:41 -04001789 return max;
Tom Tsou577cd022015-05-18 13:57:54 -07001790}
1791
Alexander Chemeris130a8002015-06-09 20:52:11 -04001792/*
1793 * RACH/Normal burst detection with clipping detection
Thomas Tsou865bca42013-08-21 20:58:00 -04001794 *
1795 * Correlation window parameters:
Alexander Chemeris130a8002015-06-09 20:52:11 -04001796 * target: Tail bits + burst length
1797 * head: Search symbols before target
1798 * tail: Search symbols after target
Thomas Tsou865bca42013-08-21 20:58:00 -04001799 */
Alexander Chemeris130a8002015-06-09 20:52:11 -04001800int detectGeneralBurst(signalVector &rxBurst,
1801 float thresh,
1802 int sps,
1803 complex &amp,
1804 float &toa,
1805 int target, int head, int tail,
1806 CorrelationSequence *sync)
Thomas Tsou865bca42013-08-21 20:58:00 -04001807{
Alexander Chemeris130a8002015-06-09 20:52:11 -04001808 int rc, start, len;
Alexander Chemeris954b1182015-06-04 15:39:41 -04001809 bool clipping = false;
Thomas Tsoub075dd22013-11-09 22:25:46 -05001810 signalVector *corr;
Thomas Tsou865bca42013-08-21 20:58:00 -04001811
1812 if ((sps != 1) && (sps != 4))
Tom Tsou577cd022015-05-18 13:57:54 -07001813 return -SIGERR_UNSUPPORTED;
1814
Alexander Chemeris954b1182015-06-04 15:39:41 -04001815 // Detect potential clipping
1816 // We still may be able to demod the burst, so we'll give it a try
1817 // and only report clipping if we can't demod.
1818 float maxAmpl = maxAmplitude(rxBurst);
1819 if (maxAmpl > CLIP_THRESH) {
1820 LOG(DEBUG) << "max burst amplitude: " << maxAmpl << " is above the clipping threshold: " << CLIP_THRESH << std::endl;
1821 clipping = true;
1822 }
Thomas Tsou865bca42013-08-21 20:58:00 -04001823
Tom Tsoud3253432016-03-06 03:08:01 -08001824 start = target - head - 1;
1825 len = head + tail;
Thomas Tsoub075dd22013-11-09 22:25:46 -05001826 corr = new signalVector(len);
Thomas Tsou865bca42013-08-21 20:58:00 -04001827
Thomas Tsoub075dd22013-11-09 22:25:46 -05001828 rc = detectBurst(rxBurst, *corr, sync,
Alexander Chemeris130a8002015-06-09 20:52:11 -04001829 thresh, sps, &amp, &toa, start, len);
Thomas Tsoub075dd22013-11-09 22:25:46 -05001830 delete corr;
1831
Thomas Tsou865bca42013-08-21 20:58:00 -04001832 if (rc < 0) {
Tom Tsou577cd022015-05-18 13:57:54 -07001833 return -SIGERR_INTERNAL;
Thomas Tsou865bca42013-08-21 20:58:00 -04001834 } else if (!rc) {
Alexander Chemeris130a8002015-06-09 20:52:11 -04001835 amp = 0.0f;
1836 toa = 0.0f;
Alexander Chemeris954b1182015-06-04 15:39:41 -04001837 return clipping?-SIGERR_CLIP:SIGERR_NONE;
dburgessb3a0ca42011-10-12 07:44:40 +00001838 }
1839
Thomas Tsou865bca42013-08-21 20:58:00 -04001840 /* Subtract forward search bits from delay */
Tom Tsoud3253432016-03-06 03:08:01 -08001841 toa -= head;
Thomas Tsou3eaae802013-08-20 19:31:14 -04001842
Thomas Tsou865bca42013-08-21 20:58:00 -04001843 return 1;
1844}
Thomas Tsou3eaae802013-08-20 19:31:14 -04001845
Alexander Chemeris130a8002015-06-09 20:52:11 -04001846
1847/*
1848 * RACH burst detection
1849 *
1850 * Correlation window parameters:
1851 * target: Tail bits + RACH length (reduced from 41 to a multiple of 4)
Alexander Chemeris78d1fc92016-03-19 21:16:22 +03001852 * head: Search 4 symbols before target
1853 * tail: Search 4 symbols + maximum expected delay
Alexander Chemeris130a8002015-06-09 20:52:11 -04001854 */
1855int detectRACHBurst(signalVector &rxBurst,
1856 float thresh,
1857 int sps,
1858 complex &amp,
Alexander Chemeris78d1fc92016-03-19 21:16:22 +03001859 float &toa,
1860 unsigned maxTOA)
Alexander Chemeris130a8002015-06-09 20:52:11 -04001861{
1862 int rc, target, head, tail;
1863 CorrelationSequence *sync;
1864
1865 target = 8 + 40;
1866 head = 4;
Alexander Chemeris78d1fc92016-03-19 21:16:22 +03001867 tail = 4 + maxTOA;
Alexander Chemeris130a8002015-06-09 20:52:11 -04001868 sync = gRACHSequence;
1869
1870 rc = detectGeneralBurst(rxBurst, thresh, sps, amp, toa,
1871 target, head, tail, sync);
1872
1873 return rc;
1874}
1875
Thomas Tsou865bca42013-08-21 20:58:00 -04001876/*
1877 * Normal burst detection
1878 *
1879 * Correlation window parameters:
1880 * target: Tail + data + mid-midamble + 1/2 remaining midamblebits
Thomas Tsoudafb3372013-09-18 16:21:26 -04001881 * head: Search 4 symbols before target
1882 * tail: Search 4 symbols + maximum expected delay
Thomas Tsou865bca42013-08-21 20:58:00 -04001883 */
1884int analyzeTrafficBurst(signalVector &rxBurst, unsigned tsc, float thresh,
Tom Tsou46569402016-03-06 01:59:38 -08001885 int sps, complex &amp, float &toa, unsigned max_toa)
Thomas Tsou865bca42013-08-21 20:58:00 -04001886{
Alexander Chemeris130a8002015-06-09 20:52:11 -04001887 int rc, target, head, tail;
Thomas Tsou865bca42013-08-21 20:58:00 -04001888 CorrelationSequence *sync;
1889
Alexander Chemeris130a8002015-06-09 20:52:11 -04001890 if ((tsc < 0) || (tsc > 7))
Tom Tsou577cd022015-05-18 13:57:54 -07001891 return -SIGERR_UNSUPPORTED;
1892
Thomas Tsou865bca42013-08-21 20:58:00 -04001893 target = 3 + 58 + 16 + 5;
Thomas Tsoudafb3372013-09-18 16:21:26 -04001894 head = 4;
1895 tail = 4 + max_toa;
Thomas Tsou865bca42013-08-21 20:58:00 -04001896 sync = gMidambles[tsc];
Thomas Tsou865bca42013-08-21 20:58:00 -04001897
Alexander Chemeris130a8002015-06-09 20:52:11 -04001898 rc = detectGeneralBurst(rxBurst, thresh, sps, amp, toa,
1899 target, head, tail, sync);
Alexander Chemeris130a8002015-06-09 20:52:11 -04001900 return rc;
dburgessb3a0ca42011-10-12 07:44:40 +00001901}
1902
Tom Tsoud3253432016-03-06 03:08:01 -08001903int detectEdgeBurst(signalVector &rxBurst, unsigned tsc, float thresh,
1904 int sps, complex &amp, float &toa, unsigned max_toa)
1905{
1906 int rc, target, head, tail;
1907 CorrelationSequence *sync;
1908
1909 if ((tsc < 0) || (tsc > 7))
1910 return -SIGERR_UNSUPPORTED;
1911
1912 target = 3 + 58 + 16 + 5;
1913 head = 5;
1914 tail = 5 + max_toa;
1915 sync = gEdgeMidambles[tsc];
1916
1917 rc = detectGeneralBurst(rxBurst, thresh, sps, amp, toa,
1918 target, head, tail, sync);
1919 return rc;
1920}
1921
1922signalVector *downsampleBurst(signalVector &burst)
1923{
1924 size_t ilen = DOWNSAMPLE_IN_LEN, olen = DOWNSAMPLE_OUT_LEN;
1925
1926 signalVector *out = new signalVector(olen);
1927 memcpy(dnsampler_in->begin(), burst.begin(), ilen * 2 * sizeof(float));
1928
1929 dnsampler->rotate((float *) dnsampler_in->begin(), ilen,
1930 (float *) out->begin(), olen);
1931 return out;
1932};
1933
Thomas Tsou94edaae2013-11-09 22:19:19 -05001934signalVector *decimateVector(signalVector &wVector, size_t factor)
dburgessb3a0ca42011-10-12 07:44:40 +00001935{
Thomas Tsou94edaae2013-11-09 22:19:19 -05001936 signalVector *dec;
dburgessb3a0ca42011-10-12 07:44:40 +00001937
Thomas Tsou94edaae2013-11-09 22:19:19 -05001938 if (factor <= 1)
1939 return NULL;
dburgessb3a0ca42011-10-12 07:44:40 +00001940
Thomas Tsou94edaae2013-11-09 22:19:19 -05001941 dec = new signalVector(wVector.size() / factor);
1942 dec->isReal(wVector.isReal());
dburgessb3a0ca42011-10-12 07:44:40 +00001943
Thomas Tsou94edaae2013-11-09 22:19:19 -05001944 signalVector::iterator itr = dec->begin();
1945 for (size_t i = 0; i < wVector.size(); i += factor)
1946 *itr++ = wVector[i];
1947
1948 return dec;
dburgessb3a0ca42011-10-12 07:44:40 +00001949}
1950
Tom Tsoud3253432016-03-06 03:08:01 -08001951/*
1952 * Soft 8-PSK decoding using Manhattan distance metric
1953 */
1954static SoftVector *softSliceEdgeBurst(signalVector &burst)
1955{
1956 size_t nsyms = 148;
1957
1958 if (burst.size() < nsyms)
1959 return NULL;
1960
1961 signalVector::iterator itr;
1962 SoftVector *bits = new SoftVector(nsyms * 3);
1963
1964 /*
1965 * Bits 0 and 1 - First and second bits of the symbol respectively
1966 */
1967 rotateBurst2(burst, -M_PI / 8.0);
1968 itr = burst.begin();
1969 for (size_t i = 0; i < nsyms; i++) {
1970 (*bits)[3 * i + 0] = -itr->imag();
1971 (*bits)[3 * i + 1] = itr->real();
1972 itr++;
1973 }
1974
1975 /*
1976 * Bit 2 - Collapse symbols into quadrant 0 (positive X and Y).
1977 * Decision area is then simplified to X=Y axis. Rotate again to
1978 * place decision boundary on X-axis.
1979 */
1980 itr = burst.begin();
1981 for (size_t i = 0; i < burst.size(); i++) {
1982 burst[i] = Complex<float>(fabs(itr->real()), fabs(itr->imag()));
1983 itr++;
1984 }
1985
1986 rotateBurst2(burst, -M_PI / 4.0);
1987 itr = burst.begin();
1988 for (size_t i = 0; i < nsyms; i++) {
1989 (*bits)[3 * i + 2] = -itr->imag();
1990 itr++;
1991 }
1992
1993 signalVector soft(bits->size());
1994 for (size_t i = 0; i < bits->size(); i++)
1995 soft[i] = (*bits)[i];
1996
1997 return bits;
1998}
1999
2000/*
Tom Tsou7fec3032016-03-06 22:33:20 -08002001 * Shared portion of GMSK and EDGE demodulators consisting of timing
2002 * recovery and single tap channel correction. For 4 SPS (if activated),
2003 * the output is downsampled prior to the 1 SPS modulation specific
2004 * stages.
2005 */
2006static signalVector *demodCommon(signalVector &burst, int sps,
2007 complex chan, float toa)
2008{
2009 signalVector *delay, *dec;
2010
2011 if ((sps != 1) && (sps != 4))
2012 return NULL;
2013
2014 scaleVector(burst, (complex) 1.0 / chan);
2015 delay = delayVector(&burst, NULL, -toa * (float) sps);
2016
2017 if (sps == 1)
2018 return delay;
2019
2020 dec = downsampleBurst(*delay);
2021
2022 delete delay;
2023 return dec;
2024}
2025
2026/*
Tom Tsoud3253432016-03-06 03:08:01 -08002027 * Demodulate GSMK burst. Prior to symbol rotation, operate at
2028 * 4 SPS (if activated) to minimize distortion through the fractional
2029 * delay filters. Symbol rotation and after always operates at 1 SPS.
2030 */
Thomas Tsou83e06892013-08-20 16:10:01 -04002031SoftVector *demodulateBurst(signalVector &rxBurst, int sps,
Thomas Tsou94edaae2013-11-09 22:19:19 -05002032 complex channel, float TOA)
dburgessb3a0ca42011-10-12 07:44:40 +00002033{
Thomas Tsou94edaae2013-11-09 22:19:19 -05002034 SoftVector *bits;
Tom Tsou7fec3032016-03-06 22:33:20 -08002035 signalVector *dec;
dburgessb3a0ca42011-10-12 07:44:40 +00002036
Tom Tsou7fec3032016-03-06 22:33:20 -08002037 dec = demodCommon(rxBurst, sps, channel, TOA);
2038 if (!dec)
2039 return NULL;
dburgessb3a0ca42011-10-12 07:44:40 +00002040
Tom Tsoud3253432016-03-06 03:08:01 -08002041 /* Shift up by a quarter of a frequency */
2042 GMSKReverseRotate(*dec, 1);
Thomas Tsou94edaae2013-11-09 22:19:19 -05002043 vectorSlicer(dec);
dburgessb3a0ca42011-10-12 07:44:40 +00002044
Thomas Tsou94edaae2013-11-09 22:19:19 -05002045 bits = new SoftVector(dec->size());
dburgessb3a0ca42011-10-12 07:44:40 +00002046
Thomas Tsou94edaae2013-11-09 22:19:19 -05002047 SoftVector::iterator bit_itr = bits->begin();
2048 signalVector::iterator burst_itr = dec->begin();
dburgessb3a0ca42011-10-12 07:44:40 +00002049
Thomas Tsou94edaae2013-11-09 22:19:19 -05002050 for (; burst_itr < dec->end(); burst_itr++)
2051 *bit_itr++ = burst_itr->real();
dburgessb3a0ca42011-10-12 07:44:40 +00002052
Thomas Tsou94edaae2013-11-09 22:19:19 -05002053 delete dec;
dburgessb3a0ca42011-10-12 07:44:40 +00002054
Thomas Tsou94edaae2013-11-09 22:19:19 -05002055 return bits;
dburgessb3a0ca42011-10-12 07:44:40 +00002056}
Thomas Tsou94edaae2013-11-09 22:19:19 -05002057
Tom Tsoud3253432016-03-06 03:08:01 -08002058/*
2059 * Demodulate an 8-PSK burst. Prior to symbol rotation, operate at
2060 * 4 SPS (if activated) to minimize distortion through the fractional
2061 * delay filters. Symbol rotation and after always operates at 1 SPS.
2062 *
2063 * Allow 1 SPS demodulation here, but note that other parts of the
2064 * transceiver restrict EDGE operatoin to 4 SPS - 8-PSK distortion
2065 * through the fractional delay filters at 1 SPS renders signal
2066 * nearly unrecoverable.
2067 */
2068SoftVector *demodEdgeBurst(signalVector &burst, int sps,
2069 complex chan, float toa)
2070{
2071 SoftVector *bits;
Tom Tsou7fec3032016-03-06 22:33:20 -08002072 signalVector *dec, *rot, *eq;
Tom Tsoud3253432016-03-06 03:08:01 -08002073
Tom Tsou7fec3032016-03-06 22:33:20 -08002074 dec = demodCommon(burst, sps, chan, toa);
2075 if (!dec)
Tom Tsoud3253432016-03-06 03:08:01 -08002076 return NULL;
2077
Tom Tsou7fec3032016-03-06 22:33:20 -08002078 /* Equalize and derotate */
Tom Tsoud3253432016-03-06 03:08:01 -08002079 eq = convolve(dec, GSMPulse4->c0_inv, NULL, NO_DELAY);
2080 rot = derotateEdgeBurst(*eq, 1);
2081
Tom Tsou7fec3032016-03-06 22:33:20 -08002082 /* Soft slice and normalize */
Tom Tsoud3253432016-03-06 03:08:01 -08002083 bits = softSliceEdgeBurst(*dec);
2084 vectorSlicer(bits);
2085
2086 delete dec;
2087 delete eq;
2088 delete rot;
2089
2090 return bits;
2091}
2092
Tom Tsou2079a3c2016-03-06 00:58:56 -08002093bool sigProcLibSetup()
Thomas Tsouc1f7c422013-10-11 13:49:55 -04002094{
Thomas Tsouc1f7c422013-10-11 13:49:55 -04002095 initTrigTables();
Thomas Tsou0e0e1f42013-11-09 22:08:51 -05002096 generateSincTable();
Tom Tsou2079a3c2016-03-06 00:58:56 -08002097 initGMSKRotationTables();
Thomas Tsouc1f7c422013-10-11 13:49:55 -04002098
Tom Tsou2079a3c2016-03-06 00:58:56 -08002099 GSMPulse1 = generateGSMPulse(1);
2100 GSMPulse4 = generateGSMPulse(4);
Thomas Tsouc1f7c422013-10-11 13:49:55 -04002101
Tom Tsou2079a3c2016-03-06 00:58:56 -08002102 generateRACHSequence(1);
Tom Tsoud3253432016-03-06 03:08:01 -08002103 for (int tsc = 0; tsc < 8; tsc++) {
Tom Tsou2079a3c2016-03-06 00:58:56 -08002104 generateMidamble(1, tsc);
Tom Tsoud3253432016-03-06 03:08:01 -08002105 gEdgeMidambles[tsc] = generateEdgeMidamble(tsc);
2106 }
Thomas Tsouc1f7c422013-10-11 13:49:55 -04002107
Thomas Tsouf79c4d02013-11-09 15:51:56 -06002108 generateDelayFilters();
2109
Tom Tsoud3253432016-03-06 03:08:01 -08002110 dnsampler = new Resampler(1, 4);
2111 if (!dnsampler->init()) {
2112 LOG(ALERT) << "Rx resampler failed to initialize";
2113 goto fail;
2114 }
2115
2116 dnsampler->enableHistory(false);
2117 dnsampler_in = new signalVector(DOWNSAMPLE_IN_LEN, dnsampler->len());
2118
Thomas Tsouc1f7c422013-10-11 13:49:55 -04002119 return true;
Tom Tsoud3253432016-03-06 03:08:01 -08002120
2121fail:
2122 sigProcLibDestroy();
2123 return false;
Thomas Tsouc1f7c422013-10-11 13:49:55 -04002124}