From 5d69a524f81f234b3fbc41d49ba18d6f6886baba Mon Sep 17 00:00:00 2001 From: jcorgan Date: Thu, 3 Aug 2006 04:51:51 +0000 Subject: Houston, we have a trunk. git-svn-id: http://gnuradio.org/svn/gnuradio/trunk@3122 221aa14e-8319-0410-a670-987f0aec2ac5 --- gr-radar/src/lib/Makefile.am | 118 +++++++++++ gr-radar/src/lib/eb-xambi.cc | 291 ++++++++++++++++++++++++++ gr-radar/src/lib/gen_run | 26 +++ gr-radar/src/lib/plot_rd.m | 38 ++++ gr-radar/src/lib/plot_xambi.m | 15 ++ gr-radar/src/lib/plot_xambi_to_file.m | 17 ++ gr-radar/src/lib/plot_xambi_tool.m | 8 + gr-radar/src/lib/sim-airplane.cc | 296 +++++++++++++++++++++++++++ gr-radar/src/lib/sim-airplane2.cc | 372 ++++++++++++++++++++++++++++++++++ gr-radar/src/lib/simulation.cc | 89 ++++++++ gr-radar/src/lib/simulation.h | 116 +++++++++++ gr-radar/src/lib/time_series.cc | 96 +++++++++ gr-radar/src/lib/time_series.h | 73 +++++++ gr-radar/src/lib/xambi.cc | 259 +++++++++++++++++++++++ 14 files changed, 1814 insertions(+) create mode 100644 gr-radar/src/lib/Makefile.am create mode 100644 gr-radar/src/lib/eb-xambi.cc create mode 100755 gr-radar/src/lib/gen_run create mode 100644 gr-radar/src/lib/plot_rd.m create mode 100755 gr-radar/src/lib/plot_xambi.m create mode 100755 gr-radar/src/lib/plot_xambi_to_file.m create mode 100755 gr-radar/src/lib/plot_xambi_tool.m create mode 100644 gr-radar/src/lib/sim-airplane.cc create mode 100644 gr-radar/src/lib/sim-airplane2.cc create mode 100644 gr-radar/src/lib/simulation.cc create mode 100644 gr-radar/src/lib/simulation.h create mode 100644 gr-radar/src/lib/time_series.cc create mode 100644 gr-radar/src/lib/time_series.h create mode 100644 gr-radar/src/lib/xambi.cc (limited to 'gr-radar/src/lib') diff --git a/gr-radar/src/lib/Makefile.am b/gr-radar/src/lib/Makefile.am new file mode 100644 index 000000000..8ea52aab4 --- /dev/null +++ b/gr-radar/src/lib/Makefile.am @@ -0,0 +1,118 @@ +# +# Copyright 2004,2005 Free Software Foundation, Inc. +# +# This file is part of GNU Radio +# +# GNU Radio is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2, or (at your option) +# any later version. +# +# GNU Radio is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with GNU Radio; see the file COPYING. If not, write to +# the Free Software Foundation, Inc., 59 Temple Place - Suite 330, +# Boston, MA 02111-1307, USA. +# + +include $(top_srcdir)/Makefile.common + +LIBS += $(GNURADIO_CORE_LIBS) + +# Install this stuff so that it ends up as the gnuradio.radar module +# This usually ends up at: +# ${prefix}/lib/python${python_version}/site-packages/gnuradio + +ourpythondir = $(grpythondir) +ourlibdir = $(grpyexecdir) + +INCLUDES = $(STD_DEFINES_AND_INCLUDES) $(PYTHON_CPPFLAGS) + +SWIGPYTHONARGS = $(SWIGPYTHONFLAGS) $(STD_DEFINES_AND_INCLUDES) + +ALL_IFILES = \ + $(LOCAL_IFILES) \ + $(NON_LOCAL_IFILES) + +NON_LOCAL_IFILES = \ + $(top_srcdr)/gnuradio-core/src/lib/swig/gnuradio.i + + +LOCAL_IFILES = + + +# These files are built by SWIG. The first is the C++ glue. +# The second is the python wrapper that loads the _howto shared library +# and knows how to call our extensions. + +BUILT_SOURCES = + +lib_LTLIBRARIES = libradar.la + +libradar_la_SOURCES = \ + time_series.h \ + time_series.cc \ + simulation.h \ + simulation.cc + +bin_PROGRAMS = \ + xambi \ + eb-xambi \ + sim-airplane \ + sim-airplane2 + +xambi_SOURCES = xambi.cc +xambi_LDADD = libradar.la + +eb_xambi_SOURCES = eb-xambi.cc +eb_xambi_LDADD = libradar.la + +sim_airplane_SOURCES = sim-airplane.cc +sim_airplane_LDADD = libradar.la + +sim_airplane2_SOURCES = sim-airplane2.cc +sim_airplane2_LDADD = libradar.la + +# This gets howto.py installed in the right place +# ourpython_PYTHON = + +#ourlib_LTLIBRARIES = _howto.la + +# These are the source files that go into the shared library +#_howto_la_SOURCES = \ +# howto.cc \ +# howto_square_ff.cc \ +# howto_square2_ff.cc + +# magic flags +#_howto_la_LDFLAGS = $(NO_UNDEFINED) -module -avoid-version + +# link the library against some comon swig runtime code and the +# c++ standard library +#_howto_la_LIBADD = \ +# $(PYTHON_LDFLAGS) \ +# -lstdc++ + +#howto.cc howto.py: howto.i $(ALL_IFILES) +# $(SWIG) $(SWIGPYTHONARGS) -module howto -o howto.cc $< + +# These headers get installed in ${prefix}/include/gnuradio +grinclude_HEADERS = + + +# These swig headers get installed in ${prefix}/include/gnuradio/swig +swiginclude_HEADERS = \ + $(LOCAL_IFILES) + + +MOSTLYCLEANFILES = $(BUILT_SOURCES) *.pyc + + +# Don't distribute output of swig +dist-hook: + @for file in $(BUILT_SOURCES); do echo $(RM) $(distdir)/$$file; done + @for file in $(BUILT_SOURCES); do $(RM) $(distdir)/$$file; done diff --git a/gr-radar/src/lib/eb-xambi.cc b/gr-radar/src/lib/eb-xambi.cc new file mode 100644 index 000000000..c3c01953e --- /dev/null +++ b/gr-radar/src/lib/eb-xambi.cc @@ -0,0 +1,291 @@ +/* -*- c++ -*- */ +/* + * Copyright 2005 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2, or (at your option) + * any later version. + * + * GNU Radio is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ + +#ifndef _GNU_SOURCE +#define _GNU_SOURCE +#endif + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "time_series.h" + + +gr_complex +complex_conj_dotprod(const gr_complex *a, const gr_complex *b, unsigned n) +{ + gr_complex acc = 0; + for (unsigned i = 0; i < n; i++) + acc += a[i] * conj(b[i]); + + return acc; +} + +/*! + * \brief frequency translate src by normalized_freq + * + * \param dst destination + * \param src source + * \param n length of src and dst in samples + * \param normalized_freq [-1/2, +1/2] + */ +void +freq_xlate(gr_complex *dst, const gr_complex *src, unsigned n, float normalized_freq) +{ + gr_fxpt_nco nco; + nco.set_freq(2 * M_PI * normalized_freq); + + for (unsigned int i = 0; i < n; i++){ + gr_complex phasor(nco.cos(), nco.sin()); + dst[i] = src[i] * phasor; + nco.step(); + } +} + +inline void +write_float(FILE *output, float x) +{ + fwrite(&x, sizeof(x), 1, output); +} + + +// write 8-float header +static void +write_header (FILE *output, float ndoppler_bins, float min_doppler, float max_doppler) +{ + write_float(output, ndoppler_bins); + write_float(output, min_doppler); + write_float(output, max_doppler); + write_float(output, 0); + write_float(output, 0); + write_float(output, 0); + write_float(output, 0); + write_float(output, 0); +} + + +void +main_loop(FILE *output, time_series &ref_ts, time_series &scat0_ts, + unsigned nranges, unsigned correlation_window_size, + float min_doppler, float max_doppler, int ndoppler_bins) +{ + fprintf(stderr, "ndoppler_bins = %10d\n", ndoppler_bins); + fprintf(stderr, "min_doppler = %10f\n", min_doppler); + fprintf(stderr, "max_doppler = %10f\n", max_doppler); + + // float scale_factor = 1.0/correlation_window_size; // FIXME, not sure this is right + float scale_factor = 1.0; // FIXME, not sure this is right + + boost::scoped_array shifted(new gr_complex[correlation_window_size]); + + // gr_complex shifted[correlation_window_size]; // doppler shifted reference + + float doppler_incr = 0; + if (ndoppler_bins == 1){ + min_doppler = 0; + max_doppler = 0; + } + else + doppler_incr = (max_doppler - min_doppler) / (ndoppler_bins - 1); + + write_header(output, ndoppler_bins, min_doppler, max_doppler); + + unsigned long long ro = 0; // reference offset + unsigned long long so = 0; // scatter offset + + for (unsigned int n = 0; n < nranges; n++){ + if (0){ + fprintf(stdout, "n = %6d\n", n); + fprintf(stdout, "ro = %6lld\n", ro); + fprintf(stdout, "so = %6lld\n", so); + } + const gr_complex *ref = (const gr_complex *) ref_ts.seek(ro, correlation_window_size); + const gr_complex *scat0 = (const gr_complex *) scat0_ts.seek(so, correlation_window_size); + if (ref == 0 || scat0 == 0) + return; + + for (int nd = 0; nd < ndoppler_bins; nd++){ + float fdop = min_doppler + doppler_incr * nd; + //fprintf(stderr, "fdop = %10g\n", fdop); + freq_xlate(&shifted[0], ref, correlation_window_size, fdop); // generated doppler shifted reference + + gr_complex ccor = complex_conj_dotprod(&shifted[0], scat0, correlation_window_size); + float out = norm(ccor) * scale_factor; + + // fprintf(output, "%12g\n", out); + write_float(output, out); + } + + so += 1; + } +} + +static void +usage(const char *argv0) +{ + const char *progname; + const char *t = std::strrchr(argv0, '/'); + if (t != 0) + progname = t + 1; + else + progname = argv0; + + fprintf(stderr, "usage: %s [options] ref_file scatter_file\n", progname); + fprintf(stderr, " -o OUTPUT_FILENAME [default=eb-xambi.out]\n"); + fprintf(stderr, " -m MIN_RANGE [default=0]\n"); + fprintf(stderr, " -M MAX_RANGE [default=300]\n"); + fprintf(stderr, " -w CORRELATION_WINDOW_SIZE [default=2500]\n"); + fprintf(stderr, " -s NSAMPLES_TO_SKIP [default=0]\n"); + fprintf(stderr, " -d max_doppler (normalized: [0, +1/2)) [default=.0012]\n"); + fprintf(stderr, " -n ndoppler_bins [default=31]\n"); +} + +int +main(int argc, char **argv) +{ + int ch; + int min_range = 0; + int max_range = 300; + const char *ref_filename = 0; + const char *scatter_filename = 0; + const char *output_filename = "eb-xambi.out"; + unsigned int correlation_window_size = 2500; + long long int nsamples_to_skip = 0; + double max_doppler = 0.0012; + int ndoppler_bins = 31; + + + while ((ch = getopt(argc, argv, "m:M:ho:w:s:d:n:")) != -1){ + switch (ch){ + case 'm': + min_range = strtol(optarg, 0, 0); + break; + + case 'M': + max_range = strtol(optarg, 0, 0); + break; + + case 'w': + correlation_window_size = strtol(optarg, 0, 0); + if (correlation_window_size <= 1){ + usage(argv[0]); + fprintf(stderr, " correlation_window_size must be >= 1\n"); + exit(1); + } + break; + + case 'o': + output_filename = optarg; + break; + + case 's': + nsamples_to_skip = (long long) strtof(optarg, 0); + if (nsamples_to_skip < 0){ + usage(argv[0]); + fprintf(stderr, " nsamples_to_skip must be >= 0\n"); + exit(1); + } + break; + + case 'd': + max_doppler = strtof(optarg, 0); + if (max_doppler < 0 || max_doppler >= 0.5){ + usage(argv[0]); + fprintf(stderr, " max_doppler must be in [0, 0.5)\n"); + exit(1); + } + break; + + case 'n': + ndoppler_bins = strtol(optarg, 0, 0); + if (ndoppler_bins < 1){ + usage(argv[0]); + fprintf(stderr, " ndoppler_bins must >= 1\n"); + exit(1); + } + break; + + case '?': + case 'h': + default: + usage(argv[0]); + exit(1); + } + } // while getopt + + if (argc - optind != 2){ + usage(argv[0]); + exit(1); + } + + if (max_range < min_range){ + usage(argv[0]); + fprintf(stderr, " max_range must be >= min_range\n"); + exit(1); + } + unsigned int nranges = max_range - min_range + 1; + + ref_filename = argv[optind++]; + scatter_filename = argv[optind++]; + + FILE *output = fopen(output_filename, "wb"); + if (output == 0){ + perror(output_filename); + exit(1); + } + + unsigned long long ref_starting_offset = 0; + unsigned long long scatter_starting_offset = 0; + + if (min_range < 0){ + ref_starting_offset = -min_range; + scatter_starting_offset = 0; + } + else { + ref_starting_offset = 0; + scatter_starting_offset = min_range; + } + + ref_starting_offset += nsamples_to_skip; + scatter_starting_offset += nsamples_to_skip; + + try { + time_series ref(sizeof(gr_complex), ref_filename, ref_starting_offset, 0); + time_series scat0(sizeof(gr_complex), scatter_filename, scatter_starting_offset, 0); + + main_loop(output, ref, scat0, nranges, correlation_window_size, + -max_doppler, max_doppler, ndoppler_bins); + } + catch (std::string &s){ + std::cerr << s << std::endl; + exit(1); + } + + return 0; +} + diff --git a/gr-radar/src/lib/gen_run b/gr-radar/src/lib/gen_run new file mode 100755 index 000000000..8b9549e7b --- /dev/null +++ b/gr-radar/src/lib/gen_run @@ -0,0 +1,26 @@ +#!/usr/bin/env python + +import sys + +def main(): + window = 25000 + #window = 250000 + max_range = 99 + gifs = [] + for t in range(0, 241, 5): + sim_fn = "sim.%04d" % (t,) + xref_fn = "xref.%04d" % (t,) + ppm_fn = "i.%04d.ppm" % (t,) + gif_fn = "i.%04d.gif" % (t,) + gifs.append(gif_fn) + + sys.stdout.write("./sim-airplane2 rfm -o %s -S %d\n" % (sim_fn, t)) + sys.stdout.write("./eb-xambi -o %s -M %d -w %d -d 0.0002 rfm %s\n" % ( + xref_fn, max_range, window, sim_fn)) + sys.stdout.write("./plot_xambi_tool.m %s %s\n" % (xref_fn, ppm_fn)) + sys.stdout.write("convert -resize 800%% %s %s\n" % (ppm_fn, gif_fn)) + + sys.stdout.write("gifsicle -d 50 %s >animated.gif\n" % (' '.join(gifs))) + +if __name__ == '__main__': + main() diff --git a/gr-radar/src/lib/plot_rd.m b/gr-radar/src/lib/plot_rd.m new file mode 100644 index 000000000..166da020e --- /dev/null +++ b/gr-radar/src/lib/plot_rd.m @@ -0,0 +1,38 @@ +filename="xambi.out" + +ambdata=read_float_binary(filename); + +fftsize=512 +nranges = prod(size(ambdata))/fftsize + +ambdata=reshape(ambdata,fftsize,nranges); +% ambdata=fftshift(ambdata,1); + +% colormap(cool(256)); +colormap(rainbow(256)); + +d = 10*log10(ambdata); +min(min(d)) +max(max(d)) +imagesc(d, 2.0); + + + + + + + + + + + + + + + + + + + + + diff --git a/gr-radar/src/lib/plot_xambi.m b/gr-radar/src/lib/plot_xambi.m new file mode 100755 index 000000000..62088597b --- /dev/null +++ b/gr-radar/src/lib/plot_xambi.m @@ -0,0 +1,15 @@ +function plot_xambi(amb) + [nr, nc] = size(amb); + ndoppler_bins = nr + nranges = nc + + %colormap(cool(256)); + %colormap(rainbow(256)); + colormap(gray(1024)); + + %d = 10*log10(amb); + d = amb; + min(min(d)) + max(max(d)) + imagesc(d, 4.0); +endfunction; diff --git a/gr-radar/src/lib/plot_xambi_to_file.m b/gr-radar/src/lib/plot_xambi_to_file.m new file mode 100755 index 000000000..973d112a1 --- /dev/null +++ b/gr-radar/src/lib/plot_xambi_to_file.m @@ -0,0 +1,17 @@ +function plot_xambi_to_file(amb, filename) + [nr, nc] = size(amb); + ndoppler_bins = nr + nranges = nc + + %colormap(cool(256)); + %colormap(rainbow(256)); + colormap(gray(1024)); + + %d = 10*log10(amb); + d = amb; + min(min(d)) + max(max(d)) + b = imagesc(d, 4.0); + saveimage(filename, b, "ppm") + +endfunction; diff --git a/gr-radar/src/lib/plot_xambi_tool.m b/gr-radar/src/lib/plot_xambi_tool.m new file mode 100755 index 000000000..508db2866 --- /dev/null +++ b/gr-radar/src/lib/plot_xambi_tool.m @@ -0,0 +1,8 @@ +#!/usr/bin/env octave + +xambi_filename = argv{1}; +ppm_filename = argv{2}; + +a = read_xambi(xambi_filename); +plot_xambi_to_file(a, ppm_filename); + diff --git a/gr-radar/src/lib/sim-airplane.cc b/gr-radar/src/lib/sim-airplane.cc new file mode 100644 index 000000000..9b27ff009 --- /dev/null +++ b/gr-radar/src/lib/sim-airplane.cc @@ -0,0 +1,296 @@ +/* -*- c++ -*- */ +/* + * Copyright 2005 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2, or (at your option) + * any later version. + * + * GNU Radio is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ + +#ifndef _GNU_SOURCE +#define _GNU_SOURCE +#endif + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "time_series.h" +#include "simulation.h" + +static const double C = 3e8; // sped of light, m/s + + +// ------------------------------------------------------------------------ + +class delay_line { + std::vector d_z; + const int d_mask; + int d_newest; +public: + delay_line(unsigned int max_delay) + : d_z(gr_rounduppow2(max_delay)), d_mask(d_z.size()-1), d_newest(0) + { + } + + void + push_item(gr_complex x) + { + d_newest = (d_newest - 1) & d_mask; + d_z[d_newest] = x; + } + + gr_complex + ref_item(int delay) const + { + return d_z[(d_newest + delay) & d_mask]; + } +}; + +// ------------------------------------------------------------------------ + +class my_sim : public simulation +{ + FILE *d_output; + time_series &d_ref; + unsigned long long d_pos; // position in time series + delay_line d_z; + dyn_object *d_tx; // transmitter (not moving) + dyn_object *d_rx0; // receiver (not moving) + dyn_object *d_ac0; // aircraft (linear motion) + gr_fxpt_nco d_nco0; + + double d_baseline; // length of baseline in meters + double d_last_slant_range; + double d_range_bin; // meters/range_bin + float d_tx_lambda; // wavelength of tx signals in meters + float d_sample_rate; + float d_gain; // linear scale factor + +public: + my_sim(FILE *output, time_series &ref, double timestep, float sample_rate, + double tx_freq, double gain_db); + ~my_sim(); + + bool update(); + bool run(long long nsteps); + + bool write_output(gr_complex x) + { + return fwrite(&x, sizeof(x), 1, d_output) == 1; + } +}; + +my_sim::my_sim(FILE *output, time_series &ref, double timestep, + float sample_rate, double tx_freq, double gain_db) + : simulation(timestep), + d_output(output), d_ref(ref), d_pos(0), d_z(1024), + d_range_bin(C * timestep), d_tx_lambda(C/tx_freq), + d_sample_rate(sample_rate), d_gain(exp10(gain_db/10)) +{ + d_tx = new dyn_object(point(0,0), point(0,0), "Tx"); + d_rx0 = new dyn_object(point(45e3,0), point(0,0), "Rx0"); + + //float aircraft_speed = 135; // meters/sec (~ 300 miles/hr) + float aircraft_speed = 350; // meters/sec (~ 750 miles/hr) + float aircraft_angle = 250 * M_PI/180; + //point aircraft_pos = point(55e3, 20e3); + point aircraft_pos = point(55e3-5.54e3, 20e3-15.23e3); + d_ac0 = new dyn_object(aircraft_pos, + point(aircraft_speed * cos(aircraft_angle), + aircraft_speed * sin(aircraft_angle)), + "Ac0"); + add_object(d_tx); + add_object(d_rx0); + add_object(d_ac0); + + d_baseline = dyn_object::distance(*d_tx, *d_rx0); + d_last_slant_range = + dyn_object::distance(*d_tx, *d_ac0) + dyn_object::distance(*d_ac0, *d_rx0); +} + +my_sim::~my_sim() +{ +} + +bool +my_sim::update() +{ + // std::cout << *d_ac0 << std::endl; + + // compute slant_range and slant_range' + double slant_range = + dyn_object::distance(*d_tx, *d_ac0) + dyn_object::distance(*d_ac0, *d_rx0); // meters + double delta_slant_range = slant_range - d_last_slant_range; + d_last_slant_range = slant_range; + double deriv_slant_range_wrt_time = delta_slant_range / timestep(); // m/sec + + // fprintf(stdout, "%10.3f\t%10.3f\n", slant_range, deriv_slant_range_wrt_time); + + // grab new item from input and insert it into delay line + const gr_complex *in = (const gr_complex *) d_ref.seek(d_pos++, 1); + if (in == 0) + return false; + d_z.push_item(*in); + + // FIXME, may want to interpolate between two bins. + int int_delay = lrint((slant_range - d_baseline) / d_range_bin); + + gr_complex x = d_z.ref_item(int_delay); + + x = x * d_gain; // scale amplitude (this includes everything: RCS, antenna gain, losses, etc...) + + // compute doppler and apply it + float f_doppler = -deriv_slant_range_wrt_time / d_tx_lambda; + fprintf(stdout, "f_dop: %10.3f\n", f_doppler); + + d_nco0.set_freq(f_doppler / d_sample_rate); + gr_complex phasor(d_nco0.cos(), d_nco0.sin()); + // x = x * phasor; + d_nco0.step(); + + write_output(x); + + return simulation::update(); // run generic update +} + +bool +my_sim::run(long long nsteps) +{ + //fprintf(stdout, "<%12.2f, %12.2f>\n", d_ac0->pos().x(), d_ac0->pos().y()); + //std::cout << *d_ac0 << std::endl; + bool ok = simulation::run(nsteps); + //std::cout << *d_ac0 << std::endl; + //fprintf(stdout, "<%12.2f, %12.2f>\n", d_ac0->pos().x(), d_ac0->pos().y()); + return ok; +} + +// ------------------------------------------------------------------------ + +static void +usage(const char *argv0) +{ + const char *progname; + const char *t = std::strrchr(argv0, '/'); + if (t != 0) + progname = t + 1; + else + progname = argv0; + + fprintf(stderr, "usage: %s [options] ref_file\n", progname); + fprintf(stderr, " -o OUTPUT_FILENAME [default=sim.dat]\n"); + fprintf(stderr, " -n NSAMPLES_TO_PRODUCE [default=+inf]\n"); + fprintf(stderr, " -s NSAMPLES_TO_SKIP [default=0]\n"); + fprintf(stderr, " -g reflection gain in dB (should be <= 0) [default=0]\n"); + fprintf(stderr, " -f transmitter freq in Hz [default=100MHz]\n"); + fprintf(stderr, " -r sample rate in Hz [default=250kHz]\n"); +} + +int +main(int argc, char **argv) +{ + int ch; + const char *output_filename = "sim.dat"; + const char *ref_filename = 0; + long long int nsamples_to_skip = 0; + long long int nsamples_to_produce = std::numeric_limits::max(); + double sample_rate = 250e3; + double gain_db = 0; + double tx_freq = 100e6; + + while ((ch = getopt(argc, argv, "o:s:n:g:f:")) != -1){ + switch (ch){ + case 'o': + output_filename = optarg; + break; + + case 's': + nsamples_to_skip = (long long) strtof(optarg, 0); + if (nsamples_to_skip < 0){ + usage(argv[0]); + fprintf(stderr, " nsamples_to_skip must be >= 0\n"); + exit(1); + } + break; + + case 'n': + nsamples_to_produce = (long long) strtof(optarg, 0); + if (nsamples_to_produce < 0){ + usage(argv[0]); + fprintf(stderr, " nsamples_to_produce must be >= 0\n"); + exit(1); + } + break; + + case 'g': + gain_db = strtof(optarg, 0); + break; + + case 'f': + tx_freq = strtof(optarg, 0); + break; + + case 'r': + sample_rate = strtof(optarg, 0); + break; + + case '?': + case 'h': + default: + usage(argv[0]); + exit(1); + } + } // while getopt + + if (argc - optind != 1){ + usage(argv[0]); + exit(1); + } + + ref_filename = argv[optind++]; + + double timestep = 1.0/sample_rate; + + + FILE *output = fopen(output_filename, "wb"); + if (output == 0){ + perror(output_filename); + exit(1); + } + + unsigned long long ref_starting_offset = 0; + ref_starting_offset += nsamples_to_skip; + + try { + time_series ref(sizeof(gr_complex), ref_filename, ref_starting_offset, 0); + + my_sim simulator(output, ref, timestep, sample_rate, tx_freq, gain_db); + simulator.run(nsamples_to_produce); + } + catch (std::string &s){ + std::cerr << s << std::endl; + exit(1); + } + + return 0; +} + diff --git a/gr-radar/src/lib/sim-airplane2.cc b/gr-radar/src/lib/sim-airplane2.cc new file mode 100644 index 000000000..a8948a7bb --- /dev/null +++ b/gr-radar/src/lib/sim-airplane2.cc @@ -0,0 +1,372 @@ +/* -*- c++ -*- */ +/* + * Copyright 2005 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2, or (at your option) + * any later version. + * + * GNU Radio is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ + +#ifndef _GNU_SOURCE +#define _GNU_SOURCE +#endif + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "time_series.h" +#include "simulation.h" + +static const double C = 3e8; // sped of light, m/s + + +// ------------------------------------------------------------------------ + +class delay_line { + std::vector d_z; + const int d_mask; + int d_newest; +public: + delay_line(unsigned int max_delay) + : d_z(gr_rounduppow2(max_delay)), d_mask(d_z.size()-1), d_newest(0) + { + } + + void + push_item(gr_complex x) + { + d_newest = (d_newest - 1) & d_mask; + d_z[d_newest] = x; + } + + gr_complex + ref_item(int delay) const + { + return d_z[(d_newest + delay) & d_mask]; + } +}; + +// ------------------------------------------------------------------------ + +class aux_state { +public: + dyn_object *d_obj; + double d_last_slant_range; + gr_fxpt_nco d_nco; + + aux_state(dyn_object *obj) : d_obj(obj) {} +}; + +// ------------------------------------------------------------------------ + +class my_sim : public simulation +{ + FILE *d_output; + time_series &d_ref; + unsigned long long d_pos; // position in time series + delay_line d_z; + dyn_object *d_tx; // transmitter (not moving) + dyn_object *d_rx0; // receiver (not moving) + std::vector d_target; + + double d_baseline; // length of baseline in meters + double d_range_bin; // meters/range_bin + float d_tx_lambda; // wavelength of tx signals in meters + float d_sample_rate; + float d_gain; // linear scale factor + + void adjust_for_start_time(double start_time); + + bool write_output(gr_complex x) + { + return fwrite(&x, sizeof(x), 1, d_output) == 1; + } + +public: + my_sim(FILE *output, time_series &ref, double timestep, float sample_rate, + double start_time, double tx_freq, double gain_db); + ~my_sim(); + + bool update(); + bool run(long long nsteps); +}; + + +my_sim::my_sim(FILE *output, time_series &ref, double timestep, + float sample_rate, double start_time, + double tx_freq, double gain_db) + : simulation(timestep), + d_output(output), d_ref(ref), d_pos(0), d_z(1024), + d_range_bin(C * timestep), d_tx_lambda(C/tx_freq), + d_sample_rate(sample_rate), d_gain(exp10(gain_db/10)) +{ + d_tx = new dyn_object(point(0,0), point(0,0), "Tx"); + d_rx0 = new dyn_object(point(45e3,0), point(0,0), "Rx0"); + + add_object(d_tx); + add_object(d_rx0); + d_baseline = dyn_object::distance(*d_tx, *d_rx0); + + { + // add targets + float aircraft_speed; + float aircraft_angle; + point aircraft_pos; + dyn_object *ac; + + // target 1 + aircraft_speed = 135; // m/s + aircraft_angle = 240 * M_PI/180; + aircraft_pos = point(55e3, 20e3); + + ac = new dyn_object(aircraft_pos, + point(aircraft_speed * cos(aircraft_angle), + aircraft_speed * sin(aircraft_angle)), + "Ac0"); + add_object(ac); + d_target.push_back(new aux_state(ac)); + + // target 2 + aircraft_speed = 350; // m/s + aircraft_angle = 0 * M_PI/180; + aircraft_pos = point(-20e3, 60e3); + + ac = new dyn_object(aircraft_pos, + point(aircraft_speed * cos(aircraft_angle), + aircraft_speed * sin(aircraft_angle)), + "Ac1"); + add_object(ac); + d_target.push_back(new aux_state(ac)); + } + + adjust_for_start_time(start_time); + + for (unsigned i = 0; i < d_target.size(); i++) + d_target[i]->d_last_slant_range = + (dyn_object::distance(*d_tx, *d_target[i]->d_obj) + + dyn_object::distance(*d_target[i]->d_obj, *d_rx0)); + +} + +my_sim::~my_sim() +{ +} + +void +my_sim::adjust_for_start_time(double start_time) +{ + for (unsigned i = 0; i < d_obj.size(); i++){ + // Adjust initial starting positions depending on simulation + // start time. FIXME Assumes velocity is constant + point p = d_obj[i]->pos(); + point v = d_obj[i]->vel(); + p.set_x(p.x() + v.x() * start_time); + p.set_y(p.y() + v.y() * start_time); + d_obj[i]->set_pos(p); + } +} + +bool +my_sim::update() +{ + // std::cout << *d_ac0 << std::endl; + + // grab new item from input and insert it into delay line + const gr_complex *in = (const gr_complex *) d_ref.seek(d_pos++, 1); + if (in == 0) + return false; + d_z.push_item(*in); + + gr_complex s = 0; // output sample + // FIXME ought to add in attenuated direct path input + + + // for each target, compute slant_range and slant_range' + + for (unsigned i = 0; i < d_target.size(); i++){ + aux_state *t = d_target[i]; + + double slant_range = + (dyn_object::distance(*d_tx, *t->d_obj) + + dyn_object::distance(*t->d_obj, *d_rx0)); // meters + + double delta_slant_range = slant_range - t->d_last_slant_range; + t->d_last_slant_range = slant_range; + double deriv_slant_range_wrt_time = delta_slant_range / timestep(); // m/sec + + //fprintf(stdout, "%10.3f\t%10.3f\n", slant_range, deriv_slant_range_wrt_time); + + // FIXME, may want to interpolate between two bins. + int int_delay = lrint((slant_range - d_baseline) / d_range_bin); + + gr_complex x = d_z.ref_item(int_delay); + + // scale amplitude (this includes everything: RCS, antenna gain, losses, etc...) + x = x * d_gain; + + if (1){ + // compute doppler and apply it + float f_doppler = -deriv_slant_range_wrt_time / d_tx_lambda; + + t->d_nco.set_freq(f_doppler / d_sample_rate); + gr_complex phasor(t->d_nco.cos(), t->d_nco.sin()); + x = x * phasor; + t->d_nco.step(); + } + + s += x; // add in this target's contribution + } + + write_output(s); + + return simulation::update(); // run generic update +} + +bool +my_sim::run(long long nsteps) +{ + //fprintf(stdout, "<%12.2f, %12.2f>\n", d_ac0->pos().x(), d_ac0->pos().y()); + //std::cout << *d_ac0 << std::endl; + bool ok = simulation::run(nsteps); + //std::cout << *d_ac0 << std::endl; + //fprintf(stdout, "<%12.2f, %12.2f>\n", d_ac0->pos().x(), d_ac0->pos().y()); + return ok; +} + +// ------------------------------------------------------------------------ + +static void +usage(const char *argv0) +{ + const char *progname; + const char *t = std::strrchr(argv0, '/'); + if (t != 0) + progname = t + 1; + else + progname = argv0; + + fprintf(stderr, "usage: %s [options] ref_file\n", progname); + fprintf(stderr, " -o OUTPUT_FILENAME [default=sim.dat]\n"); + fprintf(stderr, " -n NSAMPLES_TO_PRODUCE [default=+inf]\n"); + fprintf(stderr, " -s NSAMPLES_TO_SKIP [default=0]\n"); + fprintf(stderr, " -g reflection gain in dB (should be <= 0) [default=0]\n"); + fprintf(stderr, " -f transmitter freq in Hz [default=100MHz]\n"); + fprintf(stderr, " -r sample rate in Hz [default=250kHz]\n"); + fprintf(stderr, " -S simulation start time in seconds [default=0]\n"); +} + +int +main(int argc, char **argv) +{ + int ch; + const char *output_filename = "sim.dat"; + const char *ref_filename = 0; + long long int nsamples_to_skip = 0; + long long int nsamples_to_produce = std::numeric_limits::max(); + double sample_rate = 250e3; + double gain_db = 0; + double tx_freq = 100e6; + double start_time = 0; + + while ((ch = getopt(argc, argv, "o:s:n:g:f:S:")) != -1){ + switch (ch){ + case 'o': + output_filename = optarg; + break; + + case 's': + nsamples_to_skip = (long long) strtof(optarg, 0); + if (nsamples_to_skip < 0){ + usage(argv[0]); + fprintf(stderr, " nsamples_to_skip must be >= 0\n"); + exit(1); + } + break; + + case 'n': + nsamples_to_produce = (long long) strtof(optarg, 0); + if (nsamples_to_produce < 0){ + usage(argv[0]); + fprintf(stderr, " nsamples_to_produce must be >= 0\n"); + exit(1); + } + break; + + case 'g': + gain_db = strtof(optarg, 0); + break; + + case 'f': + tx_freq = strtof(optarg, 0); + break; + + case 'r': + sample_rate = strtof(optarg, 0); + break; + + case 'S': + start_time = strtof(optarg, 0); + break; + + case '?': + case 'h': + default: + usage(argv[0]); + exit(1); + } + } // while getopt + + if (argc - optind != 1){ + usage(argv[0]); + exit(1); + } + + ref_filename = argv[optind++]; + + double timestep = 1.0/sample_rate; + + + FILE *output = fopen(output_filename, "wb"); + if (output == 0){ + perror(output_filename); + exit(1); + } + + unsigned long long ref_starting_offset = 0; + ref_starting_offset += nsamples_to_skip; + + try { + time_series ref(sizeof(gr_complex), ref_filename, ref_starting_offset, 0); + + my_sim simulator(output, ref, timestep, sample_rate, start_time, + tx_freq, gain_db); + simulator.run(nsamples_to_produce); + } + catch (std::string &s){ + std::cerr << s << std::endl; + exit(1); + } + + return 0; +} + diff --git a/gr-radar/src/lib/simulation.cc b/gr-radar/src/lib/simulation.cc new file mode 100644 index 000000000..4ebc7f1f8 --- /dev/null +++ b/gr-radar/src/lib/simulation.cc @@ -0,0 +1,89 @@ +/* -*- c++ -*- */ +/* + * Copyright 2005 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2, or (at your option) + * any later version. + * + * GNU Radio is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include "simulation.h" + +bool +dyn_object::update(double delta_t) +{ + double new_x = pos().x() + delta_t * vel().x(); + double new_y = pos().y() + delta_t * vel().y(); + set_pos(point(new_x, new_y)); + return true; +} + +simulation::~simulation() +{ + for (unsigned i = 0; i < d_obj.size(); i++){ + delete d_obj[i]; + d_obj[i] = 0; + } +} + +bool +simulation::update() +{ + bool ok = true; + for (unsigned i = 0; i < d_obj.size(); i++){ + ok &= d_obj[i]->update(d_timestep); + } + d_now += d_timestep; + return ok; +} + +bool +simulation::run(long long nsteps) +{ + for (long long i = 0; i < nsteps; i++) + if (!update()) + return false; + + return true; +} + +void +simulation::add_object(dyn_object *obj) +{ + d_obj.push_back(obj); +} + +// ---------------------------------------------------------------- + +std::ostream& operator<<(std::ostream& out, const dyn_object& o) +{ + out << "<" << o.name() + << " pos: " << o.pos() + << " vel: " << o.vel() + << ">"; + return out; +} + +std::ostream& operator<<(std::ostream& out, const point& p) +{ + out << "(" << p.x() << ", " << p.y() << ")"; + return out; +} + diff --git a/gr-radar/src/lib/simulation.h b/gr-radar/src/lib/simulation.h new file mode 100644 index 000000000..ff32484ed --- /dev/null +++ b/gr-radar/src/lib/simulation.h @@ -0,0 +1,116 @@ +/* -*- c++ -*- */ +/* + * Copyright 2005 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2, or (at your option) + * any later version. + * + * GNU Radio is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ + +#ifndef INCLUDED_SIMULATION_H +#define INCLUDED_SIMULATION_H + +#include +#include +#include +#include + + +class point { + double d_x, d_y; +public: + point(double x = 0, double y = 0) : d_x(x), d_y(y) {} + + double x() const { return d_x; } + double y() const { return d_y; } + + void set_x(double x) { d_x = x; } + void set_y(double y) { d_y = y; } + + static point + add(const point &p1, const point &p2) + { + return point(p1.x() + p2.x(), + p1.y() + p2.y()); + } + + static point + sub(const point &p1, const point &p2) + { + return point(p1.x() - p2.x(), p1.y() - p2.y()); + } + + static double + distance(const point &p1, const point &p2) + { + point d = point::sub(p1, p2); + return std::sqrt(d.x()*d.x() + d.y()*d.y()); + } + +}; + + +class dyn_object { + point d_pos; + point d_vel; + std::string d_name; +public: + dyn_object(point pos=point(0,0), point vel=point(0,0), const std::string name="") + : d_pos(pos), d_vel(vel), d_name(name) {} + + virtual ~dyn_object() {} + + point pos() const { return d_pos; } + point vel() const { return d_vel; } + std::string name() const { return d_name; } + + void set_pos(point pos) { d_pos = pos; } + void set_vel(point vel) { d_vel = vel; } + + virtual bool update(double delta_t); + + static double + distance(const dyn_object &o1, const dyn_object &o2) + { + return point::distance(o1.pos(), o2.pos()); + } + +}; + + +class simulation { + double d_timestep; + double d_now; +protected: + std::vector d_obj; + +public: + simulation(double timestep = 1.0, double now = 0.0) + : d_timestep(timestep), d_now(now) {} + virtual ~simulation(); + virtual bool update(); + virtual bool run(long long nsteps); + + void add_object(dyn_object *obj); + double now() const { return d_now; } + double timestep() const { return d_timestep; } +}; + +std::ostream& operator<<(std::ostream& out, const dyn_object& o); +std::ostream& operator<<(std::ostream& out, const point& p); + +#endif /* INCLUDED_SIMULATION_H */ + diff --git a/gr-radar/src/lib/time_series.cc b/gr-radar/src/lib/time_series.cc new file mode 100644 index 000000000..2767165dc --- /dev/null +++ b/gr-radar/src/lib/time_series.cc @@ -0,0 +1,96 @@ +/* -*- c++ -*- */ +/* + * Copyright 2005 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2, or (at your option) + * any later version. + * + * GNU Radio is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ + +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include "time_series.h" +#include +#include +#include +#include +#include +#include + + +time_series::time_series(size_t itemsize, const std::string filename, + unsigned long long starting_offset, + long long nsamples_to_process) + : d_itemsize(itemsize), d_filename(filename), d_fd(-1), + d_start(starting_offset), d_buffer(0) +{ + if ((d_fd = open(d_filename.c_str(), O_RDONLY | O_LARGEFILE, 0660)) == -1){ + perror(d_filename.c_str()); + throw std::string("open failed: ") + d_filename; + } + + struct stat statbuf; + if (fstat(d_fd, &statbuf) == -1){ + perror(d_filename.c_str()); + throw std::string("fstat failed: ") + d_filename; + } + d_filesize = statbuf.st_size; + d_limit = d_filesize / d_itemsize; + + if (d_start > d_limit){ + std::string s = std::string("d_start > filesize: ") + d_filename; + std::cerr << s + << " d_start " << d_start + << " d_limit " << d_limit << std::endl; + throw s; + } + + if (nsamples_to_process > 0) + if ((d_start + nsamples_to_process) < d_limit) + d_limit = d_start + nsamples_to_process; + + d_buffer = mmap(0, d_filesize, PROT_READ, MAP_SHARED, d_fd, 0); + if (d_buffer == MAP_FAILED){ + perror("mmap"); + throw std::string("mmap failed: ") + d_filename; + } +} + +time_series::~time_series() +{ + munmap(d_buffer, d_filesize); + close(d_fd); +} + +const void * +time_series::seek(unsigned long long pos, unsigned long long blocksize) const +{ + if ((d_start + pos + blocksize) >= d_limit) + return 0; + + return (const void *)((char *)d_buffer + ((d_start + pos) * d_itemsize)); +} + +long long +time_series::nsamples_available(unsigned long long pos) const +{ + if ((d_start + pos) >= d_limit) + return 0; + + return d_limit - (d_start + pos); +} diff --git a/gr-radar/src/lib/time_series.h b/gr-radar/src/lib/time_series.h new file mode 100644 index 000000000..3764ba50d --- /dev/null +++ b/gr-radar/src/lib/time_series.h @@ -0,0 +1,73 @@ +/* -*- c++ -*- */ +/* + * Copyright 2005 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2, or (at your option) + * any later version. + * + * GNU Radio is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ +#ifndef INCLUDED_TIME_SERIES_H +#define INCLUDED_TIME_SERIES_H + +#include + +/*! + * \brief Memory mapped input for complex time series data + * + * Inspired by "iqts" by John Sahr, Univ of Washington + */ + +class time_series { + size_t d_itemsize; // user specified item size + std::string d_filename; + int d_fd; // file descriptor + unsigned long long d_filesize; // in bytes + unsigned long long d_start; // in items + unsigned long long d_limit; // in items + void *d_buffer; // points to base of file + +public: + /*! + * \brief Create read-only mapped file accessor. + * \param item_size size of item in bytes + * \param filename name of file to open + * \param starting_offset offset in file in item_size units at which to start + * \param nsamples_to_process maximum number of samples to map in starting at \p start. -1 implies no limit. + * + * \throws string on error opening file, etc. + */ + time_series(size_t item_size, const std::string filename, + unsigned long long starting_offset=0, + long long nsamples_to_process=-1); + ~time_series(); + + /*! + * \brief Return a pointer to a buffer of data at file offset pos. + * + * \param pos offset from beginning of file in itemsize units. + * \param blocksize minimum size of returned buffer in itemsize units. + * + * "Seek" to pos in file and return a pointer to the data at that + * location. The returned pointer will have at least blocksize valid + * elements. Return 0 if pos is out of bounds, or if there isn't + * at least blocksize units available in the file. + */ + const void *seek(unsigned long long pos, unsigned long long blocksize) const; + + long long nsamples_available(unsigned long long pos) const; +}; + +#endif /* INCLUDED_TIME_SERIES_H */ diff --git a/gr-radar/src/lib/xambi.cc b/gr-radar/src/lib/xambi.cc new file mode 100644 index 000000000..e5d4a9a95 --- /dev/null +++ b/gr-radar/src/lib/xambi.cc @@ -0,0 +1,259 @@ +/* -*- c++ -*- */ +/* + * Copyright 2006 Free Software Foundation, Inc. + * + * This file is part of GNU Radio + * + * GNU Radio is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2, or (at your option) + * any later version. + * + * GNU Radio is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with GNU Radio; see the file COPYING. If not, write to + * the Free Software Foundation, Inc., 59 Temple Place - Suite 330, + * Boston, MA 02111-1307, USA. + */ + +#include +#include +#include +#include + +#include +#include +#include "time_series.h" + + +using namespace std; + +/*! + * \file xambi.cc driver for computation of cross ambiguity + * + * Based on ideas liberally lifted from a version of xambi.cc + * obtained from John Sahr, that identified these people as authors: + * John Sahr jdsahr@u.washington.edu + * Frank Lind flind@haystack.mit.edu + * Chucai "Cliff" Zhou + * Melissa Meyer mgmeyer@ee.washington.edu + * + * Extensively revised since then. + */ + +static int default_decimation = 40; +static int default_fftsize = 256; +static int default_naverages = 1000000; // infinite + +/// usage for the xambi program + +static void usage(char *argv0) +{ + cerr << "usage: xambi [opts] scatterfile" << endl; + cerr << " where [opts] are" << endl; + cerr << " -x reffilename : Required" << endl; + cerr << " -o outputfilename[=xambi.out] : create or overwrite" << endl; + cerr << " -a averages[=100000] : restart accumulation; 100000 = infinite" << endl; + cerr << " -d decimationfactor[=40]" << endl; + cerr << " -f fftsize[=256]" << endl; + cerr << " -r range[=0] : first range" << endl; + cerr << " -n nranges[=1] : range count" << endl; + cerr << " -S start[=0] : starting offset in file" << endl; + cerr << " -N nsamples[=inf] : # of samples to process" << endl; + cerr << " -v : increment verbosity" << endl; + cerr << " -h : be helpful" << endl; + cerr << endl; + cerr << "The reffile and scatterfile are native-endian binary complex" << endl; + cerr << "and must be sampled at the same rate." << endl; + + exit(0); +} + +gr_complex +complex_dot_product(const gr_complex *xx, const gr_complex *yy, int nterms) +{ + gr_complex sum(0); + + for (int i = 0; i < nterms; i++) + sum += xx[i] * yy[i]; // FIXME? conj(yy[i]) + + return sum; +} + + +/// the main driver + +int +main(int argc, char *argv[]) +{ + char *ref_fname = 0; //< holds name of reference signal source + char *out_fname = 0; //< holds name of processed output + int decimation = default_decimation; + int range = 0; //< first range of range block + int nranges = 1; //< number of ranges of range block + int fftsize = default_fftsize; + int naverages = default_naverages; + int verbosity = 0; + int blocksize = 0; + int offset = 0; + unsigned long long starting_file_offset = 0; + unsigned long long nsamples_to_process = (unsigned long long) -1; + + int f; + + const gr_complex *x, *y; + const gr_complex *xx, *yy; + + int c; // used to process the command line + int r; // an index to count over ranges + int i; // an index to count through the time series + int a; + + while((c = getopt(argc,argv,"a:o:x:y:r:d:f:n:hvS:N:")) != EOF) { + switch(c) { + case 'x': ref_fname = optarg; break; + case 'o': out_fname = optarg; break; + case 'a': naverages = atoi(optarg); break; + case 'd': decimation = atoi(optarg); break; + case 'r': range = atoi(optarg); break; + case 'n': nranges = atoi(optarg); break; + case 'f': fftsize = atoi(optarg); break; + case 'S': starting_file_offset = strtoll(optarg, 0, 0); break; + case 'N': nsamples_to_process = strtoll(optarg, 0, 0); break; + case 'v': verbosity++; break; + case 'h': usage(argv[0]); break; + default: usage(argv[0]); break; + } + } + + // Wrapper for FFTW 1d forward FFT. N.B. output is not scaled by 1/fftsize + gri_fft_complex fft(fftsize, true); + gr_complex *fft_input = fft.get_inbuf(); + gr_complex *fft_output = fft.get_outbuf(); + + if(range < 0) { + cerr << "you specified -r " << range << "; must be non-negative (exit)" << endl; + exit(1); + } + + if(nranges < 1) { + cerr << "you specified -n " << nranges << "; must be positive (exit)" << endl; + exit(1); + } + + if(decimation < 1) { + cerr << "you specified -d " << decimation << "; must be positive (exit)" << endl; + exit(1); + } + + if(naverages < 1) { + cerr << "you specified -a " << naverages << "; must be positive (exit)" << endl; + exit(1); + } + + if(ref_fname == 0) { + cerr << "you must specify a reference signal with the -x option" << endl; + usage(argv[0]); + } + + if (optind >= argc) { + cerr << "you must specify a scattering signal after all other options" << endl; + usage(argv[0]); + } + + time_series X(sizeof(gr_complex), ref_fname, + starting_file_offset, nsamples_to_process); + + time_series Y(sizeof(gr_complex), argv[optind], + starting_file_offset, nsamples_to_process); // add more for interferometry ... + + float psd[fftsize*nranges]; + + if(out_fname == 0) { + char fname[200]; + snprintf(fname, sizeof(fname), "%s.out", "xambi"); + out_fname = fname; + } + + ofstream Z(out_fname); + + blocksize = fftsize*decimation + nranges; + offset = 0; + a = 0; + + // the fftsize is squared because we're using norm, not abs, + // when computing the psd + float scale_factor = 1.0 / (fftsize * fftsize); + + for(i = 0; i < nranges*fftsize; i++) + psd[i] = 0.0; + + + while(1){ // loop over data until exhausted. + if(verbosity > 1) { + cerr << " " << a; // write out the number of completed averages + cerr.flush(); + } + + x = (const gr_complex *) X.seek(offset, blocksize); + y = (const gr_complex *) Y.seek(offset + range, blocksize); + + if(!x || !y) // ran out of data; stop integrating + break; + + for(r = 0; r < nranges; r++) { // For Each Range ... + xx = x; + yy = y + r; + + for(f = 0; f < fftsize; f++) { // and for each Doppler bin ... + + // cross correlate and do a boxcar decimation + + fft_input[f] = complex_dot_product(xx, yy, decimation); + xx += decimation; + yy += decimation; + } + + fft.execute(); // input: fft_input; output: fft_output + + for(f = 0; f < fftsize; f++) { + psd[r*fftsize + f] += norm(fft_output[f]); + } + } // end range + + a++; + offset += fftsize * decimation; + + if(a >= naverages) { + if(verbosity > 0) + cerr << " dumping " << endl; + + for(i = 0; i < nranges*fftsize; i++) // normalize + psd[i] *= scale_factor; + + Z.write((const char *) psd, nranges*fftsize*sizeof(float)); + + for(i = 0; i < nranges*fftsize; i++) + psd[i] = 0.0; + + a = 0; + } + } + + if(verbosity > 1) + printf("\n"); + + if(a > 0) { + for(i = 0; i < nranges*fftsize; i++) // normalize + psd[i] *= scale_factor; + + Z.write((const char *) psd, nranges*fftsize*sizeof(float)); + } + + return 0; +} + -- cgit