git.gag.com Git - debian/gnuradio/blob - gr-atsc/src/lib/atsci_equalizer_lms.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_lms.h>
  24 #include <assert.h>
  25 #include <algorithm>
  26 #include <atsci_pnXXX.h>
  27
  28 #include <stdio.h>
  29
  30 using std::min;
  31 using std::max;
  32
  33 static const int        NTAPS = 256;
  34 static const int        NPRETAPS = (int) (NTAPS * 0.8); // probably should be either .2 or .8
  35
  36
  37 // the length of the field sync pattern that we know unequivocally
  38 static const int KNOWN_FIELD_SYNC_LENGTH = 4 + 511 + 3 * 63;
  39
  40 static const float *get_data_seg_sync_training_sequence (int offset);
  41 static int          get_field_sync_training_sequence_length (int offset);
  42 static const float *get_field_sync_training_sequence (int which_field, int offset);
  43
  44
  45 atsci_equalizer_lms::atsci_equalizer_lms () : d_taps (NTAPS)
  46 {
  47   for (int i = 0; i < NTAPS; i++) {
  48     d_taps[i] = 0.0;
  49   }
  50   trainingfile=fopen("taps.txt","w");
  51 }
  52
  53 atsci_equalizer_lms::~atsci_equalizer_lms ()
  54 {
  55 }
  56
  57 void
  58 atsci_equalizer_lms::reset ()
  59 {
  60   atsci_equalizer::reset ();            // invoke superclass
  61
  62   for (int i = 0; i < NTAPS; i++) {
  63     d_taps[i] = 0.0;
  64   }
  65 }
  66
  67 int
  68 atsci_equalizer_lms::ntaps () const
  69 {
  70   return NTAPS;
  71 }
  72
  73 int
  74 atsci_equalizer_lms::npretaps () const
  75 {
  76   return NPRETAPS;
  77 }
  78
  79 /*!
  80  * Input range is known NOT TO CONTAIN data segment syncs
  81  * or field syncs.  This should be the fast path.  In the
  82  * non decicion directed case, this just runs the input
  83  * through the filter without adapting it.
  84  *
  85  * \p input_samples has (nsamples + ntaps() - 1) valid entries.
  86  * input_samples[0] .. input_samples[nsamples - 1 + ntaps() - 1] may be
  87  * referenced to compute the output values.
  88  */
  89 void
  90 atsci_equalizer_lms::filter_normal (const float *input_samples,
  91                                    float *output_samples,
  92                                    int   nsamples)
  93 {
  94   // handle data
  95   filterN (input_samples, output_samples, nsamples);
  96 }
  97
  98
  99 /*!
 100  * Input range is known to consist of only a data segment sync or a
 101  * portion of a data segment sync.  \p nsamples will be in [1,4].
 102  * \p offset will be in [0,3].  \p offset is the offset of the input
 103  * from the beginning of the data segment sync pattern.
 104  *
 105  * \p input_samples has (nsamples + ntaps() - 1) valid entries.
 106  * input_samples[0] .. input_samples[nsamples - 1 + ntaps() - 1] may be
 107  * referenced to compute the output values.
 108  */
 109 void
 110 atsci_equalizer_lms::filter_data_seg_sync (const float *input_samples,
 111                                           float *output_samples,
 112                                           int   nsamples,
 113                                           int   offset)
 114 {
 115   // handle data
 116   //  adaptN (input_samples, get_data_seg_sync_training_sequence (offset),
 117   //      output_samples, nsamples);
 118  filterN (input_samples, output_samples, nsamples);
 119
 120  //  cerr << "Seg Sync: offset " << offset << "\tnsamples\t" << nsamples << "\t pre, 5 -5 -5 5\t" <<
 121  // output_samples[0] << "\t" << output_samples[1] << "\t" << output_samples[2] << "\t" << output_samples[3] << endl;
 122
 123 }
 124
 125
 126 /*!
 127  * Input range is known to consist of only a field sync segment or a
 128  * portion of a field sync segment.  \p nsamples will be in [1,832].
 129  * \p offset will be in [0,831].  \p offset is the offset of the input
 130  * from the beginning of the data segment sync pattern.  We consider the
 131  * 4 symbols of the immediately preceding data segment sync to be the
 132  * first symbols of the field sync segment.  \p which_field is in [0,1]
 133  * and specifies which field (duh).
 134  *
 135  * \p input_samples has (nsamples + ntaps() - 1) valid entries.
 136  * input_samples[0] .. input_samples[nsamples - 1 + ntaps() - 1] may be
 137  * referenced to compute the output values.
 138  */
 139 void
 140 atsci_equalizer_lms::filter_field_sync (const float *input_samples,
 141                                        float *output_samples,
 142                                        int   nsamples,
 143                                        int   offset,
 144                                        int   which_field)
 145 {
 146   // Only the first 4 + 511 + 3 * 63 symbols are completely defined.
 147   // Those after that the symbols are bilevel, so we could use decision feedback and use
 148   // that to train, but for now, don't train on them.
 149
 150   int   n = min (nsamples, get_field_sync_training_sequence_length (offset));
 151
 152   // handle known training sequence
 153   adaptN (input_samples, get_field_sync_training_sequence (which_field, offset),
 154                     output_samples, n);
 155
 156   // just filter any unknown portion
 157   if (nsamples > n)
 158     filterN (&input_samples[n], &output_samples[n], nsamples - n);
 159
 160   if (offset == 0 && nsamples > 0){
 161     for (int i = 0; i < NTAPS; i++)
 162       fprintf(trainingfile,"%f ",d_taps[i]);
 163
 164     fprintf (trainingfile,"\n");
 165   }
 166
 167 }
 168
 169 // ----------------------------------------------------------------
 170
 171 //
 172 // filter a single output
 173 //
 174 float
 175 atsci_equalizer_lms::filter1 (const float input[])
 176 {
 177   static const int N_UNROLL = 4;
 178
 179   float acc0 = 0;
 180   float acc1 = 0;
 181   float acc2 = 0;
 182   float acc3 = 0;
 183
 184
 185   unsigned      i = 0;
 186   unsigned      n = (NTAPS / N_UNROLL) * N_UNROLL;
 187
 188   for (i = 0; i < n; i += N_UNROLL){
 189     acc0 += d_taps[i + 0] * input[i + 0];
 190     acc1 += d_taps[i + 1] * input[i + 1];
 191     acc2 += d_taps[i + 2] * input[i + 2];
 192     acc3 += d_taps[i + 3] * input[i + 3];
 193   }
 194
 195   for (; i < (unsigned) NTAPS; i++)
 196     acc0 += d_taps[i] * input[i];
 197
 198   return (acc0 + acc1 + acc2 + acc3);
 199 }
 200
 201 //
 202 // filter and adapt a single output
 203 //
 204 float
 205 atsci_equalizer_lms::adapt1 (const float input[], float ideal_output)
 206 {
 207   static const double BETA = 0.00005;   // FIXME figure out what this ought to be
 208                                         // FIXME add gear-shifting
 209
 210   double y = filter1 (input);
 211   double e = y - ideal_output;
 212
 213   // update taps...
 214   for (int i = 0; i < NTAPS; i++){
 215     d_taps[i] = d_taps[i] - BETA * e * (double)(input[i]);
 216   }
 217
 218   return y;
 219 }
 220
 221 void
 222 atsci_equalizer_lms::filterN (const float *input_samples,
 223                              float *output_samples,
 224                              int nsamples)
 225 {
 226   for (int i = 0; i < nsamples; i++)
 227     output_samples[i] = filter1 (&input_samples[i]);
 228 }
 229
 230
 231 void
 232 atsci_equalizer_lms::adaptN (const float *input_samples,
 233                             const float *training_pattern,
 234                             float *output_samples,
 235                             int    nsamples)
 236 {
 237   for (int i = 0; i < nsamples; i++)
 238     output_samples[i] = adapt1 (&input_samples[i], training_pattern[i]);
 239 }
 240
 241 // ----------------------------------------------------------------
 242
 243 static float
 244 bin_map (int bit)
 245 {
 246   return bit ? +5 : -5;
 247 }
 248
 249 static void
 250 init_field_sync_common (float *p, int mask)
 251
 252 {
 253   int  i = 0;
 254
 255   p[i++] = bin_map (1);                 // data segment sync pulse
 256   p[i++] = bin_map (0);
 257   p[i++] = bin_map (0);
 258   p[i++] = bin_map (1);
 259
 260   for (int j = 0; j < 511; j++)         // PN511
 261     p[i++] = bin_map (atsc_pn511[j]);
 262
 263   for (int j = 0; j < 63; j++)          // PN63
 264     p[i++] = bin_map (atsc_pn63[j]);
 265
 266   for (int j = 0; j < 63; j++)          // PN63, toggled on field 2
 267     p[i++] = bin_map (atsc_pn63[j] ^ mask);
 268
 269   for (int j = 0; j < 63; j++)          // PN63
 270     p[i++] = bin_map (atsc_pn63[j]);
 271
 272   assert (i == KNOWN_FIELD_SYNC_LENGTH);
 273 }
 274
 275
 276 static const float *
 277 get_data_seg_sync_training_sequence (int offset)
 278 {
 279   static const float training_data[4] = { +5, -5, -5, +5 };
 280   return &training_data[offset];
 281 }
 282
 283 static int
 284 get_field_sync_training_sequence_length (int offset)
 285 {
 286   return max (0, KNOWN_FIELD_SYNC_LENGTH - offset);
 287 }
 288
 289 static const float *
 290 get_field_sync_training_sequence (int which_field, int offset)
 291 {
 292   static float *field_1 = 0;
 293   static float *field_2 = 0;
 294
 295   if (field_1 == 0){
 296     field_1 = new float[KNOWN_FIELD_SYNC_LENGTH];
 297     field_2 = new float[KNOWN_FIELD_SYNC_LENGTH];
 298     init_field_sync_common (field_1, 0);
 299     init_field_sync_common (field_2, 1);
 300   }
 301
 302   if (which_field == 0)
 303     return &field_1[offset];
 304   else
 305     return &field_2[offset];
 306 }