git.gag.com Git - debian/gnuradio/blob - gr-atsc/src/lib/atsci_equalizer_lms2.cc

   1 /* -*- c++ -*- */
   2 /*
   3  * Copyright 2002 Free Software Foundation, Inc.
   4  *
   5  * This file is part of GNU Radio
   6  *
   7  * GNU Radio is free software; you can redistribute it and/or modify
   8  * it under the terms of the GNU General Public License as published by
   9  * the Free Software Foundation; either version 3, or (at your option)
  10  * any later version.
  11  *
  12  * GNU Radio is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  * GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with GNU Radio; see the file COPYING.  If not, write to
  19  * the Free Software Foundation, Inc., 51 Franklin Street,
  20  * Boston, MA 02110-1301, USA.
  21  */
  22
  23 #include <atsci_equalizer_lms2.h>
  24 #include <assert.h>
  25 #include <algorithm>
  26 #include <atsci_pnXXX.h>
  27 #include <cmath>
  28 #include <stdlib.h>
  29 #include <gr_math.h>
  30 #include <stdio.h>
  31
  32 using std::min;
  33 using std::max;
  34
  35 static const int        NFFTAPS =  64;
  36 static const int        NFBTAPS = 192;
  37
  38 // the length of the field sync pattern that we know unequivocally
  39 static const int KNOWN_FIELD_SYNC_LENGTH = 4 + 511 + 3 * 63;
  40
  41 // static const float *get_data_seg_sync_training_sequence (int offset);
  42 static int          get_field_sync_training_sequence_length (int offset);
  43 static const float *get_field_sync_training_sequence (int which_field, int offset);
  44
  45 static inline int
  46 wrap (int d)
  47 {
  48   assert (d >= 0 && d <= (2 * NFBTAPS));
  49
  50   if(d >= NFBTAPS)
  51     return d - NFBTAPS;
  52   return d;
  53 }
  54
  55 static inline float
  56 slice (float d)
  57 {
  58   if (gr_isnan (d))
  59     return 0.0;
  60
  61   if (d >= 0.0){
  62     if (d >= 4.0){
  63       if (d >= 6.0)
  64         return 7.0;
  65       else
  66         return 5.0;
  67     }
  68     if (d >= 2.0)
  69       return 3.0;
  70     return 1.0;
  71   }
  72   else
  73     return -slice (-d);
  74 }
  75
  76 atsci_equalizer_lms2::atsci_equalizer_lms2 ()
  77   : d_taps_ff (NFFTAPS), d_taps_fb (NFBTAPS), d_old_output (NFBTAPS)
  78 {
  79   for (int i = 0; i < NFFTAPS; i++) {
  80     d_taps_ff[i] = 0.0;
  81   }
  82   for (int i = 0; i < NFBTAPS; i++) {
  83     d_taps_fb[i] = 0.0;
  84     d_old_output[i] = 0.0;
  85   }
  86   d_output_ptr = 0;
  87   trainingfile=fopen("taps.txt","w");
  88 }
  89
  90 atsci_equalizer_lms2::~atsci_equalizer_lms2 ()
  91 {
  92 }
  93
  94 void
  95 atsci_equalizer_lms2::reset ()
  96 {
  97   atsci_equalizer::reset ();            // invoke superclass
  98   for (int i = 0; i < NFFTAPS; i++) {
  99     d_taps_ff[i] = 0.0;
 100   }
 101   for (int i = 0; i < NFBTAPS; i++) {
 102     d_taps_fb[i] = 0.0;
 103     d_old_output[i] = 0.0;
 104   }
 105   d_output_ptr = 0;
 106 }
 107
 108
 109 int
 110 atsci_equalizer_lms2::ntaps () const
 111 {
 112   return NFFTAPS + NFBTAPS;
 113 }
 114
 115 int
 116 atsci_equalizer_lms2::npretaps () const
 117 {
 118   return NFFTAPS;
 119 }
 120
 121 /*!
 122  * Input range is known NOT TO CONTAIN data segment syncs
 123  * or field syncs.  This should be the fast path.  In the
 124  * non decicion directed case, this just runs the input
 125  * through the filter without adapting it.
 126  *
 127  * \p input_samples has (nsamples + ntaps() - 1) valid entries.
 128  * input_samples[0] .. input_samples[nsamples - 1 + ntaps() - 1] may be
 129  * referenced to compute the output values.
 130  */
 131 void
 132 atsci_equalizer_lms2::filter_normal (const float *input_samples,
 133                                    float *output_samples,
 134                                    int   nsamples)
 135 {
 136   // handle data
 137   filterN (input_samples, output_samples, nsamples);
 138 }
 139
 140
 141 /*!
 142  * Input range is known to consist of only a data segment sync or a
 143  * portion of a data segment sync.  \p nsamples will be in [1,4].
 144  * \p offset will be in [0,3].  \p offset is the offset of the input
 145  * from the beginning of the data segment sync pattern.
 146  *
 147  * \p input_samples has (nsamples + ntaps() - 1) valid entries.
 148  * input_samples[0] .. input_samples[nsamples - 1 + ntaps() - 1] may be
 149  * referenced to compute the output values.
 150  */
 151 void
 152 atsci_equalizer_lms2::filter_data_seg_sync (const float *input_samples,
 153                                           float *output_samples,
 154                                           int   nsamples,
 155                                           int   offset)
 156 {
 157   // handle data
 158   //  adaptN (input_samples, get_data_seg_sync_training_sequence (offset),
 159   //      output_samples, nsamples);
 160  filterN (input_samples, output_samples, nsamples);
 161
 162  //  cerr << "Seg Sync: offset " << offset << "\tnsamples\t" << nsamples << "\t pre, 5 -5 -5 5\t" <<
 163  // output_samples[0] << "\t" << output_samples[1] << "\t" << output_samples[2] << "\t" << output_samples[3] << endl;
 164
 165 }
 166
 167
 168 /*!
 169  * Input range is known to consist of only a field sync segment or a
 170  * portion of a field sync segment.  \p nsamples will be in [1,832].
 171  * \p offset will be in [0,831].  \p offset is the offset of the input
 172  * from the beginning of the data segment sync pattern.  We consider the
 173  * 4 symbols of the immediately preceding data segment sync to be the
 174  * first symbols of the field sync segment.  \p which_field is in [0,1]
 175  * and specifies which field (duh).
 176  *
 177  * \p input_samples has (nsamples + ntaps() - 1) valid entries.
 178  * input_samples[0] .. input_samples[nsamples - 1 + ntaps() - 1] may be
 179  * referenced to compute the output values.
 180  */
 181 void
 182 atsci_equalizer_lms2::filter_field_sync (const float *input_samples,
 183                                        float *output_samples,
 184                                        int   nsamples,
 185                                        int   offset,
 186                                        int   which_field)
 187 {
 188   // Only the first 4 + 511 + 3 * 63 symbols are completely defined.
 189   // Those after that the symbols are bilevel, so we could use decision feedback and use
 190   // that to train, but for now, don't train on them.
 191
 192   int   n = min (nsamples, get_field_sync_training_sequence_length (offset));
 193
 194   // handle known training sequence
 195   adaptN (input_samples, get_field_sync_training_sequence (which_field, offset),
 196                     output_samples, n);
 197
 198   // just filter any unknown portion
 199   if (nsamples > n)
 200     filterN (&input_samples[n], &output_samples[n], nsamples - n);
 201
 202   if (offset == 0 && nsamples > 0){
 203     for (int i = 0; i < NFFTAPS; i++)
 204       fprintf(trainingfile,"%f ",d_taps_ff[i]);
 205     for (int i = 0; i < NFBTAPS; i++)
 206       fprintf(trainingfile,"%f ",d_taps_fb[i]);
 207     fprintf (trainingfile,"\n");
 208   }
 209
 210 }
 211
 212 // ----------------------------------------------------------------
 213
 214 //
 215 // filter a single output
 216 //
 217 float
 218 atsci_equalizer_lms2::filter1 (const float input[])
 219 {
 220   static const int N_UNROLL = 4;
 221
 222   float acc0 = 0;
 223   float acc1 = 0;
 224   float acc2 = 0;
 225   float acc3 = 0;
 226   float acc = 0;
 227
 228
 229   unsigned      i = 0;
 230   unsigned      n = (NFFTAPS / N_UNROLL) * N_UNROLL;
 231
 232   for (i = 0; i < n; i += N_UNROLL){
 233     acc0 += d_taps_ff[i + 0] * input[i + 0];
 234     acc1 += d_taps_ff[i + 1] * input[i + 1];
 235     acc2 += d_taps_ff[i + 2] * input[i + 2];
 236     acc3 += d_taps_ff[i + 3] * input[i + 3];
 237   }
 238
 239   for (; i < (unsigned) NFFTAPS; i++)
 240     acc0 += d_taps_ff[i] * input[i];
 241
 242   acc = (acc0 + acc1 + acc2 + acc3);
 243
 244   d_output_ptr = wrap (d_output_ptr + 1);
 245
 246   for (int i = 0; i < NFBTAPS; i++) {
 247     acc -= d_taps_fb[i] * d_old_output[wrap(i + d_output_ptr)];
 248   }
 249
 250   if (gr_isnan (acc)){
 251     abort ();
 252   }
 253
 254   d_old_output[d_output_ptr] = slice (acc);
 255   return acc;
 256 }
 257
 258 //
 259 // filter and adapt a single output
 260 //
 261 float kludge ()
 262 {
 263   return 0.0;
 264 }
 265
 266 float
 267 atsci_equalizer_lms2::adapt1 (const float input[], float ideal_output)
 268 {
 269   static const double BETA = 0.00005;   // FIXME figure out what this ought to be
 270                                         // FIXME add gear-shifting
 271
 272   double y = filter1 (input);
 273   double e = y - ideal_output;
 274
 275   // update taps...
 276   for (int i = 0; i < NFFTAPS; i++){
 277     d_taps_ff[i] = d_taps_ff[i] - BETA * e * (double)(input[i]);
 278   }
 279
 280   for (int i = 0; i < NFBTAPS; i++){
 281     // d_taps_fb[i] = d_taps_fb[i] - BETA * e * (double)d_old_output[wrap(i+d_output_ptr)];
 282     d_taps_fb[i] = d_taps_fb[i] - kludge() * e * (double)d_old_output[wrap(i+d_output_ptr)];
 283   }
 284
 285   return y;
 286 }
 287
 288 void
 289 atsci_equalizer_lms2::filterN (const float *input_samples,
 290                              float *output_samples,
 291                              int nsamples)
 292 {
 293   for (int i = 0; i < nsamples; i++)
 294     output_samples[i] = filter1 (&input_samples[i]);
 295 }
 296
 297
 298 void
 299 atsci_equalizer_lms2::adaptN (const float *input_samples,
 300                             const float *training_pattern,
 301                             float *output_samples,
 302                             int    nsamples)
 303 {
 304   for (int i = 0; i < nsamples; i++)
 305     output_samples[i] = adapt1 (&input_samples[i], training_pattern[i]);
 306 }
 307
 308 // ----------------------------------------------------------------
 309
 310 static float
 311 bin_map (int bit)
 312 {
 313   return bit ? +5 : -5;
 314 }
 315
 316 static void
 317 init_field_sync_common (float *p, int mask)
 318
 319 {
 320   int  i = 0;
 321
 322   p[i++] = bin_map (1);                 // data segment sync pulse
 323   p[i++] = bin_map (0);
 324   p[i++] = bin_map (0);
 325   p[i++] = bin_map (1);
 326
 327   for (int j = 0; j < 511; j++)         // PN511
 328     p[i++] = bin_map (atsc_pn511[j]);
 329
 330   for (int j = 0; j < 63; j++)          // PN63
 331     p[i++] = bin_map (atsc_pn63[j]);
 332
 333   for (int j = 0; j < 63; j++)          // PN63, toggled on field 2
 334     p[i++] = bin_map (atsc_pn63[j] ^ mask);
 335
 336   for (int j = 0; j < 63; j++)          // PN63
 337     p[i++] = bin_map (atsc_pn63[j]);
 338
 339   assert (i == KNOWN_FIELD_SYNC_LENGTH);
 340 }
 341
 342 #if 0
 343 static const float *
 344 get_data_seg_sync_training_sequence (int offset)
 345 {
 346   static const float training_data[4] = { +5, -5, -5, +5 };
 347   return &training_data[offset];
 348 }
 349 #endif
 350
 351 static int
 352 get_field_sync_training_sequence_length (int offset)
 353 {
 354   return max (0, KNOWN_FIELD_SYNC_LENGTH - offset);
 355 }
 356
 357 static const float *
 358 get_field_sync_training_sequence (int which_field, int offset)
 359 {
 360   static float *field_1 = 0;
 361   static float *field_2 = 0;
 362
 363   if (field_1 == 0){
 364     field_1 = new float[KNOWN_FIELD_SYNC_LENGTH];
 365     field_2 = new float[KNOWN_FIELD_SYNC_LENGTH];
 366     init_field_sync_common (field_1, 0);
 367     init_field_sync_common (field_2, 1);
 368   }
 369
 370   if (which_field == 0)
 371     return &field_1[offset];
 372   else
 373     return &field_2[offset];
 374 }