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

   1 /* -*- c++ -*- */
   2 /*
   3  * Copyright 2006 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 #include <iostream>
  24 #include <string>
  25 #include <fstream>
  26 #include <unistd.h>
  27
  28 #include <gr_complex.h>
  29 #include <gri_fft.h>
  30 #include "time_series.h"
  31
  32
  33 using namespace std;
  34
  35 /*!
  36  * \file xambi.cc driver for computation of cross ambiguity
  37  *
  38  * Based on ideas liberally lifted from a version of xambi.cc
  39  * obtained from John Sahr, that identified these people as authors:
  40  *   John Sahr            jdsahr@u.washington.edu
  41  *   Frank Lind           flind@haystack.mit.edu
  42  *   Chucai "Cliff" Zhou
  43  *   Melissa Meyer        mgmeyer@ee.washington.edu
  44  *
  45  * Extensively revised since then.
  46  */
  47
  48 static int default_decimation = 40;
  49 static int default_fftsize    = 256;
  50 static int default_naverages   = 1000000; // infinite
  51
  52 /// usage for the xambi program
  53
  54 static void usage(char *argv0)
  55 {
  56   cerr << "usage:  xambi [opts] scatterfile"                << endl;
  57   cerr << "  where [opts] are"                              << endl;
  58   cerr << "  -x reffilename                  : Required"    << endl;
  59   cerr << "  -o outputfilename[=xambi.out]   : create or overwrite" << endl;
  60   cerr << "  -a averages[=100000]            : restart accumulation; 100000 = infinite" << endl;
  61   cerr << "  -d decimationfactor[=40]"                      << endl;
  62   cerr << "  -f fftsize[=256]"                              << endl;
  63   cerr << "  -r range[=0]                    : first range" << endl;
  64   cerr << "  -n nranges[=1]                  : range count" << endl;
  65   cerr << "  -S start[=0]                    : starting offset in file" << endl;
  66   cerr << "  -N nsamples[=inf]               : # of samples to process" << endl;
  67   cerr << "  -v                              : increment verbosity" << endl;
  68   cerr << "  -h                              : be helpful" << endl;
  69   cerr << endl;
  70   cerr << "The reffile and scatterfile are native-endian binary complex<float>" << endl;
  71   cerr << "and must be sampled at the same rate." << endl;
  72
  73   exit(0);
  74 }
  75
  76 gr_complex
  77 complex_dot_product(const gr_complex *xx, const gr_complex *yy, int nterms)
  78 {
  79   gr_complex sum(0);
  80
  81   for (int i = 0; i < nterms; i++)
  82     sum += xx[i] * yy[i];               // FIXME?   conj(yy[i])
  83
  84   return sum;
  85 }
  86
  87
  88 /// the main driver
  89
  90 int
  91 main(int argc, char *argv[])
  92 {
  93   char *ref_fname = 0;    //< holds name of reference signal source
  94   char *out_fname = 0;    //< holds name of processed output
  95   int decimation  = default_decimation;
  96   int range       = 0;    //< first range of range block
  97   int nranges     = 1;    //< number of ranges of range block
  98   int fftsize     = default_fftsize;
  99   int naverages   = default_naverages;
 100   int verbosity   = 0;
 101   int blocksize   = 0;
 102   int offset      = 0;
 103   unsigned long long starting_file_offset = 0;
 104   unsigned long long nsamples_to_process = (unsigned long long) -1;
 105
 106   int f;
 107
 108   const gr_complex     *x,  *y;
 109   const gr_complex     *xx, *yy;
 110
 111   int c;         // used to process the command line
 112   int r;         // an index to count over ranges
 113   int i;         // an index to count through the time series
 114   int a;
 115
 116   while((c = getopt(argc,argv,"a:o:x:y:r:d:f:n:hvS:N:")) != EOF) {
 117     switch(c) {
 118     case 'x': ref_fname  = optarg;            break;
 119     case 'o': out_fname  = optarg;            break;
 120     case 'a': naverages   = atoi(optarg);     break;
 121     case 'd': decimation = atoi(optarg);      break;
 122     case 'r': range      = atoi(optarg);      break;
 123     case 'n': nranges    = atoi(optarg);      break;
 124     case 'f': fftsize    = atoi(optarg);      break;
 125     case 'S': starting_file_offset = strtoll(optarg, 0, 0);  break;
 126     case 'N': nsamples_to_process = strtoll(optarg, 0, 0);   break;
 127     case 'v': verbosity++;                    break;
 128     case 'h': usage(argv[0]);                 break;
 129     default:  usage(argv[0]);                 break;
 130     }
 131   }
 132
 133   // Wrapper for FFTW 1d forward FFT.  N.B. output is not scaled by 1/fftsize
 134   gri_fft_complex       fft(fftsize, true);
 135   gr_complex            *fft_input = fft.get_inbuf();
 136   gr_complex            *fft_output = fft.get_outbuf();
 137
 138   if(range < 0) {
 139     cerr << "you specified -r " << range << "; must be non-negative (exit)" << endl;
 140     exit(1);
 141   }
 142
 143   if(nranges < 1) {
 144     cerr << "you specified -n " << nranges << "; must be positive (exit)" << endl;
 145     exit(1);
 146   }
 147
 148   if(decimation < 1) {
 149     cerr << "you specified -d " << decimation << "; must be positive (exit)" << endl;
 150     exit(1);
 151   }
 152
 153   if(naverages < 1) {
 154     cerr << "you specified -a " << naverages << "; must be positive (exit)" << endl;
 155     exit(1);
 156   }
 157
 158   if(ref_fname == 0) {
 159     cerr << "you must specify a reference signal with the -x option" << endl;
 160     usage(argv[0]);
 161   }
 162
 163   if (optind >= argc) {
 164     cerr << "you must specify a scattering signal after all other options" << endl;
 165     usage(argv[0]);
 166   }
 167
 168   time_series   X(sizeof(gr_complex), ref_fname,
 169                   starting_file_offset, nsamples_to_process);
 170
 171   time_series   Y(sizeof(gr_complex), argv[optind],
 172                   starting_file_offset, nsamples_to_process); // add more for interferometry ...
 173
 174   float  psd[fftsize*nranges];
 175
 176   if(out_fname == 0) {
 177     char fname[200];
 178     snprintf(fname, sizeof(fname), "%s.out", "xambi");
 179     out_fname = fname;
 180   }
 181
 182   ofstream Z(out_fname);
 183
 184   blocksize = fftsize*decimation + nranges;
 185   offset    = 0;
 186   a         = 0;
 187
 188   // the fftsize is squared because we're using norm, not abs,
 189   // when computing the psd
 190   float scale_factor = 1.0 / (fftsize * fftsize);
 191
 192   for(i = 0; i < nranges*fftsize; i++)
 193     psd[i] = 0.0;
 194
 195
 196   while(1){              // loop over data until exhausted.
 197     if(verbosity > 1) {
 198       cerr << " " << a;  // write out the number of completed averages
 199       cerr.flush();
 200     }
 201
 202     x = (const gr_complex *) X.seek(offset,         blocksize);
 203     y = (const gr_complex *) Y.seek(offset + range, blocksize);
 204
 205     if(!x || !y)          // ran out of data; stop integrating
 206       break;
 207
 208     for(r = 0; r < nranges; r++) {  // For Each Range ...
 209       xx = x;
 210       yy = y + r;
 211
 212       for(f = 0; f < fftsize; f++) {  // and for each Doppler bin ...
 213
 214         // cross correlate and do a boxcar decimation
 215
 216         fft_input[f] = complex_dot_product(xx, yy, decimation);
 217         xx += decimation;
 218         yy += decimation;
 219       }
 220
 221       fft.execute();    // input: fft_input; output: fft_output
 222
 223       for(f = 0; f < fftsize; f++) {
 224         psd[r*fftsize + f] += norm(fft_output[f]);
 225       }
 226     }   // end range
 227
 228     a++;
 229     offset += fftsize * decimation;
 230
 231     if(a >= naverages) {
 232       if(verbosity > 0)
 233         cerr << " dumping " << endl;
 234
 235       for(i = 0; i < nranges*fftsize; i++)      // normalize
 236         psd[i] *= scale_factor;
 237
 238       Z.write((const char *) psd, nranges*fftsize*sizeof(float));
 239
 240       for(i = 0; i < nranges*fftsize; i++)
 241         psd[i] = 0.0;
 242
 243       a = 0;
 244     }
 245   }
 246
 247   if(verbosity > 1)
 248     printf("\n");
 249
 250   if(a > 0) {
 251     for(i = 0; i < nranges*fftsize; i++)        // normalize
 252       psd[i] *= scale_factor;
 253
 254     Z.write((const char *) psd, nranges*fftsize*sizeof(float));
 255   }
 256
 257   return 0;
 258 }
 259