3 * Copyright 2002 Free Software Foundation, Inc.
5 * This file is part of GNU Radio
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)
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.
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.
23 #include <atsci_equalizer_lms.h>
26 #include <atsci_pnXXX.h>
33 static const int NTAPS = 256;
34 static const int NPRETAPS = (int) (NTAPS * 0.8); // probably should be either .2 or .8
37 // the length of the field sync pattern that we know unequivocally
38 static const int KNOWN_FIELD_SYNC_LENGTH = 4 + 511 + 3 * 63;
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);
45 atsci_equalizer_lms::atsci_equalizer_lms () : d_taps (NTAPS)
47 for (int i = 0; i < NTAPS; i++) {
50 trainingfile=fopen("taps.txt","w");
53 atsci_equalizer_lms::~atsci_equalizer_lms ()
58 atsci_equalizer_lms::reset ()
60 atsci_equalizer::reset (); // invoke superclass
62 for (int i = 0; i < NTAPS; i++) {
68 atsci_equalizer_lms::ntaps () const
74 atsci_equalizer_lms::npretaps () const
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.
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.
90 atsci_equalizer_lms::filter_normal (const float *input_samples,
91 float *output_samples,
95 filterN (input_samples, output_samples, nsamples);
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.
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.
110 atsci_equalizer_lms::filter_data_seg_sync (const float *input_samples,
111 float *output_samples,
116 // adaptN (input_samples, get_data_seg_sync_training_sequence (offset),
117 // output_samples, nsamples);
118 filterN (input_samples, output_samples, nsamples);
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;
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).
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.
140 atsci_equalizer_lms::filter_field_sync (const float *input_samples,
141 float *output_samples,
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.
150 int n = min (nsamples, get_field_sync_training_sequence_length (offset));
152 // handle known training sequence
153 adaptN (input_samples, get_field_sync_training_sequence (which_field, offset),
156 // just filter any unknown portion
158 filterN (&input_samples[n], &output_samples[n], nsamples - n);
160 if (offset == 0 && nsamples > 0){
161 for (int i = 0; i < NTAPS; i++)
162 fprintf(trainingfile,"%f ",d_taps[i]);
164 fprintf (trainingfile,"\n");
169 // ----------------------------------------------------------------
172 // filter a single output
175 atsci_equalizer_lms::filter1 (const float input[])
177 static const int N_UNROLL = 4;
186 unsigned n = (NTAPS / N_UNROLL) * N_UNROLL;
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];
195 for (; i < (unsigned) NTAPS; i++)
196 acc0 += d_taps[i] * input[i];
198 return (acc0 + acc1 + acc2 + acc3);
202 // filter and adapt a single output
205 atsci_equalizer_lms::adapt1 (const float input[], float ideal_output)
207 static const double BETA = 0.00005; // FIXME figure out what this ought to be
208 // FIXME add gear-shifting
210 double y = filter1 (input);
211 double e = y - ideal_output;
214 for (int i = 0; i < NTAPS; i++){
215 d_taps[i] = d_taps[i] - BETA * e * (double)(input[i]);
222 atsci_equalizer_lms::filterN (const float *input_samples,
223 float *output_samples,
226 for (int i = 0; i < nsamples; i++)
227 output_samples[i] = filter1 (&input_samples[i]);
232 atsci_equalizer_lms::adaptN (const float *input_samples,
233 const float *training_pattern,
234 float *output_samples,
237 for (int i = 0; i < nsamples; i++)
238 output_samples[i] = adapt1 (&input_samples[i], training_pattern[i]);
241 // ----------------------------------------------------------------
246 return bit ? +5 : -5;
250 init_field_sync_common (float *p, int mask)
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);
260 for (int j = 0; j < 511; j++) // PN511
261 p[i++] = bin_map (atsc_pn511[j]);
263 for (int j = 0; j < 63; j++) // PN63
264 p[i++] = bin_map (atsc_pn63[j]);
266 for (int j = 0; j < 63; j++) // PN63, toggled on field 2
267 p[i++] = bin_map (atsc_pn63[j] ^ mask);
269 for (int j = 0; j < 63; j++) // PN63
270 p[i++] = bin_map (atsc_pn63[j]);
272 assert (i == KNOWN_FIELD_SYNC_LENGTH);
277 get_data_seg_sync_training_sequence (int offset)
279 static const float training_data[4] = { +5, -5, -5, +5 };
280 return &training_data[offset];
285 get_field_sync_training_sequence_length (int offset)
287 return max (0, KNOWN_FIELD_SYNC_LENGTH - offset);
291 get_field_sync_training_sequence (int which_field, int offset)
293 static float *field_1 = 0;
294 static float *field_2 = 0;
297 field_1 = new float[KNOWN_FIELD_SYNC_LENGTH];
298 field_2 = new float[KNOWN_FIELD_SYNC_LENGTH];
299 init_field_sync_common (field_1, 0);
300 init_field_sync_common (field_2, 1);
303 if (which_field == 0)
304 return &field_1[offset];
306 return &field_2[offset];