updated comments

Makefile.am

@@ -4,7 +4,10 @@ include_HEADERS = harminv.h
 dist_man_MANS = harminv.1
 noinst_PROGRAMS = sines
 
-EXTRA_DIST = COPYRIGHT
+pkgconfigdir = $(libdir)/pkgconfig
+pkgconfig_DATA = harminv.pc
+
+EXTRA_DIST = COPYRIGHT harminv.pc.in
 
 libharminv_la_SOURCES = harminv.c harminv.h harminv-int.h check.h
 libharminv_la_LDFLAGS = -version-info @SHARED_VERSION_INFO@

README

@@ -1,5 +1,6 @@
 			       harminv:
 		  Harmonic Inversion of Time Signals
+		 by the Filter Diagonalizatino Method
 
 	       Steven G. Johnson <[email protected]>
 		Massachusetts Institute of Technology
@@ -9,24 +10,95 @@
 Introduction:
 =============
 
-harminv is a free program to solve the problem of "harmonic
-inversion."  Given a discrete, finite-length signal that consists of a
-sum of finitely-many sinusoids (possibly exponentially decaying), it
-determines the frequencies, decay constants, amplitudes, and phases of
-those sinusoids.
+harminv is a free program (and accompanying library) to solve the
+problem of "harmonic inversion."  Given a discrete, finite-length
+signal that consists of a sum of finitely-many sinusoids (possibly
+exponentially decaying), it determines the frequencies, decay
+constants, amplitudes, and phases of those sinusoids.
 
 It can, in principle, provide much better accuracy than
 straightforward FFT based methods, essentially because it assumes a
 specific form for the signal.  (Fourier transforms, in contrast,
-attempt to represent *any* data as a sum of sinusoidal component.)
+attempt to represent *any* data as a sum of sinusoidal components.)
 
 This kind of spectral analysis has wide applications in many areas of
-physics and engineering, and other fields.  For example, it could be
-used to extract the vibrational or "eigen" modes of a system from its
-response to some stimulus, and also their rates of decay in
-dissipative systems.
+physics and engineering, as well as other fields.  For example, it
+could be used to extract the vibrational or "eigen" modes of a system
+from its response to some stimulus, and also their rates of decay in
+dissipative systems.  This has been applied to analyze, e.g., NMR
+experimental data and numerical simulations of quantum mechanics.
 
