final/MultiSource/Benchmarks/VersaBench/beamformer/beamformer.c - external/llvm.org/test-suite - Git at Google

 /*
  * Copyright (c) 2003 David Maze
  *
  * Permission  is hereby  granted,  free  of  charge, to  any  person
  * obtaining a  copy of  this software  and  associated documentation
  * files   (the  "Software"),  to   deal  in  the   Software  without
  * restriction, including without  limitation the rights to use, copy,
  * modify, merge, publish,  distribute, sublicense, and/or sell copies
  * of  the Software,  and to  permit persons to  whom the  Software is
  * furnished to do so, subject to the following conditions:
  *
  * The above  copyright notice  and this  permission notice  shall be
  * included in all copies or substantial portions of the Software.
  *
  * THE  SOFTWARE IS  PROVIDED "AS IS",  WITHOUT WARRANTY OF  ANY KIND,
  * EXPRESS OR IMPLIED, INCLUDING  BUT NOT LIMITED TO THE WARRANTIES OF
  * MERCHANTABILITY,   FITNESS   FOR   A   PARTICULAR    PURPOSE    AND
  * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
  * BE LIABLE FOR ANY CLAIM,  DAMAGES OR OTHER LIABILITY, WHETHER IN AN
  * ACTION OF  CONTRACT, TORT OR OTHERWISE, ARISING FROM,  OUT OF OR IN
  * CONNECTION  WITH THE SOFTWARE OR  THE USE OR OTHER  DEALINGS IN THE
  * SOFTWARE.
  */

 /*
  * beamformer.c: Standalone beam-former reference implementation
  * David Maze <dmaze@cag.lcs.mit.edu>
  * $Id: beamformer.c,v 1.5 2003/11/07 08:47:00 thies Exp $
  */

 /* Modified by: Rodric M. Rabbah 06-03-04 */

 #include <stdio.h>
 #include <stdlib.h>
 #include <unistd.h>
 #include <math.h>

 /*
 This implementation is derived from the StreamIt implementation,
 rather than from the PCA VSIPL-based implementation.  It is intended
 to be easier to synchronize our reference implementation with our
 StreamIt implementation, not have dependencies on extra libraries, and
 generally be a fairer comparison with the StreamIt code (that is,
 equivalent to our other benchmarks).

 This version gets consistent output with the StreamIt version when
 compiled with RANDOM_INPUTS and RANDOM_WEIGHTS.  The algorithm should
 be equivalent to the SERIALIZED implementation.

 FLAGS: the StreamIt TemplateBeamFormer can be built with COARSE
 defined or not, with SERIALIZED defined or not, and with RANDOM_INPUTS
 defined or not.  Defining COARSE changes the algorithm to add two
 decimation stages; this presently isn't implemented.  SERIALIZED
 changes where in the graph printing happens; in the non-SERIALIZED
 version printing happens in the detector stage, but this doesn't have
 a deterministic order.  RANDOM_INPUTS and RANDOM_WEIGHTS change
 weights at several points in the code, where "RANDOM" is "something
 the author made up" rather than "determined via rand()".  The default
 flags are both RANDOM_WEIGHTS and RANDOM_INPUTS.

 RANDOM_INPUTS and RANDOM_WEIGHTS are easy to do, so they are
 implemented here as well.  We ignore SERIALIZED, though we do go
 through the detectors in a deterministic order which should give the
 right answer, and add a reordering stage equivalent to what SERIALZIED
 adds to get results in the same order.  COARSE is left for a future
 implementation.
 */

 #define RANDOM_INPUTS

 struct BeamFirData
 {
   int len, count, pos;
   float *weight, *buffer;
 };

 void begin(void);
 void InputGenerate(int channel, float *inputs, int n);
 void BeamFirSetup(struct BeamFirData *data, int n);
 void BeamFirFilter(struct BeamFirData *data,
                    int input_length, int decimation_ratio,
                    float *in, float *out);
 void BeamFormWeights(int beam, float *weights);
 void BeamForm(int beam, const float *weights, const float *input,
               float *output);
 void Magnitude(float *in, float *out, int n);
 void Detector(int beam, float *data, float *output);

 void begin_StrictFP(void);
 void InputGenerate_StrictFP(int channel, float *inputs, int n);
 void BeamFirSetup_StrictFP(struct BeamFirData *data, int n);
 void BeamFirFilter_StrictFP(struct BeamFirData *data,
                             int input_length, int decimation_ratio,
                             float *in, float *out);
 void BeamFormWeights_StrictFP(int beam, float *weights);
 void BeamForm_StrictFP(int beam, const float *weights, const float *input,
                        float *output);
 void Magnitude_StrictFP(float *in, float *out, int n);
 void Detector_StrictFP(int beam, float *data, float *output);

 #define NUM_CHANNELS 12
 #define NUM_SAMPLES 1024
 #define NUM_BEAMS 4
 /* Bill writes: under current implementation, decimation ratios must
    divide NUM_SAMPLES. */
 #define NUM_COARSE_TAPS 64
 #define NUM_FINE_TAPS 64
 #define COARSE_DECIMATION_RATIO 1
 #define FINE_DECIMATION_RATIO 2
 #define NUM_SEGMENTS 1
 #define NUM_POST_DEC_1 (NUM_SAMPLES/COARSE_DECIMATION_RATIO)
 #define NUM_POST_DEC_2 (NUM_POST_DEC_1/FINE_DECIMATION_RATIO)
 #define MF_SIZE (NUM_SEGMENTS*NUM_POST_DEC_2)
 #define PULSE_SIZE (NUM_POST_DEC_2/2)
 #define PREDEC_PULSE_SIZE \
   (PULSE_SIZE*COARSE_DECIMATION_RATIO*FINE_DECIMATION_RATIO)
 #define TARGET_BEAM (NUM_BEAMS/4)
 #define TARGET_SAMPLE (NUM_SAMPLES/4)
 #define TARGET_SAMPLE_POST_DEC \
    (TARGET_SAMPLE/COARSE_DECIMATION_RATIO/FINE_DECIMATION_RATIO)
 #define D_OVER_LAMBDA 0.5
 #define CFAR_THRESHOLD (0.95*D_OVER_LAMBDA*NUM_CHANNELS*0.5*PULSE_SIZE)

 float detector_out_StrictFP[NUM_BEAMS][NUM_POST_DEC_2];
 static int numiters = 100;

 int main(int argc, char **argv)
 {
   int option;

   while ((option = getopt(argc, argv, "i:")) != -1)
   {
     switch(option)
     {
     case 'i':
       numiters = atoi(optarg);
     }
   }

   begin_StrictFP();
   begin();
   return 0;
 }

 /* Return 0 when V1 and V2 do not match within the allowed FP_ABSTOLERANCE.  */
 static inline int
 check_FP(float V1, float V2) {
   double AbsTolerance = FP_ABSTOLERANCE;
   double Diff = fabs(V1 - V2);
   if (Diff > AbsTolerance) {
     fprintf(stderr, "A = %lf and B = %lf differ more than"
                   " FP_ABSTOLERANCE = %lf\n", V1, V2, AbsTolerance);
     return 0;
   }

   return 1;
 }

 /* What are the counts here?  Deriving:
  *
  * Detector takes TARGET_SAMPLE_POST_DEC inputs.
  * Magnitude takes twice as many.
  * The various BeamFirFilters have internal circular buffers,
  *   so their counts are only affected by the decimation rates.
  *   The last BeamFirFilter has no decimation, so its input count
  *   is also 2*TARGET_SAMPLE_POST_DEC.
  * BeamForm takes NUM_CHANNELS times as many.
  *
  * So, coming out of the first for loop, the buffer should have
  * 2*TARGET_SAMPLE_POST_DEC*NUM_CHANNELS values.  This comes from
  * a roundrobin(2) joiner in StreamIt, for NUM_CHANNELS children.
  *
  * The second BeamFirFilter has 2*TARGET_SAMPLE_POST_DEC outputs.
  * The first BeamFirFilter has
  *   2*TARGET_SAMPLE_POST_DEC*FINE_DECIMATION_RATIO outputs.
  * The input generator produces 2*TARGET_SAMPLE values.
  */

 void begin(void)
 {
 #ifdef AVOID_PRINTF
   volatile float result;
 #endif
   struct BeamFirData coarse_fir_data[NUM_CHANNELS];
   struct BeamFirData fine_fir_data[NUM_CHANNELS];
   struct BeamFirData mf_fir_data[NUM_BEAMS];
   float inputs[2*NUM_SAMPLES*NUM_CHANNELS];
   float predec[2*NUM_POST_DEC_1*NUM_CHANNELS];
   float postdec[NUM_CHANNELS][2*NUM_POST_DEC_2*NUM_CHANNELS];
   float beam_weights[NUM_BEAMS][2*NUM_CHANNELS];
   float beam_input[2*NUM_CHANNELS*NUM_POST_DEC_2];
   float beam_output[2*NUM_POST_DEC_2];
   float beam_fir_output[2*NUM_POST_DEC_2];
   float beam_fir_mag[NUM_POST_DEC_2];
   float detector_out[NUM_BEAMS][NUM_POST_DEC_2];
   int i, j;

   for (i = 0; i < NUM_CHANNELS; i++)
   {
     BeamFirSetup(&coarse_fir_data[i], NUM_COARSE_TAPS);
     BeamFirSetup(&fine_fir_data[i], NUM_FINE_TAPS);
   }
   for (i = 0; i < NUM_BEAMS; i++)
   {
     BeamFirSetup(&mf_fir_data[i], MF_SIZE);
     BeamFormWeights(i, beam_weights[i]);
   }

   /*** VERSABENCH START ***/

   while (numiters == -1 || numiters-- > 0) {
     for (i = 0; i < NUM_CHANNELS; i++) {
       for (j = 0; j < NUM_CHANNELS; j++)
         InputGenerate(i, inputs + j*NUM_SAMPLES*2,
                       NUM_SAMPLES);
       for (j = 0; j < NUM_POST_DEC_1; j++)
         BeamFirFilter(&coarse_fir_data[i],
                       NUM_SAMPLES, COARSE_DECIMATION_RATIO,
                       inputs + j * COARSE_DECIMATION_RATIO*2,
                       predec + j * 2);
       for (j = 0; j < NUM_POST_DEC_2; j++)
         BeamFirFilter(&fine_fir_data[i],
                       NUM_POST_DEC_1, FINE_DECIMATION_RATIO,
                       predec + j * FINE_DECIMATION_RATIO * 2,
                       postdec[i] + j * 2);
     }
     /* Assemble beam-forming input: */
     for (i = 0; i < NUM_CHANNELS; i++)
       for (j = 0; j < NUM_POST_DEC_2; j++)
       {
         beam_input[j*NUM_CHANNELS*2+2*i] = postdec[i][2*j];
         beam_input[j*NUM_CHANNELS*2+2*i+1] = postdec[i][2*j+1];
       }
     for (i = 0; i < NUM_BEAMS; i++)
     {
       /* Have now rearranged NUM_CHANNELS*NUM_POST_DEC_2 items.
        * BeamForm takes NUM_CHANNELS inputs, 2 outputs. */
       for (j = 0; j < NUM_POST_DEC_2; j++)
         BeamForm(i, beam_weights[i],
                  beam_input + j*NUM_CHANNELS*2,
                  beam_output + j*2);
       for (j = 0; j < NUM_POST_DEC_2; j++)
         BeamFirFilter(&mf_fir_data[i],
                       NUM_POST_DEC_2, 1,
                       beam_output+j*2, beam_fir_output+j*2);
       Magnitude(beam_fir_output, beam_fir_mag, NUM_POST_DEC_2);
       Detector(i, beam_fir_mag, detector_out[i]);
     }
     for (j = 0; j < NUM_POST_DEC_2; j++)
       for (i = 0; i < NUM_BEAMS; i++) {
         if (!check_FP(detector_out[i][j], detector_out_StrictFP[i][j]))
           exit(1);
 #ifdef AVOID_PRINTF
         result = detector_out[i][j];
 #else
         printf("%f\n", detector_out_StrictFP[i][j]);
 #endif
       }
   }

   /*** VERSABENCH END ***/

 }

 void InputGenerate(int channel, float *inputs, int n)
 {
   int i;
   for (i = 0; i < n; i++)
   {
     if (channel == TARGET_BEAM && i == TARGET_SAMPLE)
     {
 #ifdef RANDOM_INPUTS
       inputs[2*i] = sqrt(i*channel);
       inputs[2*i+1] = sqrt(i*channel)+1;
 #else
       inputs[2*i] = sqrt(CFAR_THRESHOLD);
       inputs[2*i+1] = 0;
 #endif
     } else {
 #ifdef RANDOM_INPUTS
       float root = sqrt(i*channel);
       inputs[2*i] = -root;
       inputs[2*i+1] = -(root+1);
 #else
       inputs[2*i] = 0;
       inputs[2*i+1] = 0;
 #endif
     }
   }
 }

 void BeamFirSetup(struct BeamFirData *data, int n)
 {
   int i, j;

   data->len = n;
   data->count = 0;
   data->pos = 0;
   data->weight = malloc(sizeof(float)*2*n);
   data->buffer = malloc(sizeof(float)*2*n);

 #ifdef RANDOM_WEIGHTS
   for (j = 0; j < n; j++) {
     int idx = j+1;
     data->weight[j*2] = sin(idx) / ((float)idx);
     data->weight[j*2+1] = cos(idx) / ((float)idx);
     data->buffer[j*2] = 0.0;
     data->buffer[j*2+1] = 0.0;
   }
 #else
   data->weight[0] = 1.0;
   data->weight[1] = 0.0;
   for (i = 1; i < 2*n; i++) {
     data->weight[i] = 0.0;
     data->buffer[i] = 0.0;
   }
 #endif
 }

 void BeamFirFilter(struct BeamFirData *data,
                    int input_length, int decimation_ratio,
                    float *in, float *out)
 {
   /* Input must be exactly 2*decimation_ratio long; output must be
    * exactly 2 long. */
   float real_curr = 0;
   float imag_curr = 0;
   int i;
   int modPos;
   int len, mask, mask2;

   len = data->len;
   mask = len - 1;
   mask2 = 2 * len - 1;
   modPos = 2*(len - 1 - data->pos);
   data->buffer[modPos] = in[0];
   data->buffer[modPos+1] = in[1];

   /* Profiling says: this is the single inner loop that matters! */
   for (i = 0; i < 2*len; i+=2) {
     float rd = data->buffer[modPos];
     float id = data->buffer[modPos+1];
     float rw = data->weight[i];
     float iw = data->weight[i+1];
     float rci = rd * rw + id * iw;
     /* sign error?  this is consistent with StreamIt --dzm */
     float ici = id * rw + rd * iw;
     real_curr += rci;
     imag_curr += ici;
     modPos = (modPos + 2) & mask2;
   }

   data->pos = (data->pos + 1) & mask;
   out[0] = real_curr;
   out[1] = imag_curr;
   data->count += decimation_ratio;
   if (data->count == input_length)
   {
     data->count = 0;
     data->pos = 0;
     for (i = 0; i < 2*data->len; i++) {
       data->buffer[i] = 0.0;
     }
   }
 }

 void BeamFormWeights(int beam, float *weights)
 {
   int i;
   for (i = 0; i < NUM_CHANNELS; i++)
   {
 #ifdef RANDOM_WEIGHTS
     int idx = i+1;
     weights[2*i] = sin(idx) / (float)(beam+idx);
     weights[2*i+1] = cos(idx) / (float)(beam+idx);
 #else
     if (i == beam) {
       weights[2*i] = 1;
       weights[2*i+1] = 0;
     } else {
       weights[2*i] = 0;
       weights[2*i+1] = 0;
     }
 #endif
   }
 }

 void BeamForm(int beam, const float *weights, const float *input,
               float *output)
 {
   /* 2*NUM_CHANNELS inputs and weights; 2 outputs. */
   float real_curr = 0;
   float imag_curr = 0;
   int i;
   for (i = 0; i < NUM_CHANNELS; i++)
   {
     real_curr += weights[2*i] * input[2*i] - weights[2*i+1] * input[2*i+1];
     imag_curr += weights[2*i] * input[2*i+1] + weights[2*i+1] * input[2*i];
   }
   output[0] = real_curr;
   output[1] = imag_curr;
 }

 void Magnitude(float *in, float *out, int n)
 {
   int i;
   /* Should be 2n inputs, n outputs. */
   for (i = 0; i < n; i++)
     out[i] = sqrt(in[2*i]*in[2*i] + in[2*i+1]*in[2*i+1]);
 }

 void Detector(int beam, float *data, float *output)
 {
   int sample;
   /* Should be exactly NUM_POST_DEC_2 samples. */
   for (sample = 0; sample < NUM_POST_DEC_2; sample++)
   {
     float outputVal;
     if (beam == TARGET_BEAM && sample == TARGET_SAMPLE) {
       if (data[sample] >= 0.1)
         outputVal = beam+1;
       else
         outputVal = 0;
     } else {
       if (data[sample] >= 0.1)
         outputVal = -(beam+1);
       else
         outputVal = 0;
     }
     outputVal = data[sample];
     output[sample]= outputVal;
   }
 }

 void begin_StrictFP(void)
 {
 #pragma STDC FP_CONTRACT OFF
   struct BeamFirData coarse_fir_data[NUM_CHANNELS];
   struct BeamFirData fine_fir_data[NUM_CHANNELS];
   struct BeamFirData mf_fir_data[NUM_BEAMS];
   float inputs[2*NUM_SAMPLES*NUM_CHANNELS];
   float predec[2*NUM_POST_DEC_1*NUM_CHANNELS];
   float postdec[NUM_CHANNELS][2*NUM_POST_DEC_2*NUM_CHANNELS];
   float beam_weights[NUM_BEAMS][2*NUM_CHANNELS];
   float beam_input[2*NUM_CHANNELS*NUM_POST_DEC_2];
   float beam_output[2*NUM_POST_DEC_2];
   float beam_fir_output[2*NUM_POST_DEC_2];
   float beam_fir_mag[NUM_POST_DEC_2];
   int i, j;

   for (i = 0; i < NUM_CHANNELS; i++)
   {
     BeamFirSetup_StrictFP(&coarse_fir_data[i], NUM_COARSE_TAPS);
     BeamFirSetup_StrictFP(&fine_fir_data[i], NUM_FINE_TAPS);
   }
   for (i = 0; i < NUM_BEAMS; i++)
   {
     BeamFirSetup_StrictFP(&mf_fir_data[i], MF_SIZE);
     BeamFormWeights_StrictFP(i, beam_weights[i]);
   }

   for (i = 0; i < NUM_CHANNELS; i++) {
     for (j = 0; j < NUM_CHANNELS; j++)
       InputGenerate_StrictFP(i, inputs + j*NUM_SAMPLES*2,
                              NUM_SAMPLES);
     for (j = 0; j < NUM_POST_DEC_1; j++)
       BeamFirFilter_StrictFP(&coarse_fir_data[i],
                              NUM_SAMPLES, COARSE_DECIMATION_RATIO,
                              inputs + j * COARSE_DECIMATION_RATIO*2,
                              predec + j * 2);
     for (j = 0; j < NUM_POST_DEC_2; j++)
       BeamFirFilter_StrictFP(&fine_fir_data[i],
                              NUM_POST_DEC_1, FINE_DECIMATION_RATIO,
                              predec + j * FINE_DECIMATION_RATIO * 2,
                              postdec[i] + j * 2);
   }
   /* Assemble beam-forming input: */
   for (i = 0; i < NUM_CHANNELS; i++)
     for (j = 0; j < NUM_POST_DEC_2; j++)
       {
         beam_input[j*NUM_CHANNELS*2+2*i] = postdec[i][2*j];
         beam_input[j*NUM_CHANNELS*2+2*i+1] = postdec[i][2*j+1];
       }
   for (i = 0; i < NUM_BEAMS; i++)
     {
       /* Have now rearranged NUM_CHANNELS*NUM_POST_DEC_2 items.
        * BeamForm takes NUM_CHANNELS inputs, 2 outputs. */
       for (j = 0; j < NUM_POST_DEC_2; j++)
         BeamForm_StrictFP(i, beam_weights[i],
                           beam_input + j*NUM_CHANNELS*2,
                           beam_output + j*2);
       for (j = 0; j < NUM_POST_DEC_2; j++)
         BeamFirFilter_StrictFP(&mf_fir_data[i],
                                NUM_POST_DEC_2, 1,
                                beam_output+j*2, beam_fir_output+j*2);
       Magnitude_StrictFP(beam_fir_output, beam_fir_mag, NUM_POST_DEC_2);
       Detector_StrictFP(i, beam_fir_mag, detector_out_StrictFP[i]);
     }
 }

 void InputGenerate_StrictFP(int channel, float *inputs, int n)
 {
 #pragma STDC FP_CONTRACT OFF
   int i;
   for (i = 0; i < n; i++)
   {
     if (channel == TARGET_BEAM && i == TARGET_SAMPLE)
     {
 #ifdef RANDOM_INPUTS
       inputs[2*i] = sqrt(i*channel);
       inputs[2*i+1] = sqrt(i*channel)+1;
 #else
       inputs[2*i] = sqrt(CFAR_THRESHOLD);
       inputs[2*i+1] = 0;
 #endif
     } else {
 #ifdef RANDOM_INPUTS
       float root = sqrt(i*channel);
       inputs[2*i] = -root;
       inputs[2*i+1] = -(root+1);
 #else
       inputs[2*i] = 0;
       inputs[2*i+1] = 0;
 #endif
     }
   }
 }

 void BeamFirSetup_StrictFP(struct BeamFirData *data, int n)
 {
 #pragma STDC FP_CONTRACT OFF
   int i, j;

   data->len = n;
   data->count = 0;
   data->pos = 0;
   data->weight = malloc(sizeof(float)*2*n);
   data->buffer = malloc(sizeof(float)*2*n);

 #ifdef RANDOM_WEIGHTS
   for (j = 0; j < n; j++) {
     int idx = j+1;
     data->weight[j*2] = sin(idx) / ((float)idx);
     data->weight[j*2+1] = cos(idx) / ((float)idx);
     data->buffer[j*2] = 0.0;
     data->buffer[j*2+1] = 0.0;
   }
 #else
   data->weight[0] = 1.0;
   data->weight[1] = 0.0;
   for (i = 1; i < 2*n; i++) {
     data->weight[i] = 0.0;
     data->buffer[i] = 0.0;
   }
 #endif
 }

 void BeamFirFilter_StrictFP(struct BeamFirData *data,
                             int input_length, int decimation_ratio,
                             float *in, float *out)
 {
 #pragma STDC FP_CONTRACT OFF
   /* Input must be exactly 2*decimation_ratio long; output must be
    * exactly 2 long. */
   float real_curr = 0;
   float imag_curr = 0;
   int i;
   int modPos;
   int len, mask, mask2;

   len = data->len;
   mask = len - 1;
   mask2 = 2 * len - 1;
   modPos = 2*(len - 1 - data->pos);
   data->buffer[modPos] = in[0];
   data->buffer[modPos+1] = in[1];

   /* Profiling says: this is the single inner loop that matters! */
   for (i = 0; i < 2*len; i+=2) {
     float rd = data->buffer[modPos];
     float id = data->buffer[modPos+1];
     float rw = data->weight[i];
     float iw = data->weight[i+1];
     float rci = rd * rw + id * iw;
     /* sign error?  this is consistent with StreamIt --dzm */
     float ici = id * rw + rd * iw;
     real_curr += rci;
     imag_curr += ici;
     modPos = (modPos + 2) & mask2;
   }

   data->pos = (data->pos + 1) & mask;
   out[0] = real_curr;
   out[1] = imag_curr;
   data->count += decimation_ratio;
   if (data->count == input_length)
   {
     data->count = 0;
     data->pos = 0;
     for (i = 0; i < 2*data->len; i++) {
       data->buffer[i] = 0.0;
     }
   }
 }

 void BeamFormWeights_StrictFP(int beam, float *weights)
 {
 #pragma STDC FP_CONTRACT OFF
   int i;
   for (i = 0; i < NUM_CHANNELS; i++)
   {
 #ifdef RANDOM_WEIGHTS
     int idx = i+1;
     weights[2*i] = sin(idx) / (float)(beam+idx);
     weights[2*i+1] = cos(idx) / (float)(beam+idx);
 #else
     if (i == beam) {
       weights[2*i] = 1;
       weights[2*i+1] = 0;
     } else {
       weights[2*i] = 0;
       weights[2*i+1] = 0;
     }
 #endif
   }
 }

 void BeamForm_StrictFP(int beam, const float *weights, const float *input,
                        float *output)
 {
 #pragma STDC FP_CONTRACT OFF
   /* 2*NUM_CHANNELS inputs and weights; 2 outputs. */
   float real_curr = 0;
   float imag_curr = 0;
   int i;
   for (i = 0; i < NUM_CHANNELS; i++)
   {
     real_curr += weights[2*i] * input[2*i] - weights[2*i+1] * input[2*i+1];
     imag_curr += weights[2*i] * input[2*i+1] + weights[2*i+1] * input[2*i];
   }
   output[0] = real_curr;
   output[1] = imag_curr;
 }

 void Magnitude_StrictFP(float *in, float *out, int n)
 {
 #pragma STDC FP_CONTRACT OFF
   int i;
   /* Should be 2n inputs, n outputs. */
   for (i = 0; i < n; i++)
     out[i] = sqrt(in[2*i]*in[2*i] + in[2*i+1]*in[2*i+1]);
 }

 void Detector_StrictFP(int beam, float *data, float *output)
 {
 #pragma STDC FP_CONTRACT OFF
   int sample;
   /* Should be exactly NUM_POST_DEC_2 samples. */
   for (sample = 0; sample < NUM_POST_DEC_2; sample++)
   {
     float outputVal;
     if (beam == TARGET_BEAM && sample == TARGET_SAMPLE) {
       if (data[sample] >= 0.1)
         outputVal = beam+1;
       else
         outputVal = 0;
     } else {
       if (data[sample] >= 0.1)
         outputVal = -(beam+1);
       else
         outputVal = 0;
     }
     outputVal = data[sample];
     output[sample]= outputVal;
   }
 }
	/*
	* Copyright (c) 2003 David Maze
	*
	* Permission is hereby granted, free of charge, to any person
	* obtaining a copy of this software and associated documentation
	* files (the "Software"), to deal in the Software without
	* restriction, including without limitation the rights to use, copy,
	* modify, merge, publish, distribute, sublicense, and/or sell copies
	* of the Software, and to permit persons to whom the Software is
	* furnished to do so, subject to the following conditions:
	*
	* The above copyright notice and this permission notice shall be
	* included in all copies or substantial portions of the Software.
	*
	* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
	* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
	* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
	* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
	* BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
	* ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
	* CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
	* SOFTWARE.
	*/

	/*
	* beamformer.c: Standalone beam-former reference implementation
	* David Maze <dmaze@cag.lcs.mit.edu>
	* $Id: beamformer.c,v 1.5 2003/11/07 08:47:00 thies Exp $
	*/

	/* Modified by: Rodric M. Rabbah 06-03-04 */

	#include <stdio.h>
	#include <stdlib.h>
	#include <unistd.h>
	#include <math.h>

	/*
	This implementation is derived from the StreamIt implementation,
	rather than from the PCA VSIPL-based implementation. It is intended
	to be easier to synchronize our reference implementation with our
	StreamIt implementation, not have dependencies on extra libraries, and
	generally be a fairer comparison with the StreamIt code (that is,
	equivalent to our other benchmarks).

	This version gets consistent output with the StreamIt version when
	compiled with RANDOM_INPUTS and RANDOM_WEIGHTS. The algorithm should
	be equivalent to the SERIALIZED implementation.

	FLAGS: the StreamIt TemplateBeamFormer can be built with COARSE
	defined or not, with SERIALIZED defined or not, and with RANDOM_INPUTS
	defined or not. Defining COARSE changes the algorithm to add two
	decimation stages; this presently isn't implemented. SERIALIZED
	changes where in the graph printing happens; in the non-SERIALIZED
	version printing happens in the detector stage, but this doesn't have
	a deterministic order. RANDOM_INPUTS and RANDOM_WEIGHTS change
	weights at several points in the code, where "RANDOM" is "something
	the author made up" rather than "determined via rand()". The default
	flags are both RANDOM_WEIGHTS and RANDOM_INPUTS.

	RANDOM_INPUTS and RANDOM_WEIGHTS are easy to do, so they are
	implemented here as well. We ignore SERIALIZED, though we do go
	through the detectors in a deterministic order which should give the
	right answer, and add a reordering stage equivalent to what SERIALZIED
	adds to get results in the same order. COARSE is left for a future
	implementation.
	*/

	#define RANDOM_INPUTS

	struct BeamFirData
	{
	int len, count, pos;
	float weight, buffer;
	};

	void begin(void);
	void InputGenerate(int channel, float *inputs, int n);
	void BeamFirSetup(struct BeamFirData *data, int n);
	void BeamFirFilter(struct BeamFirData *data,
	int input_length, int decimation_ratio,
	float in, float out);
	void BeamFormWeights(int beam, float *weights);
	void BeamForm(int beam, const float weights, const float input,
	float *output);
	void Magnitude(float in, float out, int n);
	void Detector(int beam, float data, float output);

	void begin_StrictFP(void);
	void InputGenerate_StrictFP(int channel, float *inputs, int n);
	void BeamFirSetup_StrictFP(struct BeamFirData *data, int n);
	void BeamFirFilter_StrictFP(struct BeamFirData *data,
	int input_length, int decimation_ratio,
	float in, float out);
	void BeamFormWeights_StrictFP(int beam, float *weights);
	void BeamForm_StrictFP(int beam, const float weights, const float input,
	float *output);
	void Magnitude_StrictFP(float in, float out, int n);
	void Detector_StrictFP(int beam, float data, float output);

	#define NUM_CHANNELS 12
	#define NUM_SAMPLES 1024
	#define NUM_BEAMS 4
	/* Bill writes: under current implementation, decimation ratios must
	divide NUM_SAMPLES. */
	#define NUM_COARSE_TAPS 64
	#define NUM_FINE_TAPS 64
	#define COARSE_DECIMATION_RATIO 1
	#define FINE_DECIMATION_RATIO 2
	#define NUM_SEGMENTS 1
	#define NUM_POST_DEC_1 (NUM_SAMPLES/COARSE_DECIMATION_RATIO)
	#define NUM_POST_DEC_2 (NUM_POST_DEC_1/FINE_DECIMATION_RATIO)
	#define MF_SIZE (NUM_SEGMENTS*NUM_POST_DEC_2)
	#define PULSE_SIZE (NUM_POST_DEC_2/2)
	#define PREDEC_PULSE_SIZE \
	(PULSE_SIZECOARSE_DECIMATION_RATIOFINE_DECIMATION_RATIO)
	#define TARGET_BEAM (NUM_BEAMS/4)
	#define TARGET_SAMPLE (NUM_SAMPLES/4)
	#define TARGET_SAMPLE_POST_DEC \
	(TARGET_SAMPLE/COARSE_DECIMATION_RATIO/FINE_DECIMATION_RATIO)
	#define D_OVER_LAMBDA 0.5
	#define CFAR_THRESHOLD (0.95D_OVER_LAMBDANUM_CHANNELS0.5PULSE_SIZE)

	float detector_out_StrictFP[NUM_BEAMS][NUM_POST_DEC_2];
	static int numiters = 100;

	int main(int argc, char **argv)
	{
	int option;

	while ((option = getopt(argc, argv, "i:")) != -1)
	{
	switch(option)
	{
	case 'i':
	numiters = atoi(optarg);
	}
	}

	begin_StrictFP();
	begin();
	return 0;
	}

	/* Return 0 when V1 and V2 do not match within the allowed FP_ABSTOLERANCE. */
	static inline int
	check_FP(float V1, float V2) {
	double AbsTolerance = FP_ABSTOLERANCE;
	double Diff = fabs(V1 - V2);
	if (Diff > AbsTolerance) {
	fprintf(stderr, "A = %lf and B = %lf differ more than"
	" FP_ABSTOLERANCE = %lf\n", V1, V2, AbsTolerance);
	return 0;
	}

	return 1;
	}

	/* What are the counts here? Deriving:
	*
	* Detector takes TARGET_SAMPLE_POST_DEC inputs.
	* Magnitude takes twice as many.
	* The various BeamFirFilters have internal circular buffers,
	* so their counts are only affected by the decimation rates.
	* The last BeamFirFilter has no decimation, so its input count
	* is also 2*TARGET_SAMPLE_POST_DEC.
	* BeamForm takes NUM_CHANNELS times as many.
	*
	* So, coming out of the first for loop, the buffer should have
	* 2TARGET_SAMPLE_POST_DECNUM_CHANNELS values. This comes from
	* a roundrobin(2) joiner in StreamIt, for NUM_CHANNELS children.
	*
	* The second BeamFirFilter has 2*TARGET_SAMPLE_POST_DEC outputs.
	* The first BeamFirFilter has
	* 2TARGET_SAMPLE_POST_DECFINE_DECIMATION_RATIO outputs.
	* The input generator produces 2*TARGET_SAMPLE values.
	*/

	void begin(void)
	{
	#ifdef AVOID_PRINTF
	volatile float result;
	#endif
	struct BeamFirData coarse_fir_data[NUM_CHANNELS];
	struct BeamFirData fine_fir_data[NUM_CHANNELS];
	struct BeamFirData mf_fir_data[NUM_BEAMS];
	float inputs[2NUM_SAMPLESNUM_CHANNELS];
	float predec[2NUM_POST_DEC_1NUM_CHANNELS];
	float postdec[NUM_CHANNELS][2NUM_POST_DEC_2NUM_CHANNELS];
	float beam_weights[NUM_BEAMS][2*NUM_CHANNELS];
	float beam_input[2NUM_CHANNELSNUM_POST_DEC_2];
	float beam_output[2*NUM_POST_DEC_2];
	float beam_fir_output[2*NUM_POST_DEC_2];
	float beam_fir_mag[NUM_POST_DEC_2];
	float detector_out[NUM_BEAMS][NUM_POST_DEC_2];
	int i, j;

	for (i = 0; i < NUM_CHANNELS; i++)
	{
	BeamFirSetup(&coarse_fir_data[i], NUM_COARSE_TAPS);
	BeamFirSetup(&fine_fir_data[i], NUM_FINE_TAPS);
	}
	for (i = 0; i < NUM_BEAMS; i++)
	{
	BeamFirSetup(&mf_fir_data[i], MF_SIZE);
	BeamFormWeights(i, beam_weights[i]);
	}

	/* VERSABENCH START */

	while (numiters == -1 \|\| numiters-- > 0) {
	for (i = 0; i < NUM_CHANNELS; i++) {
	for (j = 0; j < NUM_CHANNELS; j++)
	InputGenerate(i, inputs + jNUM_SAMPLES2,
	NUM_SAMPLES);
	for (j = 0; j < NUM_POST_DEC_1; j++)
	BeamFirFilter(&coarse_fir_data[i],
	NUM_SAMPLES, COARSE_DECIMATION_RATIO,
	inputs + j * COARSE_DECIMATION_RATIO*2,
	predec + j * 2);
	for (j = 0; j < NUM_POST_DEC_2; j++)
	BeamFirFilter(&fine_fir_data[i],
	NUM_POST_DEC_1, FINE_DECIMATION_RATIO,
	predec + j * FINE_DECIMATION_RATIO * 2,
	postdec[i] + j * 2);
	}
	/* Assemble beam-forming input: */
	for (i = 0; i < NUM_CHANNELS; i++)
	for (j = 0; j < NUM_POST_DEC_2; j++)
	{
	beam_input[jNUM_CHANNELS2+2i] = postdec[i][2j];
	beam_input[jNUM_CHANNELS2+2i+1] = postdec[i][2j+1];
	}
	for (i = 0; i < NUM_BEAMS; i++)
	{
	/* Have now rearranged NUM_CHANNELS*NUM_POST_DEC_2 items.
	* BeamForm takes NUM_CHANNELS inputs, 2 outputs. */
	for (j = 0; j < NUM_POST_DEC_2; j++)
	BeamForm(i, beam_weights[i],
	beam_input + jNUM_CHANNELS2,
	beam_output + j*2);
	for (j = 0; j < NUM_POST_DEC_2; j++)
	BeamFirFilter(&mf_fir_data[i],
	NUM_POST_DEC_2, 1,
	beam_output+j2, beam_fir_output+j2);
	Magnitude(beam_fir_output, beam_fir_mag, NUM_POST_DEC_2);
	Detector(i, beam_fir_mag, detector_out[i]);
	}
	for (j = 0; j < NUM_POST_DEC_2; j++)
	for (i = 0; i < NUM_BEAMS; i++) {
	if (!check_FP(detector_out[i][j], detector_out_StrictFP[i][j]))
	exit(1);
	#ifdef AVOID_PRINTF
	result = detector_out[i][j];
	#else
	printf("%f\n", detector_out_StrictFP[i][j]);
	#endif
	}
	}

	/* VERSABENCH END */

	}

	void InputGenerate(int channel, float *inputs, int n)
	{
	int i;
	for (i = 0; i < n; i++)
	{
	if (channel == TARGET_BEAM && i == TARGET_SAMPLE)
	{
	#ifdef RANDOM_INPUTS
	inputs[2i] = sqrt(ichannel);
	inputs[2i+1] = sqrt(ichannel)+1;
	#else
	inputs[2*i] = sqrt(CFAR_THRESHOLD);
	inputs[2*i+1] = 0;
	#endif
	} else {
	#ifdef RANDOM_INPUTS
	float root = sqrt(i*channel);
	inputs[2*i] = -root;
	inputs[2*i+1] = -(root+1);
	#else
	inputs[2*i] = 0;
	inputs[2*i+1] = 0;
	#endif
	}
	}
	}

	void BeamFirSetup(struct BeamFirData *data, int n)
	{
	int i, j;

	data->len = n;
	data->count = 0;
	data->pos = 0;
	data->weight = malloc(sizeof(float)2n);
	data->buffer = malloc(sizeof(float)2n);

	#ifdef RANDOM_WEIGHTS
	for (j = 0; j < n; j++) {
	int idx = j+1;
	data->weight[j*2] = sin(idx) / ((float)idx);
	data->weight[j*2+1] = cos(idx) / ((float)idx);
	data->buffer[j*2] = 0.0;
	data->buffer[j*2+1] = 0.0;
	}
	#else
	data->weight[0] = 1.0;
	data->weight[1] = 0.0;
	for (i = 1; i < 2*n; i++) {
	data->weight[i] = 0.0;
	data->buffer[i] = 0.0;
	}
	#endif
	}

	void BeamFirFilter(struct BeamFirData *data,
	int input_length, int decimation_ratio,
	float in, float out)
	{
	/* Input must be exactly 2*decimation_ratio long; output must be
	* exactly 2 long. */
	float real_curr = 0;
	float imag_curr = 0;
	int i;
	int modPos;
	int len, mask, mask2;

	len = data->len;
	mask = len - 1;
	mask2 = 2 * len - 1;
	modPos = 2*(len - 1 - data->pos);
	data->buffer[modPos] = in[0];
	data->buffer[modPos+1] = in[1];

	/* Profiling says: this is the single inner loop that matters! */
	for (i = 0; i < 2*len; i+=2) {
	float rd = data->buffer[modPos];
	float id = data->buffer[modPos+1];
	float rw = data->weight[i];
	float iw = data->weight[i+1];
	float rci = rd * rw + id * iw;
	/* sign error? this is consistent with StreamIt --dzm */
	float ici = id * rw + rd * iw;
	real_curr += rci;
	imag_curr += ici;
	modPos = (modPos + 2) & mask2;
	}

	data->pos = (data->pos + 1) & mask;
	out[0] = real_curr;
	out[1] = imag_curr;
	data->count += decimation_ratio;
	if (data->count == input_length)
	{
	data->count = 0;
	data->pos = 0;
	for (i = 0; i < 2*data->len; i++) {
	data->buffer[i] = 0.0;
	}
	}
	}

	void BeamFormWeights(int beam, float *weights)
	{
	int i;
	for (i = 0; i < NUM_CHANNELS; i++)
	{
	#ifdef RANDOM_WEIGHTS
	int idx = i+1;
	weights[2*i] = sin(idx) / (float)(beam+idx);
	weights[2*i+1] = cos(idx) / (float)(beam+idx);
	#else
	if (i == beam) {
	weights[2*i] = 1;
	weights[2*i+1] = 0;
	} else {
	weights[2*i] = 0;
	weights[2*i+1] = 0;
	}
	#endif
	}
	}

	void BeamForm(int beam, const float weights, const float input,
	float *output)
	{
	/* 2NUM_CHANNELS inputs and weights; 2 outputs. /
	float real_curr = 0;
	float imag_curr = 0;
	int i;
	for (i = 0; i < NUM_CHANNELS; i++)
	{
	real_curr += weights[2i] input[2i] - weights[2i+1] * input[2*i+1];
	imag_curr += weights[2i] input[2i+1] + weights[2i+1] * input[2*i];
	}
	output[0] = real_curr;
	output[1] = imag_curr;
	}

	void Magnitude(float in, float out, int n)
	{
	int i;
	/* Should be 2n inputs, n outputs. */
	for (i = 0; i < n; i++)
	out[i] = sqrt(in[2i]in[2i] + in[2i+1]in[2i+1]);
	}

	void Detector(int beam, float data, float output)
	{
	int sample;
	/* Should be exactly NUM_POST_DEC_2 samples. */
	for (sample = 0; sample < NUM_POST_DEC_2; sample++)
	{
	float outputVal;
	if (beam == TARGET_BEAM && sample == TARGET_SAMPLE) {
	if (data[sample] >= 0.1)
	outputVal = beam+1;
	else
	outputVal = 0;
	} else {
	if (data[sample] >= 0.1)
	outputVal = -(beam+1);
	else
	outputVal = 0;
	}
	outputVal = data[sample];
	output[sample]= outputVal;
	}
	}

	void begin_StrictFP(void)
	{
	#pragma STDC FP_CONTRACT OFF
	struct BeamFirData coarse_fir_data[NUM_CHANNELS];
	struct BeamFirData fine_fir_data[NUM_CHANNELS];
	struct BeamFirData mf_fir_data[NUM_BEAMS];
	float inputs[2NUM_SAMPLESNUM_CHANNELS];
	float predec[2NUM_POST_DEC_1NUM_CHANNELS];
	float postdec[NUM_CHANNELS][2NUM_POST_DEC_2NUM_CHANNELS];
	float beam_weights[NUM_BEAMS][2*NUM_CHANNELS];
	float beam_input[2NUM_CHANNELSNUM_POST_DEC_2];
	float beam_output[2*NUM_POST_DEC_2];
	float beam_fir_output[2*NUM_POST_DEC_2];
	float beam_fir_mag[NUM_POST_DEC_2];
	int i, j;

	for (i = 0; i < NUM_CHANNELS; i++)
	{
	BeamFirSetup_StrictFP(&coarse_fir_data[i], NUM_COARSE_TAPS);
	BeamFirSetup_StrictFP(&fine_fir_data[i], NUM_FINE_TAPS);
	}
	for (i = 0; i < NUM_BEAMS; i++)
	{
	BeamFirSetup_StrictFP(&mf_fir_data[i], MF_SIZE);
	BeamFormWeights_StrictFP(i, beam_weights[i]);
	}

	for (i = 0; i < NUM_CHANNELS; i++) {
	for (j = 0; j < NUM_CHANNELS; j++)
	InputGenerate_StrictFP(i, inputs + jNUM_SAMPLES2,
	NUM_SAMPLES);
	for (j = 0; j < NUM_POST_DEC_1; j++)
	BeamFirFilter_StrictFP(&coarse_fir_data[i],
	NUM_SAMPLES, COARSE_DECIMATION_RATIO,
	inputs + j * COARSE_DECIMATION_RATIO*2,
	predec + j * 2);
	for (j = 0; j < NUM_POST_DEC_2; j++)
	BeamFirFilter_StrictFP(&fine_fir_data[i],
	NUM_POST_DEC_1, FINE_DECIMATION_RATIO,
	predec + j * FINE_DECIMATION_RATIO * 2,
	postdec[i] + j * 2);
	}
	/* Assemble beam-forming input: */
	for (i = 0; i < NUM_CHANNELS; i++)
	for (j = 0; j < NUM_POST_DEC_2; j++)
	{
	beam_input[jNUM_CHANNELS2+2i] = postdec[i][2j];
	beam_input[jNUM_CHANNELS2+2i+1] = postdec[i][2j+1];
	}
	for (i = 0; i < NUM_BEAMS; i++)
	{
	/* Have now rearranged NUM_CHANNELS*NUM_POST_DEC_2 items.
	* BeamForm takes NUM_CHANNELS inputs, 2 outputs. */
	for (j = 0; j < NUM_POST_DEC_2; j++)
	BeamForm_StrictFP(i, beam_weights[i],
	beam_input + jNUM_CHANNELS2,
	beam_output + j*2);
	for (j = 0; j < NUM_POST_DEC_2; j++)
	BeamFirFilter_StrictFP(&mf_fir_data[i],
	NUM_POST_DEC_2, 1,
	beam_output+j2, beam_fir_output+j2);
	Magnitude_StrictFP(beam_fir_output, beam_fir_mag, NUM_POST_DEC_2);
	Detector_StrictFP(i, beam_fir_mag, detector_out_StrictFP[i]);
	}
	}

	void InputGenerate_StrictFP(int channel, float *inputs, int n)
	{
	#pragma STDC FP_CONTRACT OFF
	int i;
	for (i = 0; i < n; i++)
	{
	if (channel == TARGET_BEAM && i == TARGET_SAMPLE)
	{
	#ifdef RANDOM_INPUTS
	inputs[2i] = sqrt(ichannel);
	inputs[2i+1] = sqrt(ichannel)+1;
	#else
	inputs[2*i] = sqrt(CFAR_THRESHOLD);
	inputs[2*i+1] = 0;
	#endif
	} else {
	#ifdef RANDOM_INPUTS
	float root = sqrt(i*channel);
	inputs[2*i] = -root;
	inputs[2*i+1] = -(root+1);
	#else
	inputs[2*i] = 0;
	inputs[2*i+1] = 0;
	#endif
	}
	}
	}

	void BeamFirSetup_StrictFP(struct BeamFirData *data, int n)
	{
	#pragma STDC FP_CONTRACT OFF
	int i, j;

	data->len = n;
	data->count = 0;
	data->pos = 0;
	data->weight = malloc(sizeof(float)2n);
	data->buffer = malloc(sizeof(float)2n);

	#ifdef RANDOM_WEIGHTS
	for (j = 0; j < n; j++) {
	int idx = j+1;
	data->weight[j*2] = sin(idx) / ((float)idx);
	data->weight[j*2+1] = cos(idx) / ((float)idx);
	data->buffer[j*2] = 0.0;
	data->buffer[j*2+1] = 0.0;
	}
	#else
	data->weight[0] = 1.0;
	data->weight[1] = 0.0;
	for (i = 1; i < 2*n; i++) {
	data->weight[i] = 0.0;
	data->buffer[i] = 0.0;
	}
	#endif
	}

	void BeamFirFilter_StrictFP(struct BeamFirData *data,
	int input_length, int decimation_ratio,
	float in, float out)
	{
	#pragma STDC FP_CONTRACT OFF
	/* Input must be exactly 2*decimation_ratio long; output must be
	* exactly 2 long. */
	float real_curr = 0;
	float imag_curr = 0;
	int i;
	int modPos;
	int len, mask, mask2;

	len = data->len;
	mask = len - 1;
	mask2 = 2 * len - 1;
	modPos = 2*(len - 1 - data->pos);
	data->buffer[modPos] = in[0];
	data->buffer[modPos+1] = in[1];

	/* Profiling says: this is the single inner loop that matters! */
	for (i = 0; i < 2*len; i+=2) {
	float rd = data->buffer[modPos];
	float id = data->buffer[modPos+1];
	float rw = data->weight[i];
	float iw = data->weight[i+1];
	float rci = rd * rw + id * iw;
	/* sign error? this is consistent with StreamIt --dzm */
	float ici = id * rw + rd * iw;
	real_curr += rci;
	imag_curr += ici;
	modPos = (modPos + 2) & mask2;
	}

	data->pos = (data->pos + 1) & mask;
	out[0] = real_curr;
	out[1] = imag_curr;
	data->count += decimation_ratio;
	if (data->count == input_length)
	{
	data->count = 0;
	data->pos = 0;
	for (i = 0; i < 2*data->len; i++) {
	data->buffer[i] = 0.0;
	}
	}
	}

	void BeamFormWeights_StrictFP(int beam, float *weights)
	{
	#pragma STDC FP_CONTRACT OFF
	int i;
	for (i = 0; i < NUM_CHANNELS; i++)
	{
	#ifdef RANDOM_WEIGHTS
	int idx = i+1;
	weights[2*i] = sin(idx) / (float)(beam+idx);
	weights[2*i+1] = cos(idx) / (float)(beam+idx);
	#else
	if (i == beam) {
	weights[2*i] = 1;
	weights[2*i+1] = 0;
	} else {
	weights[2*i] = 0;
	weights[2*i+1] = 0;
	}
	#endif
	}
	}

	void BeamForm_StrictFP(int beam, const float weights, const float input,
	float *output)
	{
	#pragma STDC FP_CONTRACT OFF
	/* 2NUM_CHANNELS inputs and weights; 2 outputs. /
	float real_curr = 0;
	float imag_curr = 0;
	int i;
	for (i = 0; i < NUM_CHANNELS; i++)
	{
	real_curr += weights[2i] input[2i] - weights[2i+1] * input[2*i+1];
	imag_curr += weights[2i] input[2i+1] + weights[2i+1] * input[2*i];
	}
	output[0] = real_curr;
	output[1] = imag_curr;
	}

	void Magnitude_StrictFP(float in, float out, int n)
	{
	#pragma STDC FP_CONTRACT OFF
	int i;
	/* Should be 2n inputs, n outputs. */
	for (i = 0; i < n; i++)
	out[i] = sqrt(in[2i]in[2i] + in[2i+1]in[2i+1]);
	}

	void Detector_StrictFP(int beam, float data, float output)
	{
	#pragma STDC FP_CONTRACT OFF
	int sample;
	/* Should be exactly NUM_POST_DEC_2 samples. */
	for (sample = 0; sample < NUM_POST_DEC_2; sample++)
	{
	float outputVal;
	if (beam == TARGET_BEAM && sample == TARGET_SAMPLE) {
	if (data[sample] >= 0.1)
	outputVal = beam+1;
	else
	outputVal = 0;
	} else {
	if (data[sample] >= 0.1)
	outputVal = -(beam+1);
	else
	outputVal = 0;
	}
	outputVal = data[sample];
	output[sample]= outputVal;
	}
	}