Houston, we have a trunk.

git-svn-id: http://gnuradio.org/svn/gnuradio/trunk@3122 221aa14e-8319-0410-a670-987f0aec2ac5

14 files changed, 1814 insertions, 0 deletions

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 <iostream>
+#include <string>
+#include <fstream>
+#include <unistd.h>
+#include <stdlib.h>
+#include <getopt.h>
+#include <boost/scoped_array.hpp>
+#include <gr_complex.h>
+#include <gr_fxpt_nco.h>
+#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<gr_complex>  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 <iostream>
+#include <string>
+#include <fstream>
+#include <unistd.h>
+#include <stdlib.h>
+#include <gr_complex.h>
+#include <getopt.h>
+#include <gr_misc.h>
+#include <limits>
+#include <gr_fxpt_nco.h>
+#include "time_series.h"
+#include "simulation.h"
+
+static const double C = 3e8;	// sped of light, m/s
+
+
+// ------------------------------------------------------------------------
+
+class delay_line {
+  std::vector<gr_complex>	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<long long int>::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 <iostream>
+#include <string>
+#include <fstream>
+#include <unistd.h>
+#include <stdlib.h>
+#include <gr_complex.h>
+#include <getopt.h>
+#include <gr_misc.h>
+#include <limits>
+#include <gr_fxpt_nco.h>
+#include "time_series.h"
+#include "simulation.h"
+
+static const double C = 3e8;	// sped of light, m/s
+
+
+// ------------------------------------------------------------------------
+
+class delay_line {
+  std::vector<gr_complex>	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<aux_state*> 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<long long int>::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 <string>
+#include <cmath>
+#include <vector>
+#include <iostream>
+
+
+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<dyn_object *>	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 <sys/mman.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <unistd.h>
+#include <iostream>
+
+
+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 <string>
+
+/*!
+ * \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 <iostream>
+#include <string>
+#include <fstream>
+#include <unistd.h>
+
+#include <gr_complex.h>
+#include <gri_fft.h>
+#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<float>" << 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;
+}
+