-------------------------------------------------------------------------------
+We use a low-storage "filter diagonalization method" (FDM) for
+finding the sinusoids near a given frequency interval, described in:
+
+   V. A. Mandelshtam and H. S. Taylor, "Harmonic inversion of time
+   signals," J. Chem. Phys., vol. 107, no. 17, p. 6756-6769 (Nov. 1
+   1997).  See also erratum, ibid, vol. 109, no. 10, p. 4128 (Sep. 8
+   1998).
+
+(See also comments in harminv.c.)
+
+Program Usage:
+==============
+
+The usage of the harminv program is described by its man page ('man
+harminv'), included in the installation below.  To briefly summarize,
+it takes a sequence of numbers (real or complex) from standard input
+and a range of frequencies to search and outputs the frequencies it
+finds.
+
+Library Usage:
+==============
+
+The usage of the library -lharminv is analogous to the program.  In C
+or C++, you first #include <harminv.h>, then specify the data and the
+frequency range by calling harminv_data_create, returning a
+harminv_data data structure:
+
+harminv_data harminv_data_create(int n, const harminv_complex *signal,
+                                 double fmin, double fmax, int nf);
+
+signal is a pointer to an array of n complex numbers.  In C++,
+harminv_complex is std::complex<double>.  In C, harminv_complex is a
+double[2] with the real parts in signal[i][0] and the imaginary parts
+in signal[i][1].  fmin and fmax are the frequency range to search, and
+nf is an upper bound on the number of frequencies in that range (100
+is typically a safe value).  Frequencies are in units corresponding to
+a sampling interval of 1 time unit--if your actual sampling interval
+is dt, then you should rescale your frequencies by multiplying them by
+dt.
+
+Then, you solve for the frequencies by calling:
+
+void harminv_solve(harminv_data d);
+
+Then, the frequencies and other data can be extracted from d by the
+following routines.  The number N of frequencies found is returned by:
+
+int harminv_get_num_freqs(harminv_data d);
+
+Then, for each index 0 <= k < N, the corresponding frequency and decay
+constant (as defined in 'man harminv') are returned by:
+
+double harminv_get_freq(harminv_data d, int k);
+double harminv_get_decay(harminv_data d, int k);
+
+An array of the N corresponding complex amplitudes is allocated and
+returned by:
+
+harminv_complex *harminv_compute_amplitudes(harminv_data d);
+
+When you are done with this array, you should deallocate it by calling
+free().  Similarly, an array of N "error" figures of merit is returned
+by:
+
+double *harminv_compute_frequency_errors(harminv_data d);
+
+As described in 'man harminv', this figure of merit is not simply an
+error bar, and its quantitative interpretation seems difficult.  (See
+harminv.c for a reference to the method used to compute this figure of
+merit.)  However, solutions with a much larger error figure of merit
+compared to other solutions are likely to be spurious.
 
 Installation:
 =============

autogen.sh

@@ -0,0 +1,7 @@
+#!/bin/sh
+
+autoreconf --verbose --install --symlink --force
+autoreconf --verbose --install --symlink --force
+autoreconf --verbose --install --symlink --force
+
+./configure --enable-maintainer-mode

configure.ac

@@ -28,6 +28,9 @@ if test "$GCC" = "yes" && test "$USE_MAINTAINER_MODE" = yes; then
 	CFLAGS="$CFLAGS -Wall -W -Wbad-function-cast -Wcast-qual -Wpointer-arith -Wcast-align -pedantic"
 fi
 
+# Checks for header files.
+AC_CHECK_HEADERS(unistd.h getopt.h)
+
 ###########################################################################
 
 AC_ARG_WITH(cxx, [  --with-cxx=<dir>        force use of C++ and complex<double>], with_cxx=$withval, with_cxx=no)
@@ -99,9 +102,6 @@ fi
 AC_PROG_LIBTOOL                                                                                
 ###########################################################################
 
-# Checks for header files.
-AC_CHECK_HEADERS(unistd.h getopt.h)
-
 AC_CHECK_LIB(m, sqrt)
 
 ###########################################################################
@@ -117,5 +117,5 @@ LIBS="$LAPACK_LIBS $BLAS_LIBS $LIBS $FLIBS"
 
 ###########################################################################
 
-AC_CONFIG_FILES([Makefile])
+AC_CONFIG_FILES([Makefile harminv.pc])
 AC_OUTPUT

harminv.1

@@ -33,7 +33,7 @@ extracts their decay rates as well.
 robustness than Fourier-transform methods, essentially because it
 assumes a specific form for the input.
 
-It uses a "low-storage filter-diagonalization" (LSFD) algorithm, as
+It uses a low-storage "filter-diagonalization method" (FDM), as
 described in V. A. Mandelshtam and H. S. Taylor, "Harmonic inversion
 of time signals," \fIJ. Chem. Phys.\fR \fB107\fR, 6756 (1997).  See
 also erratum, \fIibid\fR \fB109\fR, 4128 (1998).

harminv.c

@@ -34,7 +34,7 @@
    the sinusoids much more accurately than we could via taking the FFT
    and looking at the peaks, for the same number of samples.
 
-   We use an "low-storage filter diagonalization" algorithm for finding
+   We use a low-storage "filter diagonalization method" (FDM) for finding
    the sinusoids near a given frequency interval, described in:
 
    V. A. Mandelshtam and H. S. Taylor, "Harmonic inversion of time
@@ -555,7 +555,7 @@ void harminv_solve(harminv_data d)
 
 /* Returns an array (of size harminv_get_num_freqs(d)) of estimates
    for the |error| in the solution frequencies.  Solutions with
-   errors not << 1 are likely to be spurious. */
+   errors much larger than the smallest error are likely to be spurious. */
 double *harminv_compute_frequency_errors(harminv_data d)
 {
      int i, J2, one = 1;

harminv.pc.in

@@ -0,0 +1,10 @@
+prefix=@prefix@
+exec_prefix=@exec_prefix@
+libdir=@libdir@
+includedir=@includedir@
+
+Name: Harminv
+Description: harmonic inversion of time series by filter diagonalization
+Version: @VERSION@
+Libs: -L${libdir} -lharminv @LIBS@
+Cflags: -I${includedir}

mkdist.sh

@@ -0,0 +1,15 @@
+#!/bin/sh
+
+cvs update -d
+
+# hackery to build ChangeLog
+echo >ChangeLog
+rcs2log -l 72 >ChangeLog
+emacs -batch -q -no-site-file --eval \
+   '(progn (find-file "ChangeLog") (fill-region 1 1000000) (save-buffer))'
+
+sh autogen.sh
+
+make maintainer-clean
+./configure --enable-maintainer-mode
+make dist