Mercurial > hg > audiostuff

comparison spandsp-0.0.6pre17/src/v17rx.c @ 4:26cd8f1ef0b1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
import spandsp-0.0.6pre17
author Peter Meerwald <pmeerw@cosy.sbg.ac.at>
date Fri, 25 Jun 2010 15:50:58 +0200
parents
children
comparison
equal deleted inserted replaced
3:c6c5a16ce2f2 4:26cd8f1ef0b1
1 #define IAXMODEM_STUFF
2 /*
3 * SpanDSP - a series of DSP components for telephony
4 *
5 * v17rx.c - ITU V.17 modem receive part
6 *
7 * Written by Steve Underwood <steveu@coppice.org>
8 *
9 * Copyright (C) 2004, 2005, 2006, 2007 Steve Underwood
10 *
11 * All rights reserved.
12 *
13 * This program is free software; you can redistribute it and/or modify
14 * it under the terms of the GNU Lesser General Public License version 2.1,
15 * as published by the Free Software Foundation.
16 *
17 * This program is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU Lesser General Public License for more details.
21 *
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with this program; if not, write to the Free Software
24 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
25 *
26 * $Id: v17rx.c,v 1.153.4.6 2009/12/28 12:20:46 steveu Exp $
27 */
28
29 /*! \file */
30
31 #if defined(HAVE_CONFIG_H)
32 #include "config.h"
33 #endif
34
35 #include <stdlib.h>
36 #include <inttypes.h>
37 #include <string.h>
38 #include <stdio.h>
39 #if defined(HAVE_TGMATH_H)
40 #include <tgmath.h>
41 #endif
42 #if defined(HAVE_MATH_H)
43 #include <math.h>
44 #endif
45 #include "floating_fudge.h"
46
47 #include "spandsp/telephony.h"
48 #include "spandsp/logging.h"
49 #include "spandsp/complex.h"
50 #include "spandsp/vector_float.h"
51 #include "spandsp/complex_vector_float.h"
52 #include "spandsp/vector_int.h"
53 #include "spandsp/complex_vector_int.h"
54 #include "spandsp/async.h"
55 #include "spandsp/power_meter.h"
56 #include "spandsp/arctan2.h"
57 #include "spandsp/dds.h"
58 #include "spandsp/complex_filters.h"
59
60 #include "spandsp/v29rx.h"
61 #include "spandsp/v17tx.h"
62 #include "spandsp/v17rx.h"
63
64 #include "spandsp/private/logging.h"
65 #include "spandsp/private/v17rx.h"
66
67 #include "v17_v32bis_tx_constellation_maps.h"
68 #include "v17_v32bis_rx_constellation_maps.h"
69 #if defined(SPANDSP_USE_FIXED_POINT)
70 #include "v17_v32bis_rx_fixed_rrc.h"
71 #else
72 #include "v17_v32bis_rx_floating_rrc.h"
73 #endif
74
75 /*! The nominal frequency of the carrier, in Hertz */
76 #define CARRIER_NOMINAL_FREQ 1800.0f
77 /*! The nominal baud or symbol rate */
78 #define BAUD_RATE 2400
79 /*! The adaption rate coefficient for the equalizer during initial training */
80 #define EQUALIZER_DELTA 0.21f
81 /*! The adaption rate coefficient for the equalizer during continuous fine tuning */
82 #define EQUALIZER_SLOW_ADAPT_RATIO 0.1f
83
84 /* Segments of the training sequence */
85 /*! The length of training segment 1, in symbols */
86 #define V17_TRAINING_SEG_1_LEN 256
87 /*! The length of training segment 2 in long training mode, in symbols */
88 #define V17_TRAINING_SEG_2_LEN 2976
89 /*! The length of training segment 2 in short training mode, in symbols */
90 #define V17_TRAINING_SHORT_SEG_2_LEN 38
91 /*! The length of training segment 3, in symbols */
92 #define V17_TRAINING_SEG_3_LEN 64
93 /*! The length of training segment 4A, in symbols */
94 #define V17_TRAINING_SEG_4A_LEN 15
95 /*! The length of training segment 4, in symbols */
96 #define V17_TRAINING_SEG_4_LEN 48
97
98 /*! The 16 bit pattern used in the bridge section of the training sequence */
99 #define V17_BRIDGE_WORD 0x8880
100
101 /*! The length of the equalizer buffer */
102 #define V17_EQUALIZER_LEN (V17_EQUALIZER_PRE_LEN + 1 + V17_EQUALIZER_POST_LEN)
103
104 enum
105 {
106 TRAINING_STAGE_NORMAL_OPERATION = 0,
107 TRAINING_STAGE_SYMBOL_ACQUISITION,
108 TRAINING_STAGE_LOG_PHASE,
109 TRAINING_STAGE_SHORT_WAIT_FOR_CDBA,
110 TRAINING_STAGE_WAIT_FOR_CDBA,
111 TRAINING_STAGE_COARSE_TRAIN_ON_CDBA,
112 TRAINING_STAGE_FINE_TRAIN_ON_CDBA,
113 TRAINING_STAGE_SHORT_TRAIN_ON_CDBA_AND_TEST,
114 TRAINING_STAGE_TRAIN_ON_CDBA_AND_TEST,
115 TRAINING_STAGE_BRIDGE,
116 TRAINING_STAGE_TCM_WINDUP,
117 TRAINING_STAGE_TEST_ONES,
118 TRAINING_STAGE_PARKED
119 };
120
121 /* Coefficients for the band edge symbol timing synchroniser (alpha = 0.99) */
122 /* low_edge = 2.0f*M_PI*(CARRIER_NOMINAL_FREQ - BAUD_RATE/2.0f)/SAMPLE_RATE; */
123 /* high_edge = 2.0f*M_PI*(CARRIER_NOMINAL_FREQ + BAUD_RATE/2.0f)/SAMPLE_RATE; */
124 #define SIN_LOW_BAND_EDGE 0.453990499f
125 #define COS_LOW_BAND_EDGE 0.891006542f
126 #define SIN_HIGH_BAND_EDGE 0.707106781f
127 #define COS_HIGH_BAND_EDGE -0.707106781f
128 #define ALPHA 0.99f
129
130 #if defined(SPANDSP_USE_FIXED_POINTx)
131 #define SYNC_LOW_BAND_EDGE_COEFF_0 ((int)(FP_FACTOR*(2.0f*ALPHA*COS_LOW_BAND_EDGE)))
132 #define SYNC_LOW_BAND_EDGE_COEFF_1 ((int)(FP_FACTOR*(-ALPHA*ALPHA)))
133 #define SYNC_LOW_BAND_EDGE_COEFF_2 ((int)(FP_FACTOR*(-ALPHA*SIN_LOW_BAND_EDGE)))
134 #define SYNC_HIGH_BAND_EDGE_COEFF_0 ((int)(FP_FACTOR*(2.0f*ALPHA*COS_HIGH_BAND_EDGE)))
135 #define SYNC_HIGH_BAND_EDGE_COEFF_1 ((int)(FP_FACTOR*(-ALPHA*ALPHA)))
136 #define SYNC_HIGH_BAND_EDGE_COEFF_2 ((int)(FP_FACTOR*(-ALPHA*SIN_HIGH_BAND_EDGE)))
137 #define SYNC_MIXED_EDGES_COEFF_3 ((int)(FP_FACTOR*(-ALPHA*ALPHA*(SIN_HIGH_BAND_EDGE*COS_LOW_BAND_EDGE - SIN_LOW_BAND_EDGE*COS_HIGH_BAND_EDGE))))
138 #else
139 #define SYNC_LOW_BAND_EDGE_COEFF_0 (2.0f*ALPHA*COS_LOW_BAND_EDGE)
140 #define SYNC_LOW_BAND_EDGE_COEFF_1 (-ALPHA*ALPHA)
141 #define SYNC_LOW_BAND_EDGE_COEFF_2 (-ALPHA*SIN_LOW_BAND_EDGE)
142 #define SYNC_HIGH_BAND_EDGE_COEFF_0 (2.0f*ALPHA*COS_HIGH_BAND_EDGE)
143 #define SYNC_HIGH_BAND_EDGE_COEFF_1 (-ALPHA*ALPHA)
144 #define SYNC_HIGH_BAND_EDGE_COEFF_2 (-ALPHA*SIN_HIGH_BAND_EDGE)
145 #define SYNC_MIXED_EDGES_COEFF_3 (-ALPHA*ALPHA*(SIN_HIGH_BAND_EDGE*COS_LOW_BAND_EDGE - SIN_LOW_BAND_EDGE*COS_HIGH_BAND_EDGE))
146 #endif
147
148 #if defined(SPANDSP_USE_FIXED_POINTx)
149 static const int constellation_spacing[4] =
150 {
151 ((int)(FP_FACTOR*1.414f),
152 ((int)(FP_FACTOR*2.0f)},
153 ((int)(FP_FACTOR*2.828f)},
154 ((int)(FP_FACTOR*4.0f)},
155 };
156 #else
157 static const float constellation_spacing[4] =
158 {
159 1.414f,
160 2.0f,
161 2.828f,
162 4.0f
163 };
164 #endif
165
166 SPAN_DECLARE(float) v17_rx_carrier_frequency(v17_rx_state_t *s)
167 {
168 return dds_frequencyf(s->carrier_phase_rate);
169 }
170 /*- End of function --------------------------------------------------------*/
171
172 SPAN_DECLARE(float) v17_rx_symbol_timing_correction(v17_rx_state_t *s)
173 {
174 return (float) s->total_baud_timing_correction/((float) RX_PULSESHAPER_COEFF_SETS*10.0f/3.0f);
175 }
176 /*- End of function --------------------------------------------------------*/
177
178 SPAN_DECLARE(float) v17_rx_signal_power(v17_rx_state_t *s)
179 {
180 return power_meter_current_dbm0(&s->power) + 3.98f;
181 }
182 /*- End of function --------------------------------------------------------*/
183
184 SPAN_DECLARE(void) v17_rx_signal_cutoff(v17_rx_state_t *s, float cutoff)
185 {
186 /* The 0.4 factor allows for the gain of the DC blocker */
187 s->carrier_on_power = (int32_t) (power_meter_level_dbm0(cutoff + 2.5f)*0.4f);
188 s->carrier_off_power = (int32_t) (power_meter_level_dbm0(cutoff - 2.5f)*0.4f);
189 }
190 /*- End of function --------------------------------------------------------*/
191
192 static void report_status_change(v17_rx_state_t *s, int status)
193 {
194 if (s->status_handler)
195 s->status_handler(s->status_user_data, status);
196 else if (s->put_bit)
197 s->put_bit(s->put_bit_user_data, status);
198 }
199 /*- End of function --------------------------------------------------------*/
200
201 #if defined(SPANDSP_USE_FIXED_POINTx)
202 SPAN_DECLARE(int) v17_rx_equalizer_state(v17_rx_state_t *s, complexi16_t **coeffs)
203 #else
204 SPAN_DECLARE(int) v17_rx_equalizer_state(v17_rx_state_t *s, complexf_t **coeffs)
205 #endif
206 {
207 *coeffs = s->eq_coeff;
208 return V17_EQUALIZER_LEN;
209 }
210 /*- End of function --------------------------------------------------------*/
211
212 static void equalizer_save(v17_rx_state_t *s)
213 {
214 #if defined(SPANDSP_USE_FIXED_POINTx)
215 cvec_copyi16(s->eq_coeff_save, s->eq_coeff, V17_EQUALIZER_LEN);
216 #else
217 cvec_copyf(s->eq_coeff_save, s->eq_coeff, V17_EQUALIZER_LEN);
218 #endif
219 }
220 /*- End of function --------------------------------------------------------*/
221
222 static void equalizer_restore(v17_rx_state_t *s)
223 {
224 #if defined(SPANDSP_USE_FIXED_POINTx)
225 cvec_copyi16(s->eq_coeff, s->eq_coeff_save, V17_EQUALIZER_LEN);
226 cvec_zeroi16(s->eq_buf, V17_EQUALIZER_LEN);
227 s->eq_delta = 32768.0f*EQUALIZER_SLOW_ADAPT_RATIO*EQUALIZER_DELTA/V17_EQUALIZER_LEN;
228 #else
229 cvec_copyf(s->eq_coeff, s->eq_coeff_save, V17_EQUALIZER_LEN);
230 cvec_zerof(s->eq_buf, V17_EQUALIZER_LEN);
231 s->eq_delta = EQUALIZER_SLOW_ADAPT_RATIO*EQUALIZER_DELTA/V17_EQUALIZER_LEN;
232 #endif
233
234 s->eq_put_step = RX_PULSESHAPER_COEFF_SETS*10/(3*2) - 1;
235 s->eq_step = 0;
236 }
237 /*- End of function --------------------------------------------------------*/
238
239 static void equalizer_reset(v17_rx_state_t *s)
240 {
241 /* Start with an equalizer based on everything being perfect */
242 #if defined(SPANDSP_USE_FIXED_POINTx)
243 cvec_zeroi16(s->eq_coeff, V17_EQUALIZER_LEN);
244 s->eq_coeff[V17_EQUALIZER_PRE_LEN] = complex_seti16(3*FP_FACTOR, 0);
245 cvec_zeroi16(s->eq_buf, V17_EQUALIZER_LEN);
246 s->eq_delta = 32768.0f*EQUALIZER_DELTA/V17_EQUALIZER_LEN;
247 #else
248 cvec_zerof(s->eq_coeff, V17_EQUALIZER_LEN);
249 s->eq_coeff[V17_EQUALIZER_PRE_LEN] = complex_setf(3.0f, 0.0f);
250 cvec_zerof(s->eq_buf, V17_EQUALIZER_LEN);
251 s->eq_delta = EQUALIZER_DELTA/V17_EQUALIZER_LEN;
252 #endif
253
254 s->eq_put_step = RX_PULSESHAPER_COEFF_SETS*10/(3*2) - 1;
255 s->eq_step = 0;
256 }
257 /*- End of function --------------------------------------------------------*/
258
259 #if defined(SPANDSP_USE_FIXED_POINTx)
260 static __inline__ complexi16_t equalizer_get(v17_rx_state_t *s)
261 #else
262 static __inline__ complexf_t equalizer_get(v17_rx_state_t *s)
263 #endif
264 {
265 return cvec_circular_dot_prodf(s->eq_buf, s->eq_coeff, V17_EQUALIZER_LEN, s->eq_step);
266 }
267 /*- End of function --------------------------------------------------------*/
268
269 #if defined(SPANDSP_USE_FIXED_POINTx)
270 static void tune_equalizer(v17_rx_state_t *s, const complexi16_t *z, const complexi16_t *target)
271 {
272 complexi16_t err;
273
274 /* Find the x and y mismatch from the exact constellation position. */
275 err.re = target->re*FP_FACTOR - z->re;
276 err.im = target->im*FP_FACTOR - z->im;
277 //span_log(&s->logging, SPAN_LOG_FLOW, "Equalizer error %f\n", sqrt(err.re*err.re + err.im*err.im));
278 err.re = ((int32_t) err.re*(int32_t) s->eq_delta) >> 15;
279 err.im = ((int32_t) err.im*(int32_t) s->eq_delta) >> 15;
280 cvec_circular_lmsi16(s->eq_buf, s->eq_coeff, V17_EQUALIZER_LEN, s->eq_step, &err);
281 }
282 #else
283 static void tune_equalizer(v17_rx_state_t *s, const complexf_t *z, const complexf_t *target)
284 {
285 complexf_t err;
286
287 /* Find the x and y mismatch from the exact constellation position. */
288 err = complex_subf(target, z);
289 //span_log(&s->logging, SPAN_LOG_FLOW, "Equalizer error %f\n", sqrt(err.re*err.re + err.im*err.im));
290 err.re *= s->eq_delta;
291 err.im *= s->eq_delta;
292 cvec_circular_lmsf(s->eq_buf, s->eq_coeff, V17_EQUALIZER_LEN, s->eq_step, &err);
293 }
294 #endif
295
296 static int descramble(v17_rx_state_t *s, int in_bit)
297 {
298 int out_bit;
299
300 //out_bit = (in_bit ^ (s->scramble_reg >> s->scrambler_tap) ^ (s->scramble_reg >> (23 - 1))) & 1;
301 out_bit = (in_bit ^ (s->scramble_reg >> (18 - 1)) ^ (s->scramble_reg >> (23 - 1))) & 1;
302 s->scramble_reg <<= 1;
303 if (s->training_stage > TRAINING_STAGE_NORMAL_OPERATION && s->training_stage < TRAINING_STAGE_TCM_WINDUP)
304 s->scramble_reg |= out_bit;
305 else
306 s->scramble_reg |= (in_bit & 1);
307 return out_bit;
308 }
309 /*- End of function --------------------------------------------------------*/
310
311 static void track_carrier(v17_rx_state_t *s, const complexf_t *z, const complexf_t *target)
312 {
313 float error;
314
315 /* For small errors the imaginary part of the difference between the actual and the target
316 positions is proportional to the phase error, for any particular target. However, the
317 different amplitudes of the various target positions scale things. */
318 error = z->im*target->re - z->re*target->im;
319
320 s->carrier_phase_rate += (int32_t) (s->carrier_track_i*error);
321 s->carrier_phase += (int32_t) (s->carrier_track_p*error);
322 //span_log(&s->logging, SPAN_LOG_FLOW, "Im = %15.5f f = %15.5f\n", error, dds_frequencyf(s->carrier_phase_rate));
323 //printf("XXX Im = %15.5f f = %15.5f %f %f %f %f (%f %f)\n", error, dds_frequencyf(s->carrier_phase_rate), target->re, target->im, z->re, z->im, s->carrier_track_i, s->carrier_track_p);
324 }
325 /*- End of function --------------------------------------------------------*/
326
327 static __inline__ void put_bit(v17_rx_state_t *s, int bit)
328 {
329 int out_bit;
330
331 /* We need to strip the last part of the training - the test period of all 1s -
332 before we let data go to the application. */
333 if (s->training_stage == TRAINING_STAGE_NORMAL_OPERATION)
334 {
335 out_bit = descramble(s, bit);
336 s->put_bit(s->put_bit_user_data, out_bit);
337 }
338 else if (s->training_stage == TRAINING_STAGE_TEST_ONES)
339 {
340 /* The bits during the final stage of training should be all ones. However,
341 buggy modems mean you cannot rely on this. Therefore we don't bother
342 testing for ones, but just rely on a constellation mismatch measurement. */
343 out_bit = descramble(s, bit);
344 //span_log(&s->logging, SPAN_LOG_FLOW, "A 1 is really %d\n", out_bit);
345 }
346 }
347 /*- End of function --------------------------------------------------------*/
348
349 #if defined(SPANDSP_USE_FIXED_POINTx)
350 static __inline__ uint32_t dist_sq(const complexi_t *x, const complexi_t *y)
351 {
352 return (x->re - y->re)*(x->re - y->re) + (x->im - y->im)*(x->im - y->im);
353 }
354 /*- End of function --------------------------------------------------------*/
355 #else
356 static __inline__ float dist_sq(const complexf_t *x, const complexf_t *y)
357 {
358 return (x->re - y->re)*(x->re - y->re) + (x->im - y->im)*(x->im - y->im);
359 }
360 /*- End of function --------------------------------------------------------*/
361 #endif
362
363 static int decode_baud(v17_rx_state_t *s, complexf_t *z)
364 {
365 static const uint8_t v32bis_4800_differential_decoder[4][4] =
366 {
367 {2, 3, 0, 1},
368 {0, 2, 1, 3},
369 {3, 1, 2, 0},
370 {1, 0, 3, 2}
371 };
372 static const uint8_t v17_differential_decoder[4][4] =
373 {
374 {0, 1, 2, 3},
375 {3, 0, 1, 2},
376 {2, 3, 0, 1},
377 {1, 2, 3, 0}
378 };
379 static const uint8_t tcm_paths[8][4] =
380 {
381 {0, 6, 2, 4},
382 {6, 0, 4, 2},
383 {2, 4, 0, 6},
384 {4, 2, 6, 0},
385 {1, 3, 7, 5},
386 {5, 7, 3, 1},
387 {7, 5, 1, 3},
388 {3, 1, 5, 7}
389 };
390 int nearest;
391 int i;
392 int j;
393 int k;
394 int re;
395 int im;
396 int raw;
397 int constellation_state;
398 #if defined(SPANDSP_USE_FIXED_POINTx)
399 #define DIST_FACTOR 2048 /* Something less than sqrt(0xFFFFFFFF/10)/10 */
400 complexi_t zi;
401 uint32_t distances[8];
402 uint32_t new_distances[8];
403 uint32_t min;
404 complexi_t ci;
405 #else
406 float distances[8];
407 float new_distances[8];
408 float min;
409 #endif
410
411 re = (int) ((z->re + 9.0f)*2.0f);
412 if (re > 35)
413 re = 35;
414 else if (re < 0)
415 re = 0;
416 im = (int) ((z->im + 9.0f)*2.0f);
417 if (im > 35)
418 im = 35;
419 else if (im < 0)
420 im = 0;
421
422 if (s->bits_per_symbol == 2)
423 {
424 /* 4800bps V.32bis mode, without trellis coding */
425 nearest = constel_map_4800[re][im];
426 raw = v32bis_4800_differential_decoder[s->diff][nearest];
427 s->diff = nearest;
428 put_bit(s, raw);
429 put_bit(s, raw >> 1);
430 return nearest;
431 }
432
433 /* Find a set of 8 candidate constellation positions, that are the closest
434 to the target, with different patterns in the last 3 bits. */
435 #if defined(SPANDSP_USE_FIXED_POINTx)
436 min = 0xFFFFFFFF;
437 zi = complex_seti(z->re*DIST_FACTOR, z->im*DIST_FACTOR);
438 #else
439 min = 9999999.0f;
440 #endif
441 j = 0;
442 for (i = 0; i < 8; i++)
443 {
444 nearest = constel_maps[s->space_map][re][im][i];
445 #if defined(SPANDSP_USE_FIXED_POINTx)
446 ci = complex_seti(s->constellation[nearest].re*DIST_FACTOR,
447 s->constellation[nearest].im*DIST_FACTOR);
448 distances[i] = dist_sq(&ci, &zi);
449 #else
450 distances[i] = dist_sq(&s->constellation[nearest], z);
451 #endif
452 if (min > distances[i])
453 {
454 min = distances[i];
455 j = i;
456 }
457 }
458 /* Use the nearest of these soft-decisions as the basis for DFE */
459 constellation_state = constel_maps[s->space_map][re][im][j];
460 /* Control the equalizer, carrier tracking, etc. based on the non-trellis
461 corrected information. The trellis correct stuff comes out a bit late. */
462 track_carrier(s, z, &s->constellation[constellation_state]);
463 //tune_equalizer(s, z, &s->constellation[constellation_state]);
464
465 /* Now do the trellis decoding */
466
467 /* TODO: change to processing blocks of stored symbols here, instead of processing
468 one symbol at a time, to speed up the processing. */
469
470 /* Update the minimum accumulated distance to each of the 8 states */
471 if (++s->trellis_ptr >= V17_TRELLIS_STORAGE_DEPTH)
472 s->trellis_ptr = 0;
473 for (i = 0; i < 4; i++)
474 {
475 min = distances[tcm_paths[i][0]] + s->distances[0];
476 k = 0;
477 for (j = 1; j < 4; j++)
478 {
479 if (min > distances[tcm_paths[i][j]] + s->distances[j << 1])
480 {
481 min = distances[tcm_paths[i][j]] + s->distances[j << 1];
482 k = j;
483 }
484 }
485 /* Use an elementary IIR filter to track the distance to date. */
486 #if defined(SPANDSP_USE_FIXED_POINTx)
487 new_distances[i] = s->distances[k << 1]*9/10 + distances[tcm_paths[i][k]]*1/10;
488 #else
489 new_distances[i] = s->distances[k << 1]*0.9f + distances[tcm_paths[i][k]]*0.1f;
490 #endif
491 s->full_path_to_past_state_locations[s->trellis_ptr][i] = constel_maps[s->space_map][re][im][tcm_paths[i][k]];
492 s->past_state_locations[s->trellis_ptr][i] = k << 1;
493 }
494 for (i = 4; i < 8; i++)
495 {
496 min = distances[tcm_paths[i][0]] + s->distances[1];
497 k = 0;
498 for (j = 1; j < 4; j++)
499 {
500 if (min > distances[tcm_paths[i][j]] + s->distances[(j << 1) + 1])
501 {
502 min = distances[tcm_paths[i][j]] + s->distances[(j << 1) + 1];
503 k = j;
504 }
505 }
506 #if defined(SPANDSP_USE_FIXED_POINTx)
507 new_distances[i] = s->distances[(k << 1) + 1]*9/10 + distances[tcm_paths[i][k]]*1/10;
508 #else
509 new_distances[i] = s->distances[(k << 1) + 1]*0.9f + distances[tcm_paths[i][k]]*0.1f;
510 #endif
511 s->full_path_to_past_state_locations[s->trellis_ptr][i] = constel_maps[s->space_map][re][im][tcm_paths[i][k]];
512 s->past_state_locations[s->trellis_ptr][i] = (k << 1) + 1;
513 }
514 memcpy(s->distances, new_distances, sizeof(s->distances));
515
516 /* Find the minimum distance to date. This is the start of the path back to the result. */
517 min = s->distances[0];
518 k = 0;
519 for (i = 1; i < 8; i++)
520 {
521 if (min > s->distances[i])
522 {
523 min = s->distances[i];
524 k = i;
525 }
526 }
527 /* Trace back through every time step, starting with the current one, and find the
528 state from which the path came one step before. At the end of this search, the
529 last state found also points to the constellation point at that state. This is the
530 output of the trellis. */
531 for (i = 0, j = s->trellis_ptr; i < V17_TRELLIS_LOOKBACK_DEPTH - 1; i++)
532 {
533 k = s->past_state_locations[j][k];
534 if (--j < 0)
535 j = V17_TRELLIS_STORAGE_DEPTH - 1;
536 }
537 nearest = s->full_path_to_past_state_locations[j][k] >> 1;
538
539 /* Differentially decode */
540 raw = (nearest & 0x3C) | v17_differential_decoder[s->diff][nearest & 0x03];
541 s->diff = nearest & 0x03;
542 for (i = 0; i < s->bits_per_symbol; i++)
543 {
544 put_bit(s, raw);
545 raw >>= 1;
546 }
547 return constellation_state;
548 }
549 /*- End of function --------------------------------------------------------*/
550
551 static __inline__ void symbol_sync(v17_rx_state_t *s)
552 {
553 int i;
554 #if defined(SPANDSP_USE_FIXED_POINTx)
555 int32_t v;
556 int32_t p;
557 #else
558 float v;
559 float p;
560 #endif
561
562 /* This routine adapts the position of the half baud samples entering the equalizer. */
563
564 /* This symbol sync scheme is based on the technique first described by Dominique Godard in
565 Passband Timing Recovery in an All-Digital Modem Receiver
566 IEEE TRANSACTIONS ON COMMUNICATIONS, VOL. COM-26, NO. 5, MAY 1978 */
567
568 /* This is slightly rearranged for figure 3b of the Godard paper, as this saves a couple of
569 maths operations */
570 #if defined(SPANDSP_USE_FIXED_POINTx)
571 /* TODO: The scalings used here need more thorough evaluation, to see if overflows are possible. */
572 /* Cross correlate */
573 v = (((s->symbol_sync_low[1] >> 5)*(s->symbol_sync_high[0] >> 4)) >> 15)*SYNC_LOW_BAND_EDGE_COEFF_2
574 - (((s->symbol_sync_low[0] >> 5)*(s->symbol_sync_high[1] >> 4)) >> 15)*SYNC_HIGH_BAND_EDGE_COEFF_2
575 + (((s->symbol_sync_low[1] >> 5)*(s->symbol_sync_high[1] >> 4)) >> 15)*SYNC_MIXED_EDGES_COEFF_3;
576 /* Filter away any DC component */
577 p = v - s->symbol_sync_dc_filter[1];
578 s->symbol_sync_dc_filter[1] = s->symbol_sync_dc_filter[0];
579 s->symbol_sync_dc_filter[0] = v;
580 /* A little integration will now filter away much of the HF noise */
581 s->baud_phase -= p;
582 if (abs(s->baud_phase) > 100*FP_FACTOR)
583 {
584 if (s->baud_phase > 0)
585 i = (s->baud_phase > 1000*FP_FACTOR) ? 15 : 1;
586 else
587 i = (s->baud_phase < -1000*FP_FACTOR) ? -15 : -1;
588 //printf("v = %10.5f %5d - %f %f %d %d\n", v, i, p, s->baud_phase, s->total_baud_timing_correction);
589 s->eq_put_step += i;
590 s->total_baud_timing_correction += i;
591 }
592 #else
593 /* Cross correlate */
594 v = s->symbol_sync_low[1]*s->symbol_sync_high[0]*SYNC_LOW_BAND_EDGE_COEFF_2
595 - s->symbol_sync_low[0]*s->symbol_sync_high[1]*SYNC_HIGH_BAND_EDGE_COEFF_2
596 + s->symbol_sync_low[1]*s->symbol_sync_high[1]*SYNC_MIXED_EDGES_COEFF_3;
597 /* Filter away any DC component */
598 p = v - s->symbol_sync_dc_filter[1];
599 s->symbol_sync_dc_filter[1] = s->symbol_sync_dc_filter[0];
600 s->symbol_sync_dc_filter[0] = v;
601 /* A little integration will now filter away much of the HF noise */
602 s->baud_phase -= p;
603 if (fabsf(s->baud_phase) > 100.0f)
604 {
605 if (s->baud_phase > 0.0f)
606 i = (s->baud_phase > 1000.0f) ? 15 : 1;
607 else
608 i = (s->baud_phase < -1000.0f) ? -15 : -1;
609 //printf("v = %10.5f %5d - %f %f %d\n", v, i, p, s->baud_phase, s->total_baud_timing_correction);
610 s->eq_put_step += i;
611 s->total_baud_timing_correction += i;
612 }
613 #endif
614 }
615 /*- End of function --------------------------------------------------------*/
616
617 static void process_half_baud(v17_rx_state_t *s, const complexf_t *sample)
618 {
619 static const complexf_t cdba[4] =
620 {
621 { 6.0f, 2.0f},
622 {-2.0f, 6.0f},
623 { 2.0f, -6.0f},
624 {-6.0f, -2.0f}
625 };
626 complexf_t z;
627 complexf_t zz;
628 #if defined(SPANDSP_USE_FIXED_POINTx)
629 const complexi_t *target;
630 static const complexi16_t zero = {0, 0};
631 #else
632 const complexf_t *target;
633 static const complexf_t zero = {0, 0};
634 #endif
635 float p;
636 int bit;
637 int i;
638 int j;
639 int32_t angle;
640 int32_t ang;
641 int constellation_state;
642
643 /* This routine processes every half a baud, as we put things into the equalizer at the T/2 rate. */
644
645 /* Add a sample to the equalizer's circular buffer, but don't calculate anything at this time. */
646 s->eq_buf[s->eq_step] = *sample;
647 if (++s->eq_step >= V17_EQUALIZER_LEN)
648 s->eq_step = 0;
649
650 /* On alternate insertions we have a whole baud and must process it. */
651 if ((s->baud_half ^= 1))
652 return;
653
654 /* Symbol timing synchronisation */
655 symbol_sync(s);
656
657 z = equalizer_get(s);
658
659 constellation_state = 0;
660 switch (s->training_stage)
661 {
662 case TRAINING_STAGE_NORMAL_OPERATION:
663 /* Normal operation. */
664 constellation_state = decode_baud(s, &z);
665 target = &s->constellation[constellation_state];
666 break;
667 case TRAINING_STAGE_SYMBOL_ACQUISITION:
668 /* Allow time for the symbol synchronisation to settle the symbol timing. */
669 target = &zero;
670 if (++s->training_count >= 100)
671 {
672 /* Record the current phase angle */
673 s->angles[0] =
674 s->start_angles[0] = arctan2(z.im, z.re);
675 s->training_stage = TRAINING_STAGE_LOG_PHASE;
676 if (s->agc_scaling_save == 0.0f)
677 s->agc_scaling_save = s->agc_scaling;
678 }
679 break;
680 case TRAINING_STAGE_LOG_PHASE:
681 /* Record the current alternate phase angle */
682 target = &zero;
683 angle = arctan2(z.im, z.re);
684 s->training_count = 1;
685 if (s->short_train)
686 {
687 /* We should already know the accurate carrier frequency. All we need to sort
688 out is the phase. */
689 /* Check if we just saw A or B */
690 if ((uint32_t) (angle - s->start_angles[0]) < 0x80000000U)
691 {
692 angle = s->start_angles[0];
693 s->angles[0] = 0xC0000000 + 219937506;
694 s->angles[1] = 0x80000000 + 219937506;
695 }
696 else
697 {
698 s->angles[0] = 0x80000000 + 219937506;
699 s->angles[1] = 0xC0000000 + 219937506;
700 }
701 /* Make a step shift in the phase, to pull it into line. We need to rotate the equalizer
702 buffer, as well as the carrier phase, for this to play out nicely. */
703 /* angle is now the difference between where A is, and where it should be */
704 p = 3.14159f + angle*2.0f*3.14159f/(65536.0f*65536.0f) - 0.321751f;
705 span_log(&s->logging, SPAN_LOG_FLOW, "Spin (short) by %.5f rads\n", p);
706 zz = complex_setf(cosf(p), -sinf(p));
707 for (i = 0; i < V17_EQUALIZER_LEN; i++)
708 s->eq_buf[i] = complex_mulf(&s->eq_buf[i], &zz);
709 s->carrier_phase += (0x80000000 + angle - 219937506);
710
711 s->carrier_track_p = 500000.0f;
712
713 s->training_stage = TRAINING_STAGE_SHORT_WAIT_FOR_CDBA;
714 }
715 else
716 {
717 s->angles[1] =
718 s->start_angles[1] = angle;
719 s->training_stage = TRAINING_STAGE_WAIT_FOR_CDBA;
720 }
721 break;
722 case TRAINING_STAGE_WAIT_FOR_CDBA:
723 target = &zero;
724 angle = arctan2(z.im, z.re);
725 /* Look for the initial ABAB sequence to display a phase reversal, which will
726 signal the start of the scrambled CDBA segment */
727 ang = angle - s->angles[(s->training_count - 1) & 0xF];
728 s->angles[(s->training_count + 1) & 0xF] = angle;
729
730 /* Do a coarse frequency adjustment about half way through the reversals, as if we wait until
731 the end, we might have rotated too far to correct properly. */
732 if (s->training_count == 100)
733 {
734 i = s->training_count;
735 /* Avoid the possibility of a divide by zero */
736 if (i)
737 {
738 j = i & 0xF;
739 ang = (s->angles[j] - s->start_angles[0])/i
740 + (s->angles[j | 0x1] - s->start_angles[1])/i;
741 s->carrier_phase_rate += 3*(ang/20);
742 //span_log(&s->logging, SPAN_LOG_FLOW, "Angles %x, %x, %x, %x, dist %d\n", s->angles[j], s->start_angles[0], s->angles[j | 0x1], s->start_angles[1], i);
743
744 s->start_angles[0] = s->angles[j];
745 s->start_angles[1] = s->angles[j | 0x1];
746 }
747 //span_log(&s->logging, SPAN_LOG_FLOW, "%d %d %d %d %d\n", s->angles[s->training_count & 0xF], s->start_angles[0], s->angles[(s->training_count | 0x1) & 0xF], s->start_angles[1], s->training_count);
748 span_log(&s->logging, SPAN_LOG_FLOW, "First coarse carrier frequency %7.2f (%d)\n", dds_frequencyf(s->carrier_phase_rate), s->training_count);
749
750 }
751 if ((ang > 0x40000000 || ang < -0x40000000) && s->training_count >= 13)
752 {
753 span_log(&s->logging, SPAN_LOG_FLOW, "We seem to have a reversal at symbol %d\n", s->training_count);
754 /* We seem to have a phase reversal */
755 /* Slam the carrier frequency into line, based on the total phase drift over the last
756 section. Use the shift from the odd bits and the shift from the even bits to get
757 better jitter suppression. */
758 /* TODO: We are supposed to deal with frequancy errors up to +-8Hz. Over 200+
759 symbols that is more than half a cycle. We get confused an do crazy things.
760 We can only cope with errors up to 5Hz right now. We need to implement
761 greater tolerance to be compliant, although it doesn't really matter much
762 these days. */
763 /* Step back a few symbols so we don't get ISI distorting things. */
764 i = (s->training_count - 8) & ~1;
765 /* Avoid the possibility of a divide by zero */
766 if (i - 100 + 8)
767 {
768 j = i & 0xF;
769 ang = (s->angles[j] - s->start_angles[0])/(i - 100 + 8)
770 + (s->angles[j | 0x1] - s->start_angles[1])/(i - 100 + 8);
771 s->carrier_phase_rate += 3*(ang/20);
772 span_log(&s->logging, SPAN_LOG_FLOW, "Angles %x, %x, %x, %x, dist %d\n", s->angles[j], s->start_angles[0], s->angles[j | 0x1], s->start_angles[1], i);
773 }
774 //span_log(&s->logging, SPAN_LOG_FLOW, "%d %d %d %d %d\n", s->angles[s->training_count & 0xF], s->start_angles[0], s->angles[(s->training_count | 0x1) & 0xF], s->start_angles[1], s->training_count);
775 span_log(&s->logging, SPAN_LOG_FLOW, "Second coarse carrier frequency %7.2f (%d)\n", dds_frequencyf(s->carrier_phase_rate), s->training_count);
776 /* Check if the carrier frequency is plausible */
777 if (s->carrier_phase_rate < dds_phase_ratef(CARRIER_NOMINAL_FREQ - 20.0f)
778 ||
779 s->carrier_phase_rate > dds_phase_ratef(CARRIER_NOMINAL_FREQ + 20.0f))
780 {
781 span_log(&s->logging, SPAN_LOG_FLOW, "Training failed (sequence failed)\n");
782 /* Park this modem */
783 s->agc_scaling_save = 0.0f;
784 s->training_stage = TRAINING_STAGE_PARKED;
785 report_status_change(s, SIG_STATUS_TRAINING_FAILED);
786 break;
787 }
788
789 /* Make a step shift in the phase, to pull it into line. We need to rotate the equalizer buffer,
790 as well as the carrier phase, for this to play out nicely. */
791 /* angle is now the difference between where C is, and where it should be */
792 p = angle*2.0f*3.14159f/(65536.0f*65536.0f) - 0.321751f;
793 span_log(&s->logging, SPAN_LOG_FLOW, "Spin (long) by %.5f rads\n", p);
794 zz = complex_setf(cosf(p), -sinf(p));
795 for (i = 0; i < V17_EQUALIZER_LEN; i++)
796 s->eq_buf[i] = complex_mulf(&s->eq_buf[i], &zz);
797 s->carrier_phase += (angle - 219937506);
798
799 /* We have just seen the first symbol of the scrambled sequence, so skip it. */
800 bit = descramble(s, 1);
801 bit = (bit << 1) | descramble(s, 1);
802 target = &cdba[bit];
803 s->training_count = 1;
804 s->training_stage = TRAINING_STAGE_COARSE_TRAIN_ON_CDBA;
805 report_status_change(s, SIG_STATUS_TRAINING_IN_PROGRESS);
806 break;
807 }
808 if (++s->training_count > V17_TRAINING_SEG_1_LEN)
809 {
810 /* This is bogus. There are not this many bits in this section
811 of a real training sequence. Note that this might be TEP. */
812 span_log(&s->logging, SPAN_LOG_FLOW, "Training failed (sequence failed)\n");
813 /* Park this modem */
814 s->agc_scaling_save = 0.0f;
815 s->training_stage = TRAINING_STAGE_PARKED;
816 report_status_change(s, SIG_STATUS_TRAINING_FAILED);
817 }
818 break;
819 case TRAINING_STAGE_COARSE_TRAIN_ON_CDBA:
820 /* Train on the scrambled CDBA section. */
821 bit = descramble(s, 1);
822 bit = (bit << 1) | descramble(s, 1);
823 target = &cdba[bit];
824 track_carrier(s, &z, target);
825 tune_equalizer(s, &z, target);
826 #if defined(IAXMODEM_STUFF)
827 zz = complex_subf(&z, target);
828 s->training_error = powerf(&zz);
829 if (++s->training_count == V17_TRAINING_SEG_2_LEN - 2000 || s->training_error < 1.0f || s->training_error > 200.0f)
830 #else
831 if (++s->training_count == V17_TRAINING_SEG_2_LEN - 2000)
832 #endif
833 {
834 /* Now the equaliser adaption should be getting somewhere, slow it down, or it will never
835 tune very well on a noisy signal. */
836 s->eq_delta *= EQUALIZER_SLOW_ADAPT_RATIO;
837 s->carrier_track_i = 1000.0f;
838 s->training_stage = TRAINING_STAGE_FINE_TRAIN_ON_CDBA;
839 }
840 break;
841 case TRAINING_STAGE_FINE_TRAIN_ON_CDBA:
842 /* Train on the scrambled CDBA section. */
843 bit = descramble(s, 1);
844 bit = (bit << 1) | descramble(s, 1);
845 target = &cdba[bit];
846 /* By this point the training should be comming into focus. */
847 track_carrier(s, &z, target);
848 tune_equalizer(s, &z, target);
849 if (++s->training_count >= V17_TRAINING_SEG_2_LEN - 48)
850 {
851 s->training_error = 0.0f;
852 s->carrier_track_i = 100.0f;
853 s->carrier_track_p = 500000.0f;
854 s->training_stage = TRAINING_STAGE_TRAIN_ON_CDBA_AND_TEST;
855 }
856 break;
857 case TRAINING_STAGE_TRAIN_ON_CDBA_AND_TEST:
858 /* Continue training on the scrambled CDBA section, but measure the quality of training too. */
859 bit = descramble(s, 1);
860 bit = (bit << 1) | descramble(s, 1);
861 target = &cdba[bit];
862 //span_log(&s->logging, SPAN_LOG_FLOW, "%5d [%15.5f, %15.5f] [%15.5f, %15.5f]\n", s->training_count, z.re, z.im, cdba[bit].re, cdba[bit].im);
863 /* We ignore the last few symbols because it seems some modems do not end this
864 part properly, and it throws things off. */
865 if (++s->training_count < V17_TRAINING_SEG_2_LEN - 20)
866 {
867 track_carrier(s, &z, target);
868 tune_equalizer(s, &z, target);
869 /* Measure the training error */
870 zz = complex_subf(&z, &cdba[bit]);
871 s->training_error += powerf(&zz);
872 }
873 else if (s->training_count >= V17_TRAINING_SEG_2_LEN)
874 {
875 span_log(&s->logging, SPAN_LOG_FLOW, "Long training error %f\n", s->training_error);
876 if (s->training_error < 20.0f*1.414f*constellation_spacing[s->space_map])
877 {
878 s->training_count = 0;
879 s->training_error = 0.0f;
880 s->training_stage = TRAINING_STAGE_BRIDGE;
881 }
882 else
883 {
884 span_log(&s->logging, SPAN_LOG_FLOW, "Training failed (convergence failed)\n");
885 /* Park this modem */
886 s->agc_scaling_save = 0.0f;
887 s->training_stage = TRAINING_STAGE_PARKED;
888 report_status_change(s, SIG_STATUS_TRAINING_FAILED);
889 }
890 }
891 break;
892 case TRAINING_STAGE_BRIDGE:
893 descramble(s, V17_BRIDGE_WORD >> ((s->training_count & 0x7) << 1));
894 descramble(s, V17_BRIDGE_WORD >> (((s->training_count & 0x7) << 1) + 1));
895 target = &z;
896 if (++s->training_count >= V17_TRAINING_SEG_3_LEN)
897 {
898 s->training_count = 0;
899 s->training_error = 0.0f;
900 if (s->bits_per_symbol == 2)
901 {
902 /* Restart the differential decoder */
903 /* There is no trellis, so go straight to processing decoded data */
904 s->diff = (s->short_train) ? 0 : 1;
905 s->training_stage = TRAINING_STAGE_TEST_ONES;
906 }
907 else
908 {
909 /* Wait for the trellis to wind up */
910 s->training_stage = TRAINING_STAGE_TCM_WINDUP;
911 }
912 }
913 break;
914 case TRAINING_STAGE_SHORT_WAIT_FOR_CDBA:
915 /* Look for the initial ABAB sequence to display a phase reversal, which will
916 signal the start of the scrambled CDBA segment */
917 angle = arctan2(z.im, z.re);
918 ang = angle - s->angles[s->training_count & 1];
919 if (ang > 0x40000000 || ang < -0x40000000)
920 {
921 /* We seem to have a phase reversal */
922 /* We have just seen the first symbol of the scrambled sequence, so skip it. */
923 bit = descramble(s, 1);
924 bit = (bit << 1) | descramble(s, 1);
925 target = &cdba[bit];
926 s->training_count = 1;
927 s->training_error = 0.0f;
928 s->training_stage = TRAINING_STAGE_SHORT_TRAIN_ON_CDBA_AND_TEST;
929 break;
930 }
931 target = &cdba[(s->training_count & 1) + 2];
932 track_carrier(s, &z, target);
933 if (++s->training_count > V17_TRAINING_SEG_1_LEN)
934 {
935 /* This is bogus. There are not this many bits in this section
936 of a real training sequence. Note that this might be TEP. */
937 span_log(&s->logging, SPAN_LOG_FLOW, "Training failed (sequence failed)\n");
938 /* Park this modem */
939 s->training_stage = TRAINING_STAGE_PARKED;
940 report_status_change(s, SIG_STATUS_TRAINING_FAILED);
941 }
942 break;
943 case TRAINING_STAGE_SHORT_TRAIN_ON_CDBA_AND_TEST:
944 /* Short retrain on the scrambled CDBA section, but measure the quality of training too. */
945 bit = descramble(s, 1);
946 bit = (bit << 1) | descramble(s, 1);
947 //span_log(&s->logging, SPAN_LOG_FLOW, "%5d [%15.5f, %15.5f] [%15.5f, %15.5f] %d\n", s->training_count, z.re, z.im, cdba[bit].re, cdba[bit].im, arctan2(z.im, z.re));
948 target = &cdba[bit];
949 track_carrier(s, &z, target);
950 //tune_equalizer(s, &z, target);
951 /* Measure the training error */
952 if (s->training_count > 8)
953 {
954 zz = complex_subf(&z, &cdba[bit]);
955 s->training_error += powerf(&zz);
956 }
957 if (++s->training_count >= V17_TRAINING_SHORT_SEG_2_LEN)
958 {
959 span_log(&s->logging, SPAN_LOG_FLOW, "Short training error %f\n", s->training_error);
960 s->carrier_track_i = 100.0f;
961 s->carrier_track_p = 500000.0f;
962 /* TODO: This was increased by a factor of 10 after studying real world failures.
963 However, it is not clear why this is an improvement, If something gives
964 a huge training error, surely it shouldn't decode too well? */
965 if (s->training_error < (V17_TRAINING_SHORT_SEG_2_LEN - 8)*4.0f*constellation_spacing[s->space_map])
966 {
967 s->training_count = 0;
968 if (s->bits_per_symbol == 2)
969 {
970 /* There is no trellis, so go straight to processing decoded data */
971 /* Restart the differential decoder */
972 s->diff = (s->short_train) ? 0 : 1;
973 s->training_error = 0.0f;
974 s->training_stage = TRAINING_STAGE_TEST_ONES;
975 }
976 else
977 {
978 /* Wait for the trellis to wind up */
979 s->training_stage = TRAINING_STAGE_TCM_WINDUP;
980 }
981 report_status_change(s, SIG_STATUS_TRAINING_IN_PROGRESS);
982 }
983 else
984 {
985 span_log(&s->logging, SPAN_LOG_FLOW, "Short training failed (convergence failed)\n");
986 /* Park this modem */
987 s->training_stage = TRAINING_STAGE_PARKED;
988 report_status_change(s, SIG_STATUS_TRAINING_FAILED);
989 }
990 }
991 break;
992 case TRAINING_STAGE_TCM_WINDUP:
993 /* We need to wait 15 bauds while the trellis fills up. */
994 //span_log(&s->logging, SPAN_LOG_FLOW, "%5d %15.5f, %15.5f\n", s->training_count, z.re, z.im);
995 constellation_state = decode_baud(s, &z);
996 target = &s->constellation[constellation_state];
997 /* Measure the training error */
998 zz = complex_subf(&z, target);
999 s->training_error += powerf(&zz);
1000 if (++s->training_count >= V17_TRAINING_SEG_4A_LEN)
1001 {
1002 s->training_count = 0;
1003 s->training_error = 0.0f;
1004 /* Restart the differential decoder */
1005 s->diff = (s->short_train) ? 0 : 1;
1006 s->training_stage = TRAINING_STAGE_TEST_ONES;
1007 }
1008 break;
1009 case TRAINING_STAGE_TEST_ONES:
1010 /* We are in the test phase, where we check that we can receive reliably.
1011 We should get a run of 1's, 48 symbols long. */
1012 //span_log(&s->logging, SPAN_LOG_FLOW, "%5d %15.5f, %15.5f\n", s->training_count, z.re, z.im);
1013 constellation_state = decode_baud(s, &z);
1014 target = &s->constellation[constellation_state];
1015 /* Measure the training error */
1016 zz = complex_subf(&z, target);
1017 s->training_error += powerf(&zz);
1018 if (++s->training_count >= V17_TRAINING_SEG_4_LEN)
1019 {
1020 if (s->training_error < V17_TRAINING_SEG_4_LEN*constellation_spacing[s->space_map])
1021 {
1022 /* We are up and running */
1023 span_log(&s->logging, SPAN_LOG_FLOW, "Training succeeded at %dbps (constellation mismatch %f)\n", s->bit_rate, s->training_error);
1024 report_status_change(s, SIG_STATUS_TRAINING_SUCCEEDED);
1025 /* Apply some lag to the carrier off condition, to ensure the last few bits get pushed through
1026 the processing. */
1027 s->signal_present = 60;
1028 equalizer_save(s);
1029 s->carrier_phase_rate_save = s->carrier_phase_rate;
1030 s->short_train = TRUE;
1031 s->training_stage = TRAINING_STAGE_NORMAL_OPERATION;
1032 }
1033 else
1034 {
1035 /* Training has failed */
1036 span_log(&s->logging, SPAN_LOG_FLOW, "Training failed (constellation mismatch %f)\n", s->training_error);
1037 /* Park this modem */
1038 if (!s->short_train)
1039 s->agc_scaling_save = 0.0f;
1040 s->training_stage = TRAINING_STAGE_PARKED;
1041 report_status_change(s, SIG_STATUS_TRAINING_FAILED);
1042 }
1043 }
1044 break;
1045 case TRAINING_STAGE_PARKED:
1046 default:
1047 /* We failed to train! */
1048 /* Park here until the carrier drops. */
1049 target = &zero;
1050 break;
1051 }
1052 if (s->qam_report)
1053 s->qam_report(s->qam_user_data, &z, target, constellation_state);
1054 }
1055 /*- End of function --------------------------------------------------------*/
1056
1057 static __inline__ int signal_detect(v17_rx_state_t *s, int16_t amp)
1058 {
1059 int16_t diff;
1060 int16_t x;
1061 int32_t power;
1062
1063 /* There should be no DC in the signal, but sometimes there is.
1064 We need to measure the power with the DC blocked, but not using
1065 a slow to respond DC blocker. Use the most elementary HPF. */
1066 x = amp >> 1;
1067 /* There could be overflow here, but it isn't a problem in practice */
1068 diff = x - s->last_sample;
1069 s->last_sample = x;
1070 power = power_meter_update(&(s->power), diff);
1071 #if defined(IAXMODEM_STUFF)
1072 /* Quick power drop fudge */
1073 diff = abs(diff);
1074 if (10*diff < s->high_sample)
1075 {
1076 if (++s->low_samples > 120)
1077 {
1078 power_meter_init(&(s->power), 4);
1079 s->high_sample = 0;
1080 s->low_samples = 0;
1081 }
1082 }
1083 else
1084 {
1085 s->low_samples = 0;
1086 if (diff > s->high_sample)
1087 s->high_sample = diff;
1088 }
1089 #endif
1090 if (s->signal_present > 0)
1091 {
1092 /* Look for power below turn-off threshold to turn the carrier off */
1093 #if defined(IAXMODEM_STUFF)
1094 if (s->carrier_drop_pending || power < s->carrier_off_power)
1095 #else
1096 if (power < s->carrier_off_power)
1097 #endif
1098 {
1099 if (--s->signal_present <= 0)
1100 {
1101 /* Count down a short delay, to ensure we push the last
1102 few bits through the filters before stopping. */
1103 v17_rx_restart(s, s->bit_rate, s->short_train);
1104 report_status_change(s, SIG_STATUS_CARRIER_DOWN);
1105 return 0;
1106 }
1107 #if defined(IAXMODEM_STUFF)
1108 /* Carrier has dropped, but the put_bit is pending the signal_present delay. */
1109 s->carrier_drop_pending = TRUE;
1110 #endif
1111 }
1112 }
1113 else
1114 {
1115 /* Look for power exceeding turn-on threshold to turn the carrier on */
1116 if (power < s->carrier_on_power)
1117 return 0;
1118 s->signal_present = 1;
1119 #if defined(IAXMODEM_STUFF)
1120 s->carrier_drop_pending = FALSE;
1121 #endif
1122 report_status_change(s, SIG_STATUS_CARRIER_UP);
1123 }
1124 return power;
1125 }
1126 /*- End of function --------------------------------------------------------*/
1127
1128 SPAN_DECLARE_NONSTD(int) v17_rx(v17_rx_state_t *s, const int16_t amp[], int len)
1129 {
1130 int i;
1131 int step;
1132 complexf_t z;
1133 complexf_t zz;
1134 complexf_t sample;
1135 #if defined(SPANDSP_USE_FIXED_POINT)
1136 int32_t vi;
1137 #endif
1138 #if defined(SPANDSP_USE_FIXED_POINTx)
1139 int32_t v;
1140 #else
1141 float v;
1142 #endif
1143 int32_t power;
1144
1145 for (i = 0; i < len; i++)
1146 {
1147 s->rrc_filter[s->rrc_filter_step] = amp[i];
1148 if (++s->rrc_filter_step >= V17_RX_FILTER_STEPS)
1149 s->rrc_filter_step = 0;
1150
1151 if ((power = signal_detect(s, amp[i])) == 0)
1152 continue;
1153 if (s->training_stage == TRAINING_STAGE_PARKED)
1154 continue;
1155 /* Only spend effort processing this data if the modem is not
1156 parked, after training failure. */
1157 s->eq_put_step -= RX_PULSESHAPER_COEFF_SETS;
1158 step = -s->eq_put_step;
1159 if (step > RX_PULSESHAPER_COEFF_SETS - 1)
1160 step = RX_PULSESHAPER_COEFF_SETS - 1;
1161 if (step < 0)
1162 step += RX_PULSESHAPER_COEFF_SETS;
1163 #if defined(SPANDSP_USE_FIXED_POINT)
1164 vi = vec_circular_dot_prodi16(s->rrc_filter, rx_pulseshaper_re[step], V17_RX_FILTER_STEPS, s->rrc_filter_step);
1165 //sample.re = (vi*(int32_t) s->agc_scaling) >> 15;
1166 sample.re = vi*s->agc_scaling;
1167 #else
1168 v = vec_circular_dot_prodf(s->rrc_filter, rx_pulseshaper_re[step], V17_RX_FILTER_STEPS, s->rrc_filter_step);
1169 sample.re = v*s->agc_scaling;
1170 #endif
1171 /* Symbol timing synchronisation band edge filters */
1172 /* Low Nyquist band edge filter */
1173 v = s->symbol_sync_low[0]*SYNC_LOW_BAND_EDGE_COEFF_0 + s->symbol_sync_low[1]*SYNC_LOW_BAND_EDGE_COEFF_1 + sample.re;
1174 s->symbol_sync_low[1] = s->symbol_sync_low[0];
1175 s->symbol_sync_low[0] = v;
1176 /* High Nyquist band edge filter */
1177 v = s->symbol_sync_high[0]*SYNC_HIGH_BAND_EDGE_COEFF_0 + s->symbol_sync_high[1]*SYNC_HIGH_BAND_EDGE_COEFF_1 + sample.re;
1178 s->symbol_sync_high[1] = s->symbol_sync_high[0];
1179 s->symbol_sync_high[0] = v;
1180
1181 /* Put things into the equalization buffer at T/2 rate. The symbol sync.
1182 will fiddle the step to align this with the symbols. */
1183 if (s->eq_put_step <= 0)
1184 {
1185 /* Only AGC until we have locked down the setting. */
1186 if (s->agc_scaling_save == 0.0f)
1187 s->agc_scaling = (1.0f/RX_PULSESHAPER_GAIN)*2.17f/sqrtf(power);
1188 /* Pulse shape while still at the carrier frequency, using a quadrature
1189 pair of filters. This results in a properly bandpass filtered complex
1190 signal, which can be brought directly to baseband by complex mixing.
1191 No further filtering, to remove mixer harmonics, is needed. */
1192 step = -s->eq_put_step;
1193 if (step > RX_PULSESHAPER_COEFF_SETS - 1)
1194 step = RX_PULSESHAPER_COEFF_SETS - 1;
1195 s->eq_put_step += RX_PULSESHAPER_COEFF_SETS*10/(3*2);
1196 #if defined(SPANDSP_USE_FIXED_POINT)
1197 vi = vec_circular_dot_prodi16(s->rrc_filter, rx_pulseshaper_im[step], V17_RX_FILTER_STEPS, s->rrc_filter_step);
1198 //sample.im = (vi*(int32_t) s->agc_scaling) >> 15;
1199 sample.im = vi*s->agc_scaling;
1200 z = dds_lookup_complexf(s->carrier_phase);
1201 zz.re = sample.re*z.re - sample.im*z.im;
1202 zz.im = -sample.re*z.im - sample.im*z.re;
1203 #else
1204 v = vec_circular_dot_prodf(s->rrc_filter, rx_pulseshaper_im[step], V17_RX_FILTER_STEPS, s->rrc_filter_step);
1205 sample.im = v*s->agc_scaling;
1206 z = dds_lookup_complexf(s->carrier_phase);
1207 zz.re = sample.re*z.re - sample.im*z.im;
1208 zz.im = -sample.re*z.im - sample.im*z.re;
1209 #endif
1210 process_half_baud(s, &zz);
1211 }
1212 #if defined(SPANDSP_USE_FIXED_POINT)
1213 dds_advance(&s->carrier_phase, s->carrier_phase_rate);
1214 #else
1215 dds_advancef(&s->carrier_phase, s->carrier_phase_rate);
1216 #endif
1217 }
1218 return 0;
1219 }
1220 /*- End of function --------------------------------------------------------*/
1221
1222 SPAN_DECLARE(int) v17_rx_fillin(v17_rx_state_t *s, int len)
1223 {
1224 int i;
1225
1226 /* We want to sustain the current state (i.e carrier on<->carrier off), and
1227 try to sustain the carrier phase. We should probably push the filters, as well */
1228 span_log(&s->logging, SPAN_LOG_FLOW, "Fill-in %d samples\n", len);
1229 if (s->signal_present <= 0)
1230 return 0;
1231 if (s->training_stage == TRAINING_STAGE_PARKED)
1232 return 0;
1233 for (i = 0; i < len; i++)
1234 {
1235 #if defined(SPANDSP_USE_FIXED_POINT)
1236 dds_advance(&s->carrier_phase, s->carrier_phase_rate);
1237 #else
1238 dds_advancef(&s->carrier_phase, s->carrier_phase_rate);
1239 #endif
1240 /* Advance the symbol phase the appropriate amount */
1241 s->eq_put_step -= RX_PULSESHAPER_COEFF_SETS;
1242 if (s->eq_put_step <= 0)
1243 s->eq_put_step += RX_PULSESHAPER_COEFF_SETS*10/(3*2);
1244 /* TODO: Should we rotate any buffers */
1245 }
1246 return 0;
1247 }
1248 /*- End of function --------------------------------------------------------*/
1249
1250 SPAN_DECLARE(void) v17_rx_set_put_bit(v17_rx_state_t *s, put_bit_func_t put_bit, void *user_data)
1251 {
1252 s->put_bit = put_bit;
1253 s->put_bit_user_data = user_data;
1254 }
1255 /*- End of function --------------------------------------------------------*/
1256
1257 SPAN_DECLARE(void) v17_rx_set_modem_status_handler(v17_rx_state_t *s, modem_tx_status_func_t handler, void *user_data)
1258 {
1259 s->status_handler = handler;
1260 s->status_user_data = user_data;
1261 }
1262 /*- End of function --------------------------------------------------------*/
1263
1264 SPAN_DECLARE(logging_state_t *) v17_rx_get_logging_state(v17_rx_state_t *s)
1265 {
1266 return &s->logging;
1267 }
1268 /*- End of function --------------------------------------------------------*/
1269
1270 SPAN_DECLARE(int) v17_rx_restart(v17_rx_state_t *s, int bit_rate, int short_train)
1271 {
1272 int i;
1273
1274 span_log(&s->logging, SPAN_LOG_FLOW, "Restarting V.17, %dbps, %s training\n", bit_rate, (short_train) ? "short" : "long");
1275 switch (bit_rate)
1276 {
1277 case 14400:
1278 s->constellation = v17_v32bis_14400_constellation;
1279 s->space_map = 0;
1280 s->bits_per_symbol = 6;
1281 break;
1282 case 12000:
1283 s->constellation = v17_v32bis_12000_constellation;
1284 s->space_map = 1;
1285 s->bits_per_symbol = 5;
1286 break;
1287 case 9600:
1288 s->constellation = v17_v32bis_9600_constellation;
1289 s->space_map = 2;
1290 s->bits_per_symbol = 4;
1291 break;
1292 case 7200:
1293 s->constellation = v17_v32bis_7200_constellation;
1294 s->space_map = 3;
1295 s->bits_per_symbol = 3;
1296 break;
1297 case 4800:
1298 /* This does not exist in the V.17 spec as a valid mode of operation.
1299 However, it does exist in V.32bis, so it is here for completeness. */
1300 s->constellation = v17_v32bis_4800_constellation;
1301 s->space_map = 0;
1302 s->bits_per_symbol = 2;
1303 break;
1304 default:
1305 return -1;
1306 }
1307 s->bit_rate = bit_rate;
1308 #if defined(SPANDSP_USE_FIXED_POINT)
1309 vec_zeroi16(s->rrc_filter, sizeof(s->rrc_filter)/sizeof(s->rrc_filter[0]));
1310 #else
1311 vec_zerof(s->rrc_filter, sizeof(s->rrc_filter)/sizeof(s->rrc_filter[0]));
1312 #endif
1313 s->rrc_filter_step = 0;
1314
1315 s->diff = 1;
1316 s->scramble_reg = 0x2ECDD5;
1317 s->training_stage = TRAINING_STAGE_SYMBOL_ACQUISITION;
1318 s->training_count = 0;
1319 s->training_error = 0.0f;
1320 s->signal_present = 0;
1321 #if defined(IAXMODEM_STUFF)
1322 s->high_sample = 0;
1323 s->low_samples = 0;
1324 s->carrier_drop_pending = FALSE;
1325 #endif
1326 if (short_train != 2)
1327 s->short_train = short_train;
1328 memset(s->start_angles, 0, sizeof(s->start_angles));
1329 memset(s->angles, 0, sizeof(s->angles));
1330
1331 /* Initialise the TCM decoder parameters. */
1332 /* The accumulated distance vectors are set so state zero starts
1333 at a value of zero, and all others start larger. This forces the
1334 initial paths to merge at the zero states. */
1335 for (i = 0; i < 8; i++)
1336 #if defined(SPANDSP_USE_FIXED_POINTx)
1337 s->distances[i] = 99*DIST_FACTOR*DIST_FACTOR;
1338 #else
1339 s->distances[i] = 99.0f;
1340 #endif
1341 memset(s->full_path_to_past_state_locations, 0, sizeof(s->full_path_to_past_state_locations));
1342 memset(s->past_state_locations, 0, sizeof(s->past_state_locations));
1343 s->distances[0] = 0;
1344 s->trellis_ptr = 14;
1345
1346 span_log(&s->logging, SPAN_LOG_FLOW, "Phase rates %f %f\n", dds_frequencyf(s->carrier_phase_rate), dds_frequencyf(s->carrier_phase_rate_save));
1347 s->carrier_phase = 0;
1348 power_meter_init(&(s->power), 4);
1349
1350 if (s->short_train)
1351 {
1352 s->carrier_phase_rate = s->carrier_phase_rate_save;
1353 s->agc_scaling = s->agc_scaling_save;
1354 equalizer_restore(s);
1355 /* Don't allow any frequency correction at all, until we start to pull the phase in. */
1356 #if defined(SPANDSP_USE_FIXED_POINTx)
1357 s->carrier_track_i = 0;
1358 s->carrier_track_p = 40000;
1359 #else
1360 s->carrier_track_i = 0.0f;
1361 s->carrier_track_p = 40000.0f;
1362 #endif
1363 }
1364 else
1365 {
1366 s->carrier_phase_rate = dds_phase_ratef(CARRIER_NOMINAL_FREQ);
1367 equalizer_reset(s);
1368 #if defined(SPANDSP_USE_FIXED_POINTx)
1369 s->agc_scaling_save = 0;
1370 s->agc_scaling = (float) FP_FACTOR*32768.0f*0.0017f/RX_PULSESHAPER_GAIN;
1371 s->carrier_track_i = 5000;
1372 s->carrier_track_p = 40000;
1373 #else
1374 s->agc_scaling_save = 0.0f;
1375 s->agc_scaling = 0.0017f/RX_PULSESHAPER_GAIN;
1376 s->carrier_track_i = 5000.0f;
1377 s->carrier_track_p = 40000.0f;
1378 #endif
1379 }
1380 s->last_sample = 0;
1381
1382 /* Initialise the working data for symbol timing synchronisation */
1383 #if defined(SPANDSP_USE_FIXED_POINTx)
1384 for (i = 0; i < 2; i++)
1385 {
1386 s->symbol_sync_low[i] = 0;
1387 s->symbol_sync_high[i] = 0;
1388 s->symbol_sync_dc_filter[i] = 0;
1389 }
1390 s->baud_phase = 0;
1391 #else
1392 for (i = 0; i < 2; i++)
1393 {
1394 s->symbol_sync_low[i] = 0.0f;
1395 s->symbol_sync_high[i] = 0.0f;
1396 s->symbol_sync_dc_filter[i] = 0.0f;
1397 }
1398 s->baud_phase = 0.0f;
1399 #endif
1400 s->baud_half = 0;
1401
1402 s->total_baud_timing_correction = 0;
1403
1404 return 0;
1405 }
1406 /*- End of function --------------------------------------------------------*/
1407
1408 SPAN_DECLARE(v17_rx_state_t *) v17_rx_init(v17_rx_state_t *s, int bit_rate, put_bit_func_t put_bit, void *user_data)
1409 {
1410 switch (bit_rate)
1411 {
1412 case 14400:
1413 case 12000:
1414 case 9600:
1415 case 7200:
1416 case 4800:
1417 /* 4800 is an extension of V.17, to provide full converage of the V.32bis modes */
1418 break;
1419 default:
1420 return NULL;
1421 }
1422 if (s == NULL)
1423 {
1424 if ((s = (v17_rx_state_t *) malloc(sizeof(*s))) == NULL)
1425 return NULL;
1426 }
1427 memset(s, 0, sizeof(*s));
1428 span_log_init(&s->logging, SPAN_LOG_NONE, NULL);
1429 span_log_set_protocol(&s->logging, "V.17 RX");
1430 s->put_bit = put_bit;
1431 s->put_bit_user_data = user_data;
1432 s->short_train = FALSE;
1433 //s->scrambler_tap = 18 - 1;
1434 v17_rx_signal_cutoff(s, -45.5f);
1435 s->agc_scaling = 0.0017f/RX_PULSESHAPER_GAIN;
1436 s->agc_scaling_save = 0.0f;
1437 s->carrier_phase_rate_save = dds_phase_ratef(CARRIER_NOMINAL_FREQ);
1438 v17_rx_restart(s, bit_rate, s->short_train);
1439 return s;
1440 }
1441 /*- End of function --------------------------------------------------------*/
1442
1443 SPAN_DECLARE(int) v17_rx_release(v17_rx_state_t *s)
1444 {
1445 return 0;
1446 }
1447 /*- End of function --------------------------------------------------------*/
1448
1449 SPAN_DECLARE(int) v17_rx_free(v17_rx_state_t *s)
1450 {
1451 free(s);
1452 return 0;
1453 }
1454 /*- End of function --------------------------------------------------------*/
1455
1456 SPAN_DECLARE(void) v17_rx_set_qam_report_handler(v17_rx_state_t *s, qam_report_handler_t handler, void *user_data)
1457 {
1458 s->qam_report = handler;
1459 s->qam_user_data = user_data;
1460 }
1461 /*- End of function --------------------------------------------------------*/
1462 /*- End of file ------------------------------------------------------------*/

Repositories maintained by Peter Meerwald, pmeerw@pmeerw.net.