git.gag.com Git - debian/gnuradio/blob - gr-radar/src/lib/eb-xambi.cc

   1 /* -*- c++ -*- */
   2 /*
   3  * Copyright 2005 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 2, 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 #ifndef _GNU_SOURCE
  24 #define _GNU_SOURCE
  25 #endif
  26
  27 #include <iostream>
  28 #include <string>
  29 #include <fstream>
  30 #include <unistd.h>
  31 #include <stdlib.h>
  32 #include <getopt.h>
  33 #include <boost/scoped_array.hpp>
  34 #include <gr_complex.h>
  35 #include <gr_fxpt_nco.h>
  36 #include "time_series.h"
  37
  38
  39 gr_complex
  40 complex_conj_dotprod(const gr_complex *a, const gr_complex *b, unsigned n)
  41 {
  42   gr_complex acc = 0;
  43   for (unsigned i = 0; i < n; i++)
  44     acc += a[i] * conj(b[i]);
  45
  46   return acc;
  47 }
  48
  49 /*!
  50  * \brief frequency translate src by normalized_freq
  51  *
  52  * \param dst                   destination
  53  * \param src                   source
  54  * \param n                     length of src and dst in samples
  55  * \param normalized_freq       [-1/2, +1/2]
  56  */
  57 void
  58 freq_xlate(gr_complex *dst, const gr_complex *src, unsigned n, float normalized_freq)
  59 {
  60   gr_fxpt_nco   nco;
  61   nco.set_freq(2 * M_PI * normalized_freq);
  62
  63   for (unsigned int i = 0; i < n; i++){
  64     gr_complex  phasor(nco.cos(), nco.sin());
  65     dst[i] = src[i] * phasor;
  66     nco.step();
  67   }
  68 }
  69
  70 inline void
  71 write_float(FILE *output, float x)
  72 {
  73   fwrite(&x, sizeof(x), 1, output);
  74 }
  75
  76
  77 // write 8-float header
  78 static void
  79 write_header (FILE *output, float ndoppler_bins, float min_doppler, float max_doppler)
  80 {
  81   write_float(output, ndoppler_bins);
  82   write_float(output, min_doppler);
  83   write_float(output, max_doppler);
  84   write_float(output, 0);
  85   write_float(output, 0);
  86   write_float(output, 0);
  87   write_float(output, 0);
  88   write_float(output, 0);
  89 }
  90
  91
  92 void
  93 main_loop(FILE *output, time_series &ref_ts, time_series &scat0_ts,
  94           unsigned nranges, unsigned correlation_window_size,
  95           float min_doppler, float max_doppler, int ndoppler_bins)
  96 {
  97   fprintf(stderr, "ndoppler_bins = %10d\n", ndoppler_bins);
  98   fprintf(stderr, "min_doppler   = %10f\n", min_doppler);
  99   fprintf(stderr, "max_doppler   = %10f\n", max_doppler);
 100
 101   // float scale_factor = 1.0/correlation_window_size;  // FIXME, not sure this is right
 102   float scale_factor = 1.0;                             // FIXME, not sure this is right
 103
 104   boost::scoped_array<gr_complex>  shifted(new gr_complex[correlation_window_size]);
 105
 106   // gr_complex shifted[correlation_window_size];               // doppler shifted reference
 107
 108   float doppler_incr = 0;
 109   if (ndoppler_bins == 1){
 110     min_doppler = 0;
 111     max_doppler = 0;
 112   }
 113   else
 114     doppler_incr = (max_doppler - min_doppler) / (ndoppler_bins - 1);
 115
 116   write_header(output, ndoppler_bins, min_doppler, max_doppler);
 117
 118   unsigned long long ro = 0;    // reference offset
 119   unsigned long long so = 0;    // scatter offset
 120
 121   for (unsigned int n = 0; n < nranges; n++){
 122     if (0){
 123       fprintf(stdout, "n =  %6d\n", n);
 124       fprintf(stdout, "ro = %6lld\n", ro);
 125       fprintf(stdout, "so = %6lld\n", so);
 126     }
 127     const gr_complex *ref = (const gr_complex *) ref_ts.seek(ro, correlation_window_size);
 128     const gr_complex *scat0 = (const gr_complex *) scat0_ts.seek(so, correlation_window_size);
 129     if (ref == 0 || scat0 == 0)
 130       return;
 131
 132     for (int nd = 0; nd < ndoppler_bins; nd++){
 133       float fdop = min_doppler + doppler_incr * nd;
 134       //fprintf(stderr, "fdop = %10g\n", fdop);
 135       freq_xlate(&shifted[0], ref, correlation_window_size, fdop);      // generated doppler shifted reference
 136
 137       gr_complex ccor = complex_conj_dotprod(&shifted[0], scat0, correlation_window_size);
 138       float out = norm(ccor) * scale_factor;
 139
 140       // fprintf(output, "%12g\n", out);
 141       write_float(output, out);
 142     }
 143
 144     so += 1;
 145   }
 146 }
 147
 148 static void
 149 usage(const char *argv0)
 150 {
 151   const char *progname;
 152   const char *t = std::strrchr(argv0, '/');
 153   if (t != 0)
 154     progname = t + 1;
 155   else
 156     progname = argv0;
 157
 158   fprintf(stderr, "usage: %s [options] ref_file scatter_file\n", progname);
 159   fprintf(stderr, "    -o OUTPUT_FILENAME [default=eb-xambi.out]\n");
 160   fprintf(stderr, "    -m MIN_RANGE [default=0]\n");
 161   fprintf(stderr, "    -M MAX_RANGE [default=300]\n");
 162   fprintf(stderr, "    -w CORRELATION_WINDOW_SIZE [default=2500]\n");
 163   fprintf(stderr, "    -s NSAMPLES_TO_SKIP [default=0]\n");
 164   fprintf(stderr, "    -d max_doppler (normalized: [0, +1/2)) [default=.0012]\n");
 165   fprintf(stderr, "    -n ndoppler_bins [default=31]\n");
 166 }
 167
 168 int
 169 main(int argc, char **argv)
 170 {
 171   int   ch;
 172   int min_range =    0;
 173   int max_range =  300;
 174   const char *ref_filename = 0;
 175   const char *scatter_filename = 0;
 176   const char *output_filename = "eb-xambi.out";
 177   unsigned int correlation_window_size = 2500;
 178   long long int nsamples_to_skip = 0;
 179   double max_doppler = 0.0012;
 180   int ndoppler_bins = 31;
 181
 182
 183   while ((ch = getopt(argc, argv, "m:M:ho:w:s:d:n:")) != -1){
 184     switch (ch){
 185     case 'm':
 186       min_range = strtol(optarg, 0, 0);
 187       break;
 188
 189     case 'M':
 190       max_range = strtol(optarg, 0, 0);
 191       break;
 192
 193     case 'w':
 194       correlation_window_size = strtol(optarg, 0, 0);
 195       if (correlation_window_size <= 1){
 196         usage(argv[0]);
 197         fprintf(stderr, "    correlation_window_size must be >= 1\n");
 198         exit(1);
 199       }
 200       break;
 201
 202     case 'o':
 203       output_filename = optarg;
 204       break;
 205
 206     case 's':
 207       nsamples_to_skip = (long long) strtod(optarg, 0);
 208       if (nsamples_to_skip < 0){
 209         usage(argv[0]);
 210         fprintf(stderr, "    nsamples_to_skip must be >= 0\n");
 211         exit(1);
 212       }
 213       break;
 214
 215     case 'd':
 216       max_doppler = strtod(optarg, 0);
 217       if (max_doppler < 0 || max_doppler >= 0.5){
 218         usage(argv[0]);
 219         fprintf(stderr, "    max_doppler must be in [0, 0.5)\n");
 220         exit(1);
 221       }
 222       break;
 223
 224     case 'n':
 225       ndoppler_bins = strtol(optarg, 0, 0);
 226       if (ndoppler_bins < 1){
 227         usage(argv[0]);
 228         fprintf(stderr, "    ndoppler_bins must >= 1\n");
 229         exit(1);
 230       }
 231       break;
 232
 233     case '?':
 234     case 'h':
 235     default:
 236       usage(argv[0]);
 237       exit(1);
 238     }
 239   } // while getopt
 240
 241   if (argc - optind != 2){
 242     usage(argv[0]);
 243     exit(1);
 244   }
 245
 246   if (max_range < min_range){
 247     usage(argv[0]);
 248     fprintf(stderr, "    max_range must be >= min_range\n");
 249     exit(1);
 250   }
 251   unsigned int nranges = max_range - min_range + 1;
 252
 253   ref_filename = argv[optind++];
 254   scatter_filename = argv[optind++];
 255
 256   FILE *output = fopen(output_filename, "wb");
 257   if (output == 0){
 258     perror(output_filename);
 259     exit(1);
 260   }
 261
 262   unsigned long long ref_starting_offset = 0;
 263   unsigned long long scatter_starting_offset = 0;
 264
 265   if (min_range < 0){
 266     ref_starting_offset = -min_range;
 267     scatter_starting_offset = 0;
 268   }
 269   else {
 270     ref_starting_offset = 0;
 271     scatter_starting_offset = min_range;
 272   }
 273
 274   ref_starting_offset += nsamples_to_skip;
 275   scatter_starting_offset += nsamples_to_skip;
 276
 277   try {
 278     time_series ref(sizeof(gr_complex), ref_filename, ref_starting_offset, 0);
 279     time_series scat0(sizeof(gr_complex), scatter_filename, scatter_starting_offset, 0);
 280
 281     main_loop(output, ref, scat0, nranges, correlation_window_size,
 282               -max_doppler, max_doppler, ndoppler_bins);
 283   }
 284   catch (std::string &s){
 285     std::cerr << s << std::endl;
 286     exit(1);
 287   }
 288
 289   return 0;
 290 }
 291