3 * Copyright 2006 Free Software Foundation, Inc.
5 * This file is part of GNU Radio
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)
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.
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.
28 #include <gr_complex.h>
30 #include "time_series.h"
36 * \file xambi.cc driver for computation of cross ambiguity
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
43 * Melissa Meyer mgmeyer@ee.washington.edu
45 * Extensively revised since then.
48 static int default_decimation = 40;
49 static int default_fftsize = 256;
50 static int default_naverages = 1000000; // infinite
52 /// usage for the xambi program
54 static void usage(char *argv0)
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;
70 cerr << "The reffile and scatterfile are native-endian binary complex<float>" << endl;
71 cerr << "and must be sampled at the same rate." << endl;
77 complex_dot_product(const gr_complex *xx, const gr_complex *yy, int nterms)
81 for (int i = 0; i < nterms; i++)
82 sum += xx[i] * yy[i]; // FIXME? conj(yy[i])
91 main(int argc, char *argv[])
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;
103 unsigned long long starting_file_offset = 0;
104 unsigned long long nsamples_to_process = (unsigned long long) -1;
108 const gr_complex *x, *y;
109 const gr_complex *xx, *yy;
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
116 while((c = getopt(argc,argv,"a:o:x:y:r:d:f:n:hvS:N:")) != EOF) {
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;
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();
139 cerr << "you specified -r " << range << "; must be non-negative (exit)" << endl;
144 cerr << "you specified -n " << nranges << "; must be positive (exit)" << endl;
149 cerr << "you specified -d " << decimation << "; must be positive (exit)" << endl;
154 cerr << "you specified -a " << naverages << "; must be positive (exit)" << endl;
159 cerr << "you must specify a reference signal with the -x option" << endl;
163 if (optind >= argc) {
164 cerr << "you must specify a scattering signal after all other options" << endl;
168 time_series X(sizeof(gr_complex), ref_fname,
169 starting_file_offset, nsamples_to_process);
171 time_series Y(sizeof(gr_complex), argv[optind],
172 starting_file_offset, nsamples_to_process); // add more for interferometry ...
174 float psd[fftsize*nranges];
178 snprintf(fname, sizeof(fname), "%s.out", "xambi");
182 ofstream Z(out_fname);
184 blocksize = fftsize*decimation + nranges;
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);
192 for(i = 0; i < nranges*fftsize; i++)
196 while(1){ // loop over data until exhausted.
198 cerr << " " << a; // write out the number of completed averages
202 x = (const gr_complex *) X.seek(offset, blocksize);
203 y = (const gr_complex *) Y.seek(offset + range, blocksize);
205 if(!x || !y) // ran out of data; stop integrating
208 for(r = 0; r < nranges; r++) { // For Each Range ...
212 for(f = 0; f < fftsize; f++) { // and for each Doppler bin ...
214 // cross correlate and do a boxcar decimation
216 fft_input[f] = complex_dot_product(xx, yy, decimation);
221 fft.execute(); // input: fft_input; output: fft_output
223 for(f = 0; f < fftsize; f++) {
224 psd[r*fftsize + f] += norm(fft_output[f]);
229 offset += fftsize * decimation;
233 cerr << " dumping " << endl;
235 for(i = 0; i < nranges*fftsize; i++) // normalize
236 psd[i] *= scale_factor;
238 Z.write((const char *) psd, nranges*fftsize*sizeof(float));
240 for(i = 0; i < nranges*fftsize; i++)
251 for(i = 0; i < nranges*fftsize; i++) // normalize
252 psd[i] *= scale_factor;
254 Z.write((const char *) psd, nranges*fftsize*sizeof(float));