git.gag.com Git - debian/gnuradio/blob - gr-radar/src/lib/sim-airplane.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 my_sim : public simulation
  72 {
  73   FILE                  *d_output;
  74   time_series           &d_ref;
  75   unsigned long long     d_pos;         // position in time series
  76   delay_line             d_z;
  77   dyn_object            *d_tx;          // transmitter (not moving)
  78   dyn_object            *d_rx0;         // receiver (not moving)
  79   dyn_object            *d_ac0;         // aircraft (linear motion)
  80   gr_fxpt_nco            d_nco0;
  81
  82   double                d_baseline;             // length of baseline in meters
  83   double                d_last_slant_range;
  84   double                d_range_bin;            // meters/range_bin
  85   float                 d_tx_lambda;            // wavelength of tx signals in meters
  86   float                 d_sample_rate;
  87   float                 d_gain;                 // linear scale factor
  88
  89 public:
  90   my_sim(FILE *output, time_series &ref, double timestep, float sample_rate,
  91          double tx_freq, double gain_db);
  92   ~my_sim();
  93
  94   bool update();
  95   bool run(long long nsteps);
  96
  97   bool write_output(gr_complex x)
  98   {
  99     return fwrite(&x, sizeof(x), 1, d_output) == 1;
 100   }
 101 };
 102
 103 my_sim::my_sim(FILE *output, time_series &ref, double timestep,
 104                float sample_rate, double tx_freq, double gain_db)
 105   : simulation(timestep),
 106     d_output(output), d_ref(ref), d_pos(0), d_z(1024),
 107     d_range_bin(C * timestep), d_tx_lambda(C/tx_freq),
 108     d_sample_rate(sample_rate), d_gain(exp10(gain_db/10))
 109 {
 110   d_tx = new dyn_object(point(0,0), point(0,0), "Tx");
 111   d_rx0 = new dyn_object(point(45e3,0), point(0,0), "Rx0");
 112
 113   //float aircraft_speed =  135;         // meters/sec (~ 300 miles/hr)
 114   float aircraft_speed =  350;         // meters/sec (~ 750 miles/hr)
 115   float aircraft_angle =  250 * M_PI/180;
 116   //point aircraft_pos = point(55e3, 20e3);
 117   point aircraft_pos = point(55e3-5.54e3, 20e3-15.23e3);
 118   d_ac0 = new dyn_object(aircraft_pos,
 119                          point(aircraft_speed * cos(aircraft_angle),
 120                                aircraft_speed * sin(aircraft_angle)),
 121                          "Ac0");
 122   add_object(d_tx);
 123   add_object(d_rx0);
 124   add_object(d_ac0);
 125
 126   d_baseline = dyn_object::distance(*d_tx, *d_rx0);
 127   d_last_slant_range =
 128     dyn_object::distance(*d_tx, *d_ac0) + dyn_object::distance(*d_ac0, *d_rx0);
 129 }
 130
 131 my_sim::~my_sim()
 132 {
 133 }
 134
 135 bool
 136 my_sim::update()
 137 {
 138   // std::cout << *d_ac0 << std::endl;
 139
 140   // compute slant_range and slant_range'
 141   double slant_range =
 142     dyn_object::distance(*d_tx, *d_ac0) + dyn_object::distance(*d_ac0, *d_rx0); // meters
 143   double delta_slant_range = slant_range - d_last_slant_range;
 144   d_last_slant_range = slant_range;
 145   double deriv_slant_range_wrt_time = delta_slant_range / timestep();           // m/sec
 146
 147   // fprintf(stdout, "%10.3f\t%10.3f\n", slant_range, deriv_slant_range_wrt_time);
 148
 149   // grab new item from input and insert it into delay line
 150   const gr_complex *in = (const gr_complex *) d_ref.seek(d_pos++, 1);
 151   if (in == 0)
 152     return false;
 153   d_z.push_item(*in);
 154
 155   // FIXME, may want to interpolate between two bins.
 156   int int_delay = lrint((slant_range - d_baseline) / d_range_bin);
 157
 158   gr_complex x = d_z.ref_item(int_delay);
 159
 160   x = x * d_gain;               // scale amplitude (this includes everything: RCS, antenna gain, losses, etc...)
 161
 162   // compute doppler and apply it
 163   float f_doppler = -deriv_slant_range_wrt_time / d_tx_lambda;
 164   fprintf(stdout, "f_dop: %10.3f\n", f_doppler);
 165
 166   d_nco0.set_freq(f_doppler / d_sample_rate);
 167   gr_complex phasor(d_nco0.cos(), d_nco0.sin());
 168   // x = x * phasor;
 169   d_nco0.step();
 170
 171   write_output(x);
 172
 173   return simulation::update();          // run generic update
 174 }
 175
 176 bool
 177 my_sim::run(long long nsteps)
 178 {
 179   //fprintf(stdout, "<%12.2f, %12.2f>\n", d_ac0->pos().x(), d_ac0->pos().y());
 180   //std::cout << *d_ac0 << std::endl;
 181   bool ok = simulation::run(nsteps);
 182   //std::cout << *d_ac0 << std::endl;
 183   //fprintf(stdout, "<%12.2f, %12.2f>\n", d_ac0->pos().x(), d_ac0->pos().y());
 184   return ok;
 185 }
 186
 187 // ------------------------------------------------------------------------
 188
 189 static void
 190 usage(const char *argv0)
 191 {
 192   const char *progname;
 193   const char *t = std::strrchr(argv0, '/');
 194   if (t != 0)
 195     progname = t + 1;
 196   else
 197     progname = argv0;
 198
 199   fprintf(stderr, "usage: %s [options] ref_file\n", progname);
 200   fprintf(stderr, "    -o OUTPUT_FILENAME [default=sim.dat]\n");
 201   fprintf(stderr, "    -n NSAMPLES_TO_PRODUCE [default=+inf]\n");
 202   fprintf(stderr, "    -s NSAMPLES_TO_SKIP [default=0]\n");
 203   fprintf(stderr, "    -g reflection gain in dB (should be <= 0) [default=0]\n");
 204   fprintf(stderr, "    -f transmitter freq in Hz [default=100MHz]\n");
 205   fprintf(stderr, "    -r sample rate in Hz [default=250kHz]\n");
 206 }
 207
 208 int
 209 main(int argc, char **argv)
 210 {
 211   int   ch;
 212   const char *output_filename = "sim.dat";
 213   const char *ref_filename = 0;
 214   long long int nsamples_to_skip = 0;
 215   long long int nsamples_to_produce = std::numeric_limits<long long int>::max();
 216   double sample_rate = 250e3;
 217   double gain_db = 0;
 218   double tx_freq = 100e6;
 219
 220   while ((ch = getopt(argc, argv, "o:s:n:g:f:")) != -1){
 221     switch (ch){
 222     case 'o':
 223       output_filename = optarg;
 224       break;
 225
 226     case 's':
 227       nsamples_to_skip = (long long) strtof(optarg, 0);
 228       if (nsamples_to_skip < 0){
 229         usage(argv[0]);
 230         fprintf(stderr, "    nsamples_to_skip must be >= 0\n");
 231         exit(1);
 232       }
 233       break;
 234
 235     case 'n':
 236       nsamples_to_produce = (long long) strtof(optarg, 0);
 237       if (nsamples_to_produce < 0){
 238         usage(argv[0]);
 239         fprintf(stderr, "    nsamples_to_produce must be >= 0\n");
 240         exit(1);
 241       }
 242       break;
 243
 244     case 'g':
 245       gain_db = strtof(optarg, 0);
 246       break;
 247
 248     case 'f':
 249       tx_freq = strtof(optarg, 0);
 250       break;
 251
 252     case 'r':
 253       sample_rate = strtof(optarg, 0);
 254       break;
 255
 256     case '?':
 257     case 'h':
 258     default:
 259       usage(argv[0]);
 260       exit(1);
 261     }
 262   } // while getopt
 263
 264   if (argc - optind != 1){
 265     usage(argv[0]);
 266     exit(1);
 267   }
 268
 269   ref_filename = argv[optind++];
 270
 271   double timestep = 1.0/sample_rate;
 272
 273
 274   FILE *output = fopen(output_filename, "wb");
 275   if (output == 0){
 276     perror(output_filename);
 277     exit(1);
 278   }
 279
 280   unsigned long long ref_starting_offset = 0;
 281   ref_starting_offset += nsamples_to_skip;
 282
 283   try {
 284     time_series ref(sizeof(gr_complex), ref_filename, ref_starting_offset, 0);
 285
 286     my_sim      simulator(output, ref, timestep, sample_rate, tx_freq, gain_db);
 287     simulator.run(nsamples_to_produce);
 288   }
 289   catch (std::string &s){
 290     std::cerr << s << std::endl;
 291     exit(1);
 292   }
 293
 294   return 0;
 295 }
 296