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