git.gag.com Git - debian/gnuradio/blob - gr-radar/src/lib/sim-airplane2.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., 59 Temple Place - Suite 330,
  20  * Boston, MA 02111-1307, 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 <gr_complex.h>
  33 #include <getopt.h>
  34 #include <gr_misc.h>
  35 #include <limits>
  36 #include <gr_fxpt_nco.h>
  37 #include "time_series.h"
  38 #include "simulation.h"
  39
  40 static const double C = 3e8;    // sped of light, m/s
  41
  42
  43 // ------------------------------------------------------------------------
  44
  45 class delay_line {
  46   std::vector<gr_complex>       d_z;
  47   const int                     d_mask;
  48   int                           d_newest;
  49 public:
  50   delay_line(unsigned int max_delay)
  51     : d_z(gr_rounduppow2(max_delay)), d_mask(d_z.size()-1), d_newest(0)
  52   {
  53   }
  54
  55   void
  56   push_item(gr_complex x)
  57   {
  58     d_newest = (d_newest - 1) & d_mask;
  59     d_z[d_newest] = x;
  60   }
  61
  62   gr_complex
  63   ref_item(int delay) const
  64   {
  65     return d_z[(d_newest + delay) & d_mask];
  66   }
  67 };
  68
  69 // ------------------------------------------------------------------------
  70
  71 class aux_state {
  72 public:
  73   dyn_object    *d_obj;
  74   double         d_last_slant_range;
  75   gr_fxpt_nco    d_nco;
  76
  77   aux_state(dyn_object *obj) : d_obj(obj) {}
  78 };
  79
  80 // ------------------------------------------------------------------------
  81
  82 class my_sim : public simulation
  83 {
  84   FILE                   *d_output;
  85   time_series            &d_ref;
  86   unsigned long long      d_pos;                // position in time series
  87   delay_line              d_z;
  88   dyn_object             *d_tx;         // transmitter (not moving)
  89   dyn_object             *d_rx0;                // receiver (not moving)
  90   std::vector<aux_state*> d_target;
  91
  92   double                  d_baseline;           // length of baseline in meters
  93   double                  d_range_bin;          // meters/range_bin
  94   float                   d_tx_lambda;          // wavelength of tx signals in meters
  95   float                   d_sample_rate;
  96   float                   d_gain;               // linear scale factor
  97
  98   void adjust_for_start_time(double start_time);
  99
 100   bool write_output(gr_complex x)
 101   {
 102     return fwrite(&x, sizeof(x), 1, d_output) == 1;
 103   }
 104
 105 public:
 106   my_sim(FILE *output, time_series &ref, double timestep, float sample_rate,
 107          double start_time, double tx_freq, double gain_db);
 108   ~my_sim();
 109
 110   bool update();
 111   bool run(long long nsteps);
 112 };
 113
 114
 115 my_sim::my_sim(FILE *output, time_series &ref, double timestep,
 116                float sample_rate, double start_time,
 117                double tx_freq, double gain_db)
 118   : simulation(timestep),
 119     d_output(output), d_ref(ref), d_pos(0), d_z(1024),
 120     d_range_bin(C * timestep), d_tx_lambda(C/tx_freq),
 121     d_sample_rate(sample_rate), d_gain(exp10(gain_db/10))
 122 {
 123   d_tx = new dyn_object(point(0,0), point(0,0), "Tx");
 124   d_rx0 = new dyn_object(point(45e3,0), point(0,0), "Rx0");
 125
 126   add_object(d_tx);
 127   add_object(d_rx0);
 128   d_baseline = dyn_object::distance(*d_tx, *d_rx0);
 129
 130   {
 131     // add targets
 132     float aircraft_speed;
 133     float aircraft_angle;
 134     point aircraft_pos;
 135     dyn_object  *ac;
 136
 137     // target 1
 138     aircraft_speed = 135;                       // m/s
 139     aircraft_angle = 240 * M_PI/180;
 140     aircraft_pos = point(55e3, 20e3);
 141
 142     ac = new dyn_object(aircraft_pos,
 143                         point(aircraft_speed * cos(aircraft_angle),
 144                               aircraft_speed * sin(aircraft_angle)),
 145                         "Ac0");
 146     add_object(ac);
 147     d_target.push_back(new aux_state(ac));
 148
 149     // target 2
 150     aircraft_speed = 350;                       // m/s
 151     aircraft_angle = 0 * M_PI/180;
 152     aircraft_pos = point(-20e3, 60e3);
 153
 154     ac = new dyn_object(aircraft_pos,
 155                         point(aircraft_speed * cos(aircraft_angle),
 156                               aircraft_speed * sin(aircraft_angle)),
 157                         "Ac1");
 158     add_object(ac);
 159     d_target.push_back(new aux_state(ac));
 160   }
 161
 162   adjust_for_start_time(start_time);
 163
 164   for (unsigned i = 0; i < d_target.size(); i++)
 165     d_target[i]->d_last_slant_range =
 166       (dyn_object::distance(*d_tx, *d_target[i]->d_obj)
 167        + dyn_object::distance(*d_target[i]->d_obj, *d_rx0));
 168
 169 }
 170
 171 my_sim::~my_sim()
 172 {
 173 }
 174
 175 void
 176 my_sim::adjust_for_start_time(double start_time)
 177 {
 178   for (unsigned i = 0; i < d_obj.size(); i++){
 179     // Adjust initial starting positions depending on simulation
 180     // start time.  FIXME Assumes velocity is constant
 181     point p = d_obj[i]->pos();
 182     point v = d_obj[i]->vel();
 183     p.set_x(p.x() + v.x() * start_time);
 184     p.set_y(p.y() + v.y() * start_time);
 185     d_obj[i]->set_pos(p);
 186   }
 187 }
 188
 189 bool
 190 my_sim::update()
 191 {
 192   // std::cout << *d_ac0 << std::endl;
 193
 194   // grab new item from input and insert it into delay line
 195   const gr_complex *in = (const gr_complex *) d_ref.seek(d_pos++, 1);
 196   if (in == 0)
 197     return false;
 198   d_z.push_item(*in);
 199
 200   gr_complex s = 0;     // output sample
 201   // FIXME ought to add in attenuated direct path input
 202
 203
 204   // for each target, compute slant_range and slant_range'
 205
 206   for (unsigned i = 0; i < d_target.size(); i++){
 207     aux_state *t = d_target[i];
 208
 209     double slant_range =
 210       (dyn_object::distance(*d_tx, *t->d_obj)
 211        + dyn_object::distance(*t->d_obj, *d_rx0));                        // meters
 212
 213     double delta_slant_range = slant_range - t->d_last_slant_range;
 214     t->d_last_slant_range = slant_range;
 215     double deriv_slant_range_wrt_time = delta_slant_range / timestep();   // m/sec
 216
 217     //fprintf(stdout, "%10.3f\t%10.3f\n", slant_range, deriv_slant_range_wrt_time);
 218
 219     // FIXME, may want to interpolate between two bins.
 220     int int_delay = lrint((slant_range - d_baseline) / d_range_bin);
 221
 222     gr_complex x = d_z.ref_item(int_delay);
 223
 224     // scale amplitude (this includes everything: RCS, antenna gain, losses, etc...)
 225     x = x * d_gain;
 226
 227     if (1){
 228       // compute doppler and apply it
 229       float f_doppler = -deriv_slant_range_wrt_time / d_tx_lambda;
 230
 231       t->d_nco.set_freq(f_doppler / d_sample_rate);
 232       gr_complex phasor(t->d_nco.cos(), t->d_nco.sin());
 233       x = x * phasor;
 234       t->d_nco.step();
 235     }
 236
 237     s += x;             // add in this target's contribution
 238   }
 239
 240   write_output(s);
 241
 242   return simulation::update();          // run generic update
 243 }
 244
 245 bool
 246 my_sim::run(long long nsteps)
 247 {
 248   //fprintf(stdout, "<%12.2f, %12.2f>\n", d_ac0->pos().x(), d_ac0->pos().y());
 249   //std::cout << *d_ac0 << std::endl;
 250   bool ok = simulation::run(nsteps);
 251   //std::cout << *d_ac0 << std::endl;
 252   //fprintf(stdout, "<%12.2f, %12.2f>\n", d_ac0->pos().x(), d_ac0->pos().y());
 253   return ok;
 254 }
 255
 256 // ------------------------------------------------------------------------
 257
 258 static void
 259 usage(const char *argv0)
 260 {
 261   const char *progname;
 262   const char *t = std::strrchr(argv0, '/');
 263   if (t != 0)
 264     progname = t + 1;
 265   else
 266     progname = argv0;
 267
 268   fprintf(stderr, "usage: %s [options] ref_file\n", progname);
 269   fprintf(stderr, "    -o OUTPUT_FILENAME [default=sim.dat]\n");
 270   fprintf(stderr, "    -n NSAMPLES_TO_PRODUCE [default=+inf]\n");
 271   fprintf(stderr, "    -s NSAMPLES_TO_SKIP [default=0]\n");
 272   fprintf(stderr, "    -g reflection gain in dB (should be <= 0) [default=0]\n");
 273   fprintf(stderr, "    -f transmitter freq in Hz [default=100MHz]\n");
 274   fprintf(stderr, "    -r sample rate in Hz [default=250kHz]\n");
 275   fprintf(stderr, "    -S simulation start time in seconds [default=0]\n");
 276 }
 277
 278 int
 279 main(int argc, char **argv)
 280 {
 281   int   ch;
 282   const char *output_filename = "sim.dat";
 283   const char *ref_filename = 0;
 284   long long int nsamples_to_skip = 0;
 285   long long int nsamples_to_produce = std::numeric_limits<long long int>::max();
 286   double sample_rate = 250e3;
 287   double gain_db = 0;
 288   double tx_freq = 100e6;
 289   double start_time = 0;
 290
 291   while ((ch = getopt(argc, argv, "o:s:n:g:f:S:")) != -1){
 292     switch (ch){
 293     case 'o':
 294       output_filename = optarg;
 295       break;
 296
 297     case 's':
 298       nsamples_to_skip = (long long) strtof(optarg, 0);
 299       if (nsamples_to_skip < 0){
 300         usage(argv[0]);
 301         fprintf(stderr, "    nsamples_to_skip must be >= 0\n");
 302         exit(1);
 303       }
 304       break;
 305
 306     case 'n':
 307       nsamples_to_produce = (long long) strtof(optarg, 0);
 308       if (nsamples_to_produce < 0){
 309         usage(argv[0]);
 310         fprintf(stderr, "    nsamples_to_produce must be >= 0\n");
 311         exit(1);
 312       }
 313       break;
 314
 315     case 'g':
 316       gain_db = strtof(optarg, 0);
 317       break;
 318
 319     case 'f':
 320       tx_freq = strtof(optarg, 0);
 321       break;
 322
 323     case 'r':
 324       sample_rate = strtof(optarg, 0);
 325       break;
 326
 327     case 'S':
 328       start_time = strtof(optarg, 0);
 329       break;
 330
 331     case '?':
 332     case 'h':
 333     default:
 334       usage(argv[0]);
 335       exit(1);
 336     }
 337   } // while getopt
 338
 339   if (argc - optind != 1){
 340     usage(argv[0]);
 341     exit(1);
 342   }
 343
 344   ref_filename = argv[optind++];
 345
 346   double timestep = 1.0/sample_rate;
 347
 348
 349   FILE *output = fopen(output_filename, "wb");
 350   if (output == 0){
 351     perror(output_filename);
 352     exit(1);
 353   }
 354
 355   unsigned long long ref_starting_offset = 0;
 356   ref_starting_offset += nsamples_to_skip;
 357
 358   try {
 359     time_series ref(sizeof(gr_complex), ref_filename, ref_starting_offset, 0);
 360
 361     my_sim      simulator(output, ref, timestep, sample_rate, start_time,
 362                           tx_freq, gain_db);
 363     simulator.run(nsamples_to_produce);
 364   }
 365   catch (std::string &s){
 366     std::cerr << s << std::endl;
 367     exit(1);
 368   }
 369
 370   return 0;
 371 }
 372