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