[xiph-commits] r17952 - trunk/chirptest

xiphmont at svn.xiph.org xiphmont at svn.xiph.org
Tue May 3 13:34:45 PDT 2011


Author: xiphmont
Date: 2011-05-03 13:34:45 -0700 (Tue, 03 May 2011)
New Revision: 17952

Added:
   trunk/chirptest/AUTHORS
   trunk/chirptest/COPYING
   trunk/chirptest/ChangeLog
   trunk/chirptest/INSTALL
   trunk/chirptest/Makefile.am
   trunk/chirptest/NEWS
   trunk/chirptest/README
   trunk/chirptest/TODO
   trunk/chirptest/autogen.sh
   trunk/chirptest/chirp.c
   trunk/chirptest/chirp.h
   trunk/chirptest/chirpgraph.c
   trunk/chirptest/chirpgraph.h
   trunk/chirptest/chirptest.c
   trunk/chirptest/configure.ac
   trunk/chirptest/scales.h
   trunk/chirptest/window.c
   trunk/chirptest/window.h
Log:
Add source files / build system


Copied: trunk/chirptest/AUTHORS (from rev 13251, trunk/ghost/AUTHORS)
===================================================================
--- trunk/chirptest/AUTHORS	                        (rev 0)
+++ trunk/chirptest/AUTHORS	2011-05-03 20:34:45 UTC (rev 17952)
@@ -0,0 +1,3 @@
+Benchmarking, graphing: Monty <monty at xiph.org>
+Original estimator code: Jean-Marc Valin <jm at xiph.org>
+Additional optimizations and windows: Greg Maxwell <greg at xiph.org>
\ No newline at end of file

Copied: trunk/chirptest/COPYING (from rev 13251, trunk/ghost/COPYING)
===================================================================
--- trunk/chirptest/COPYING	                        (rev 0)
+++ trunk/chirptest/COPYING	2011-05-03 20:34:45 UTC (rev 17952)
@@ -0,0 +1,28 @@
+Copyright (c) 2007-2011 Xiph.org Foundation
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions
+are met:
+
+- Redistributions of source code must retain the above copyright
+notice, this list of conditions and the following disclaimer.
+
+- Redistributions in binary form must reproduce the above copyright
+notice, this list of conditions and the following disclaimer in the
+documentation and/or other materials provided with the distribution.
+
+- Neither the name of the Xiph.org Foundation nor the names of its
+contributors may be used to endorse or promote products derived from
+this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION
+OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

Copied: trunk/chirptest/INSTALL (from rev 13251, trunk/ghost/INSTALL)
===================================================================
--- trunk/chirptest/INSTALL	                        (rev 0)
+++ trunk/chirptest/INSTALL	2011-05-03 20:34:45 UTC (rev 17952)
@@ -0,0 +1,18 @@
+This code builds a single, standalone graphing application 'chirp',
+which puts the linearized sinusoidal/chirp estimation code through its
+paces with varied inputs and graphs the results.
+
+REQUIREMENTS:
+
+A GNU-ish build system (automake, autoconf, GNU make)
+Cairo development libraries (libcairo and headers)
+A reasonably serious amount of CPU
+
+BUILD (from SVN):
+
+./autogen.sh
+make
+
+RUN:
+
+./chirp

Copied: trunk/chirptest/Makefile.am (from rev 13251, trunk/ghost/Makefile.am)
===================================================================
--- trunk/chirptest/Makefile.am	                        (rev 0)
+++ trunk/chirptest/Makefile.am	2011-05-03 20:34:45 UTC (rev 17952)
@@ -0,0 +1,17 @@
+## Process this file with automake to produce Makefile.in. -*-Makefile-*-
+
+# To disable automatic dependency tracking if using other tools than
+# gcc and gmake, add the option 'no-dependencies'
+AUTOMAKE_OPTIONS = 1.6
+
+EXTRA_DIST = COPYING autogen.sh
+
+bin_PROGRAMS = chirp
+
+chirp_SOURCES = chirp.c chirp.h chirpgraph.c chirpgraph.h window.c window.h chirptest.c 
+
+debug:
+        $(MAKE) all CFLAGS="@DEBUG@"
+
+profile:
+        $(MAKE) all CFLAGS="@PROFILE@"

Copied: trunk/chirptest/README (from rev 13251, trunk/ghost/README)
===================================================================
--- trunk/chirptest/README	                        (rev 0)
+++ trunk/chirptest/README	2011-05-03 20:34:45 UTC (rev 17952)
@@ -0,0 +1,5 @@
+This is basic benchmarking code for the algorithm presented in "An
+Iterative Linearised Solution to the Sinusoidal Parameter Estimation
+Problem".  The intent is to provide known good results for comparison
+with independent implementations and future improvements.
+

Copied: trunk/chirptest/autogen.sh (from rev 13251, trunk/ghost/autogen.sh)
===================================================================
--- trunk/chirptest/autogen.sh	                        (rev 0)
+++ trunk/chirptest/autogen.sh	2011-05-03 20:34:45 UTC (rev 17952)
@@ -0,0 +1,128 @@
+#!/bin/sh
+# Run this to set up the build system: configure, makefiles, etc.
+# (based on the version in enlightenment's cvs)
+
+package="chirptest"
+
+ACLOCAL_FLAGS=""
+
+olddir=`pwd`
+srcdir=`dirname $0`
+test -z "$srcdir" && srcdir=.
+
+cd "$srcdir"
+DIE=0
+
+/bin/echo "checking for autoconf... "
+(autoconf --version) < /dev/null > /dev/null 2>&1 || {
+        echo
+        echo "You must have autoconf installed to compile $package."
+        echo "Download the appropriate package for your distribution,"
+        echo "or get the source tarball at ftp://ftp.gnu.org/pub/gnu/"
+        DIE=1
+}
+
+VERSIONGREP="sed -e s/.*[^0-9\.]\([0-9][0-9]*\.[0-9][0-9]*\).*/\1/"
+VERSIONMKMAJ="sed -e s/\([0-9][0-9]*\)[^0-9].*/\\1/"
+VERSIONMKMIN="sed -e s/.*[0-9][0-9]*\.//"
+
+# do we need automake?
+if test -r Makefile.am; then
+  AM_OPTIONS=`fgrep AUTOMAKE_OPTIONS Makefile.am`
+  AM_NEEDED=`echo $AM_OPTIONS | $VERSIONGREP`
+  if test x"$AM_NEEDED" = "x$AM_OPTIONS"; then
+    AM_NEEDED=""
+  fi
+  if test -z $AM_NEEDED; then
+    /bin/echo -n "checking for automake... "
+    AUTOMAKE=automake
+    ACLOCAL=aclocal
+    if ($AUTOMAKE --version < /dev/null > /dev/null 2>&1); then
+      /bin/echo "yes"
+    else
+      /bin/echo "no"
+      AUTOMAKE=
+    fi
+  else
+    /bin/echo -n "checking for automake $AM_NEEDED or later... "
+    majneeded=`echo $AM_NEEDED | $VERSIONMKMAJ`
+    minneeded=`echo $AM_NEEDED | $VERSIONMKMIN`
+    for am in automake-$AM_NEEDED automake$AM_NEEDED \
+	automake automake-1.7 automake-1.8 automake-1.9 \
+        automake-1.10 automake-1.11; do
+      ($am --version < /dev/null > /dev/null 2>&1) || continue
+      ver=`$am --version < /dev/null | head -n 1 | $VERSIONGREP`
+      maj=`echo $ver | $VERSIONMKMAJ`
+      min=`echo $ver | $VERSIONMKMIN`
+      if test $maj -eq $majneeded -a $min -ge $minneeded; then
+        AUTOMAKE=$am
+        /bin/echo $AUTOMAKE
+        break
+      fi
+    done
+    test -z $AUTOMAKE &&  /bin/echo "no"
+    /bin/echo -n "checking for aclocal $AM_NEEDED or later... "
+    for ac in aclocal-$AM_NEEDED aclocal$AM_NEEDED \
+	aclocal aclocal-1.7 aclocal-1.8 aclocal-1.9 aclocal-1.10 aclocal-1.11; do
+      ($ac --version < /dev/null > /dev/null 2>&1) || continue
+      ver=`$ac --version < /dev/null | head -n 1 | $VERSIONGREP`
+      maj=`echo $ver | $VERSIONMKMAJ`
+      min=`echo $ver | $VERSIONMKMIN`
+      if test $maj -eq $majneeded -a $min -ge $minneeded; then
+        ACLOCAL=$ac
+        /bin/echo $ACLOCAL
+        break
+      fi
+    done
+    test -z $ACLOCAL && /bin/echo "no"
+  fi
+  test -z $AUTOMAKE || test -z $ACLOCAL && {
+        echo
+        echo "You must have automake installed to compile $package."
+        echo "Download the appropriate package for your distribution,"
+        echo "or get the source tarball at ftp://ftp.gnu.org/pub/gnu/"
+        exit 1
+  }
+fi
+
+/bin/echo -n "checking for libtool... "
+for LIBTOOLIZE in libtoolize glibtoolize nope; do
+  ($LIBTOOLIZE --version) < /dev/null > /dev/null 2>&1 && break
+done
+if test x$LIBTOOLIZE = xnope; then
+  /bin/echo "nope."
+  LIBTOOLIZE=libtoolize
+else
+  /bin/echo $LIBTOOLIZE
+fi
+($LIBTOOLIZE --version) < /dev/null > /dev/null 2>&1 || {
+	echo
+	echo "You must have libtool installed to compile $package."
+	echo "Download the appropriate package for your system,"
+	echo "or get the source from one of the GNU ftp sites"
+	echo "listed in http://www.gnu.org/order/ftp.html"
+	DIE=1
+}
+
+if test "$DIE" -eq 1; then
+        exit 1
+fi
+
+if test -z "$*"; then
+        echo "I am going to run ./configure with no arguments - if you wish "
+        echo "to pass any to it, please specify them on the $0 command line."
+fi
+
+/bin/echo "Generating configuration files for $package, please wait...."
+
+/bin/echo "  $ACLOCAL $ACLOCAL_FLAGS"
+$ACLOCAL $ACLOCAL_FLAGS || exit 1
+/bin/echo "  $LIBTOOLIZE --automake --force"
+$LIBTOOLIZE --automake --force || exit 1
+/bin/echo "  $AUTOMAKE --add-missing $AUTOMAKE_FLAGS"
+$AUTOMAKE --add-missing $AUTOMAKE_FLAGS || exit 1
+/bin/echo "  autoconf"
+autoconf || exit 1
+
+cd $olddir
+$srcdir/configure "$@" && /bin/echo

Copied: trunk/chirptest/chirp.c (from rev 17948, trunk/ghost/monty/chirp/chirp.c)
===================================================================
--- trunk/chirptest/chirp.c	                        (rev 0)
+++ trunk/chirptest/chirp.c	2011-05-03 20:34:45 UTC (rev 17952)
@@ -0,0 +1,1019 @@
+/********************************************************************
+ *                                                                  *
+ * THIS FILE IS PART OF THE OggGhost SOFTWARE CODEC SOURCE CODE.    *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
+ * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
+ * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
+ *                                                                  *
+ * THE OggGhost SOURCE CODE IS (C) COPYRIGHT 2007-2011              *
+ * by the Xiph.Org Foundation http://www.xiph.org/                  *
+ *                                                                  *
+ ********************************************************************
+
+ function: research-grade chirp estimation code
+ last mod: $Id$
+
+ Provides a set of chirp estimation algorithm variants for comparison
+ and characterization.  This is canned code meant for testing and
+ exploration of sinusoidal estimation.  It's designed to be held still
+ and used for reference and comparison against later improvement.
+ Please, don't optimize this code.  Test optimizations elsewhere in
+ code to be compared to this code.  Roll new technique improvements
+ into this reference set only after they have been completed.
+
+ ********************************************************************/
+
+#define _GNU_SOURCE
+#include <math.h>
+#include "chirp.h"
+#include "scales.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+static int cmp_descending_A(const void *p1, const void *p2){
+  chirp *A = (chirp *)p1;
+  chirp *B = (chirp *)p2;
+
+  return (A->A<B->A) - (B->A<A->A);
+}
+static int cmp_ascending_W(const void *p1, const void *p2){
+  chirp *A = (chirp *)p1;
+  chirp *B = (chirp *)p2;
+
+  return (B->W<A->W) - (A->W<B->W);
+}
+
+/* The stimators project a residue vecotr onto basis functions; the
+   results of the accumulation map directly to the parameters we're
+   interested in (A, P, W, dA, dW, ddA) and vice versa.  The mapping
+   equations are broken out below. */
+
+static float toA(float Ai, float Bi){
+  return hypotf(Ai,Bi);
+}
+static float toP(float Ai, float Bi){
+  return -atan2f(Bi, Ai);
+}
+static float toW(float Ai, float Bi, float Ci, float Di){
+  return (Ai*Ai || Bi*Bi) ? (Ci*Bi - Di*Ai)/(Ai*Ai + Bi*Bi) : 0;
+}
+static float todA(float Ai, float Bi, float Ci, float Di){
+  return (Ai*Ai || Bi*Bi) ? (Ci*Ai + Di*Bi)/hypotf(Ai,Bi) : 0;
+}
+static float todW(float Ai, float Bi, float Ei, float Fi){
+  return (Ai*Ai || Bi*Bi) ? (Ei*Bi - Fi*Ai)/(Ai*Ai + Bi*Bi) : 0;
+}
+static float toddA(float Ai, float Bi, float Ei, float Fi){
+  return (Ai*Ai || Bi*Bi) ? (Ei*Ai + Fi*Bi)/hypotf(Ai,Bi) : 0;
+}
+
+static float toAi(float A, float P){
+  return A * cosf(P);
+}
+static float toBi(float A, float P){
+  return A * -sinf(P);
+}
+
+static float WtoCi(float A, float P, float W){
+  return -W*A*sinf(P);
+}
+
+static float WtoDi(float A, float P, float W){
+  return -W*A*cosf(P);
+}
+
+static float dWtoEi(float A, float P, float dW){
+  return -A*dW*sin(P);
+}
+static float dWtoFi(float A, float P, float dW){
+  return -A*dW*cos(P);
+}
+
+static float dAtoCi(float P, float dA){
+  return dA*cosf(P);
+}
+static float dAtoDi(float P, float dA){
+  return -dA*sinf(P);
+}
+static float ddAtoEi(float P, float ddA){
+  return ddA*cos(P);
+}
+static float ddAtoFi(float P, float ddA){
+  return -ddA*sin(P);
+}
+
+/* fully nonlinear estimation iterator; it restarts the basis
+   functions (W and dW) and error vector for each chirp after each new
+   fit. Reconstruction (and thus error) are exact calculations. */
+static int full_nonlinear_iterate(const float *x,
+                                  const float *window,
+                                  int len,
+                                  chirp *cc,
+                                  int n,
+
+                                  float fit_limit,
+                                  int iter_limit,
+                                  int fit_gs,
+
+                                  int fitW,
+                                  int fitdA,
+                                  int fitdW,
+                                  int fitddA,
+                                  float fit_W_alpha,
+                                  float fit_dW_alpha,
+
+                                  int symm_norm,
+                                  float E0,
+                                  float E1,
+                                  float E2,
+
+                                  int bound_edges){
+  int i,j;
+  int flag=1;
+  float r[len];
+  float y[len];
+
+  /* Build a starting reconstruction based on the passed-in initial
+     chirp estimates */
+  memset(y,0,sizeof(*y)*len);
+  for(i=0;i<n;i++){
+    for(j=0;j<len;j++){
+      double jj = j-len*.5+.5;
+      float a = cc[i].A + (cc[i].dA + cc[i].ddA*jj)*jj;
+      y[j] += a*cos((cc[i].W + cc[i].dW*jj)*jj + cc[i].P);
+    }
+  }
+
+  /* outer fit iteration */
+  while(flag && iter_limit>0){
+    flag=0;
+
+    /* precompute the portion of the projection/fit estimate shared by
+       the zero, first and second order fits.  Subtracts the current
+       best fits from the input signal and windows the result. */
+    for(j=0;j<len;j++)
+      r[j]=(x[j]-y[j])*window[j];
+    memset(y,0,sizeof(*y)*len);
+
+    /* Sort chirps by descending amplitude */
+    qsort(cc, n, sizeof(*cc), cmp_descending_A);
+
+    for(i=0;i<n;i++){
+      chirp *c = cc+i;
+      float aP=0, bP=0;
+      float cP=0, dP=0;
+      float eP=0, fP=0;
+      float cP2=0, dP2=0;
+      float eP2=0, fP2=0;
+      float eP3=0, fP3=0;
+      float aE=0, bE=0;
+      float cE=0, dE=0;
+      float eE=0, fE=0;
+      float aC = cos(c->P);
+      float aS = sin(c->P);
+      float move;
+
+      for(j=0;j<len;j++){
+
+        /* no part of the nonlinear algorithm requires double
+           precision, however the largest source of noise will be from
+           a naive computation of the instantaneous basis or
+           reconstruction chirp phase.  Because dW is usually very
+           small compared to W (especially c->W*jj), (W + dW*jj)*jj
+           quickly becomes noticably short of significant bits.
+
+           We can either calculate everything mod 2PI, use a double
+           for just this calculation (and passing the result to
+           sincos) or decide we don't care.  In the production code we
+           will probably not care as most of the depth in the
+           estimator isn't needed.  Here we'll take the easiest route
+           to have a deep reference and just use doubles. */
+
+        double jj = j-len*.5+.5;
+        float jj2 = jj*jj;
+        double co,si;
+        float c2,s2;
+        float yy=r[j];
+
+        sincos((c->W + c->dW*jj)*jj,&si,&co);
+
+        si*=window[j];
+        co*=window[j];
+
+        /* add the current estimate back to the residue vector */
+        r[j] += (aC*co-aS*si) * (c->A + (c->dA + c->ddA*jj)*jj);
+
+        c2 = co*co*jj;
+        s2 = si*si*jj;
+
+        /* zero order projection */
+        aP += co*yy;
+        bP += si*yy;
+        /* first order projection */
+        cP += co*yy*jj;
+        dP += si*yy*jj;
+        /* second order projection */
+        eP += co*yy*jj2;
+        fP += si*yy*jj2;
+
+        if(fit_gs){
+          /* subtract zero order estimate from first */
+          cP2 += c2;
+          dP2 += s2;
+          /* subtract zero order estimate from second */
+          eP2 += c2*jj;
+          fP2 += s2*jj;
+          /* subtract first order estimate from second */
+          eP3 += c2*jj2;
+          fP3 += s2*jj2;
+        }
+
+        if(!symm_norm){
+          /* accumulate per-basis energy for normalization */
+          aE += co*co;
+          bE += si*si;
+          cE += c2*jj;
+          dE += s2*jj;
+          eE += c2*jj*jj2;
+          fE += s2*jj*jj2;
+        }
+      }
+
+      /* normalize, complete projections */
+      if(symm_norm){
+        aP *= E0;
+        bP *= E0;
+        cP = (cP-aP*cP2)*E1;
+        dP = (dP-bP*dP2)*E1;
+        eP = (eP-aP*eP2-cP*eP3)*E2;
+        fP = (fP-bP*fP2-dP*fP3)*E2;
+      }else{
+        aP = (aE>1e-20f?aP/aE:0.f);
+        bP = (bE>1e-20f?bP/bE:0.f);
+        cP = (cE>1e-20f?(cP-aP*cP2)/cE:0.f);
+        dP = (dE>1e-20f?(dP-bP*dP2)/dE:0.f);
+        eP = (eE>1e-20f?(eP-aP*eP2-cP*eP3)/eE:0.f);
+        fP = (fE>1e-20f?(fP-bP*fP2-dP*fP3)/fE:0.f);
+      }
+
+      if(!fitW && !fitdA)
+        cP=dP=0;
+      if(!fitdW && !fitddA)
+        eP=fP=0;
+
+      move = aP*aP + bP*bP;
+
+      /* we're fitting to the remaining error; add the fit to date
+         back in to relate our newest incremental results to the
+         global fit so far.  Note that this does not include W (or dW
+         if it is also recentered), as they're already 'subtracted'
+         when the bases are recentered each iteration */
+      aP += toAi(c->A, c->P);
+      bP += toBi(c->A, c->P);
+      cP += dAtoCi(c->P,c->dA);
+      dP += dAtoDi(c->P,c->dA);
+      eP += ddAtoEi(c->P,c->ddA);
+      fP += ddAtoFi(c->P,c->ddA);
+
+      /* guard overflow; if we're this far out, assume we're never
+         coming back. drop out now. */
+      if((aP*aP + bP*bP)>1e10 ||
+         (cP*cP + dP*dP)>1e10 ||
+         (eP*eP + fP*fP)>1e10){
+        iter_limit=0;
+        i=n;
+        break;
+      }
+
+      {
+        chirp old = *c;
+        float cc=0,dd=0;
+        float ee=0,ff=0;
+
+        /* save new estimate */
+        c->A = toA(aP,bP);
+        c->P = toP(aP,bP);
+        c->W += fit_W_alpha*(fitW ? toW(aP,bP,cP,dP) : 0);
+        c->dA = (fitdA ? todA(aP,bP,cP,dP) : 0);
+        c->dW += fit_dW_alpha*(fitdW ? todW(aP,bP,eP,fP) : 0);
+        c->ddA = (fitddA ? toddA(aP,bP,eP,fP) : 0);
+
+        /* base convergence on basis projection movement this
+           iteration, but only consider the contribution of terms we fit. */
+
+        if(fitW){
+          float W = c->W - old.W;
+          cc += WtoCi(c->A,c->P,W);
+          dd += WtoDi(c->A,c->P,W);
+        }
+        if(fitdA){
+          float dA = c->dA - old.dA;
+          cc += dAtoCi(c->P,dA);
+          dd += dAtoDi(c->P,dA);
+        }
+        if(fitdW){
+          float dW = c->dW - old.dW;
+          ee += dWtoEi(c->A,c->P,dW);
+          ff += dWtoFi(c->A,c->P,dW);
+        }
+        if(fitddA){
+          float ddA = c->ddA - old.ddA;
+          ee += ddAtoEi(c->P,ddA);
+          ff += ddAtoFi(c->P,ddA);
+        }
+        move = move + cc*cc + dd*dd + ee*ee + ff*ff;
+        if(fit_limit>0 && move>fit_limit*fit_limit)flag=1;
+        if(fit_limit<0 && move>1e-14)flag=1;
+      }
+
+      if(bound_edges){
+        /* do not allow negative or trans-Nyquist center frequencies */
+        if(c->W<0){
+          c->W = 0; /* clamp frequency to 0 (DC) */
+          c->dW = 0; /* if W is 0, the chirp rate must also be 0 to
+                                 avoid negative frequencies */
+        }
+        if(c->W>M_PI){
+          c->W = M_PI; /* clamp frequency to Nyquist */
+          c->dW = 0; /* if W is at Nyquist, the chirp rate must
+                           also be 0 to avoid trans-Nyquist
+                           frequencies */
+        }
+        /* Just like with the center frequency, don't allow the
+           chirp rate (dW) parameter to cross over to a negative
+           instantaneous frequency */
+        if(c->W+c->dW*len < 0)
+          c->dW = -c->W/len;
+        if(c->W-c->dW*len < 0)
+          c->dW = c->W/len;
+        /* ...or exceed Nyquist */
+        if(c->W + c->dW*len > M_PI)
+          c->dW = M_PI/len - c->W/len;
+        if(c->W - c->dW*len > M_PI)
+          c->dW = c->W/len - M_PI/len;
+      }
+
+      /* update the reconstruction/residue vectors with new fit */
+      {
+        for(j=0;j<len;j++){
+          double jj = j-len*.5+.5;
+          float a = c->A + (c->dA + c->ddA*jj)*jj;
+          float v = a*cos(c->P + (c->W + c->dW*jj)*jj);
+          r[j] -= v*window[j];
+          y[j] += v;
+        }
+      }
+    }
+    if(flag)iter_limit--;
+  }
+  return iter_limit;
+}
+
+
+/* partial-nonlinear estimation iterator; it restarts the basis
+   functions (W only).  Reconstruction and error are approximated from
+   basis (as done in the original paper). Basis function is not
+   chirped, regardless of initial dW estimate */
+static int partial_nonlinear_iterate(const float *x,
+                                     const float *window,
+                                     int len,
+                                     chirp *ch,
+                                     int n,
+
+                                     float fit_limit,
+                                     int   iter_limit,
+                                     int   fit_gs,
+
+                                     int fitW,
+                                     int fitdA,
+                                     int fitdW,
+                                     int fitddA,
+                                     float fit_W_alpha,
+
+                                     int symm_norm,
+                                     float E0,
+                                     float E1,
+                                     float E2){
+  int i,j;
+  float y[len];
+  int flag=1;
+
+  float *cos_table[n];
+  float *sin_table[n];
+  float *tcos_table[n];
+  float *tsin_table[n];
+  float *ttcos_table[n];
+  float *ttsin_table[n];
+  float cosE[n], sinE[n];
+  float tcosE[n], tsinE[n];
+  float ttcosE[n], ttsinE[n];
+  float tcosC2[n], tsinC2[n];
+  float ttcosC2[n], ttsinC2[n];
+  float ttcosC3[n], ttsinC3[n];
+
+  float ai[n], bi[n];
+  float ci[n], di[n];
+  float ei[n], fi[n];
+
+  for (i=0;i<n;i++){
+    cos_table[i]=calloc(len,sizeof(**cos_table));
+    sin_table[i]=calloc(len,sizeof(**sin_table));
+    tcos_table[i]=calloc(len,sizeof(**tcos_table));
+    tsin_table[i]=calloc(len,sizeof(**tsin_table));
+    ttcos_table[i]=calloc(len,sizeof(**ttcos_table));
+    ttsin_table[i]=calloc(len,sizeof(**ttsin_table));
+  }
+
+  /* outer fit iteration */
+  while(flag && iter_limit>0){
+    flag=0;
+    memset(y,0,sizeof(y));
+
+    /* Sort chirps by descending amplitude */
+    qsort(ch, n, sizeof(*ch), cmp_descending_A);
+
+    /* Calculate new basis functions */
+    for (i=0;i<n;i++){
+      float tmpa=0;
+      float tmpb=0;
+      float tmpc=0;
+      float tmpd=0;
+      float tmpe=0;
+      float tmpf=0;
+
+      float tmpc2=0;
+      float tmpd2=0;
+      float tmpe3=0;
+      float tmpf3=0;
+
+      /* chirp param -> basis acc */
+      ai[i] = toAi(ch[i].A, ch[i].P);
+      bi[i] = toBi(ch[i].A, ch[i].P);
+      ci[i] = dAtoCi(ch[i].P, ch[i].dA);
+      di[i] = dAtoDi(ch[i].P, ch[i].dA);
+      ei[i] = ddAtoEi(ch[i].P, ch[i].ddA) + dWtoEi(ch[i].A, ch[i].P, ch[i].dW);
+      fi[i] = ddAtoFi(ch[i].P, ch[i].ddA) + dWtoFi(ch[i].A, ch[i].P, ch[i].dW);
+
+      for (j=0;j<len;j++){
+        double jj = j-len/2.+.5;
+        double c,s;
+        sincos(ch[i].W*jj,&s,&c);
+
+        /* save basis funcs */
+        /* calc reconstruction vector */
+        y[j] +=
+          ai[i]*(cos_table[i][j] = c) +
+          bi[i]*(sin_table[i][j] = s) +
+          ci[i]*(tcos_table[i][j] = jj*c) +
+          di[i]*(tsin_table[i][j] = jj*s) +
+          ei[i]*(ttcos_table[i][j] = jj*jj*c) +
+          fi[i]*(ttsin_table[i][j] = jj*jj*s);
+
+        /* The modulation terms */
+        tmpc += tcos_table[i][j]*tcos_table[i][j]*window[j]*window[j];
+        tmpd += tsin_table[i][j]*tsin_table[i][j]*window[j]*window[j];
+
+        if(!symm_norm){
+          /* The sinusoidal terms */
+          tmpa += cos_table[i][j]*cos_table[i][j]*window[j]*window[j];
+          tmpb += sin_table[i][j]*sin_table[i][j]*window[j]*window[j];
+
+          /* The second order modulations */
+          tmpe += ttcos_table[i][j]*ttcos_table[i][j]*window[j]*window[j];
+          tmpf += ttsin_table[i][j]*ttsin_table[i][j]*window[j]*window[j];
+        }
+
+        /* gs fit terms */
+        if(fit_gs){
+          tmpc2 += cos_table[i][j]*tcos_table[i][j]*window[j]*window[j];
+          tmpd2 += sin_table[i][j]*tsin_table[i][j]*window[j]*window[j];
+          tmpe3 += tcos_table[i][j]*ttcos_table[i][j]*window[j]*window[j];
+          tmpf3 += tsin_table[i][j]*ttsin_table[i][j]*window[j]*window[j];
+        }
+      }
+
+      /* Set basis normalizations */
+      if(symm_norm){
+        cosE[i] = sinE[i] = E0;
+        tcosE[i] = tsinE[i] = E1;
+        ttcosE[i] = ttsinE[i] = E2;
+      }else{
+        cosE[i] = (tmpa>1e-20f ? 1.f/tmpa : 0);
+        sinE[i] = (tmpb>1e-20f ? 1.f/tmpb : 0);
+        tcosE[i] = (tmpc>1e-20f ? 1.f/tmpc : 0);
+        tsinE[i] = (tmpd>1e-20f ? 1.f/tmpd : 0);
+        ttcosE[i] = (tmpe>1e-20f ? 1.f/tmpe : 0);
+        ttsinE[i] = (tmpf>1e-20f ? 1.f/tmpf : 0);
+      }
+
+      if(fit_gs){
+        /* set gs fit terms */
+        tcosC2[i] = tmpc2;
+        tsinC2[i] = tmpd2;
+        ttcosC2[i] = tmpc;
+        ttsinC2[i] = tmpd;
+        ttcosC3[i] = tmpe3;
+        ttsinC3[i] = tmpf3;
+      }
+    }
+
+    for (i=0;i<n;i++){
+
+      float tmpa=0, tmpb=0;
+      float tmpc=0, tmpd=0;
+      float tmpe=0, tmpf=0;
+
+      /* (Sort of) project the residual on the four basis functions */
+      for (j=0;j<len;j++){
+        float r = (x[j]-y[j])*window[j]*window[j];
+        tmpa += r*cos_table[i][j];
+        tmpb += r*sin_table[i][j];
+        tmpc += r*tcos_table[i][j];
+        tmpd += r*tsin_table[i][j];
+        tmpe += r*ttcos_table[i][j];
+        tmpf += r*ttsin_table[i][j];
+      }
+
+      tmpa*=cosE[i];
+      tmpb*=sinE[i];
+
+      if(fit_gs){
+        tmpc -= tmpa*tcosC2[i];
+        tmpd -= tmpb*tsinC2[i];
+      }
+
+      tmpc*=tcosE[i];
+      tmpd*=tsinE[i];
+
+      if(fit_gs){
+        tmpe -= tmpa*ttcosC2[i] + tmpc*ttcosC3[i];
+        tmpf -= tmpb*ttsinC2[i] + tmpd*ttsinC3[i];
+      }
+
+      tmpe*=ttcosE[i];
+      tmpf*=ttsinE[i];
+
+      if(!fitW && !fitdA)
+        tmpc=tmpd=0;
+      if(!fitdW && !fitddA)
+        tmpe=tmpf=0;
+
+      /* Update the signal approximation for the next iteration */
+      for (j=0;j<len;j++){
+        y[j] +=
+          tmpa*cos_table[i][j]+
+          tmpb*sin_table[i][j]+
+          tmpc*tcos_table[i][j]+
+          tmpd*tsin_table[i][j]+
+          tmpe*ttcos_table[i][j]+
+          tmpf*ttsin_table[i][j];
+      }
+
+      ai[i] += tmpa;
+      bi[i] += tmpb;
+      ci[i] += tmpc;
+      di[i] += tmpd;
+      ei[i] += tmpe;
+      fi[i] += tmpf;
+
+      /* guard overflow; if we're this far out, assume we're never
+         coming back. drop out now. */
+      if((ai[i]*ai[i] + bi[i]*bi[i])>1e10 ||
+         (ci[i]*ci[i] + di[i]*di[i])>1e10 ||
+         (ei[i]*ei[i] + fi[i]*fi[i])>1e10){
+        iter_limit=0;
+        i=n;
+      }
+
+      /* save new estimate */
+      ch[i].A = toA(ai[i],bi[i]);
+      ch[i].P = toP(ai[i],bi[i]);
+      ch[i].W += fit_W_alpha*(fitW ? toW(ai[i],bi[i],ci[i],di[i]) : 0);
+      ch[i].dA = (fitdA ? todA(ai[i],bi[i],ci[i],di[i]) : 0);
+      ch[i].dW = (fitdW ? todW(ai[i],bi[i],ei[i],fi[i]) : 0);
+      ch[i].ddA = (fitddA ? toddA(ai[i],bi[i],ei[i],fi[i]) : 0);
+
+      /* base convergence on basis projection movement this
+         iteration, but only consider the contribution of terms we fit. */
+      {
+        float move;
+        float cc=0,dd=0;
+        float ee=0,ff=0;
+        if(fitW){
+          float W = toW(ai[i],bi[i],ci[i],di[i]);
+          cc += WtoCi(ch[i].A,ch[i].P,W);
+          dd += WtoDi(ch[i].A,ch[i].P,W);
+        }
+        if(fitdA){
+          float dA = todA(ai[i],bi[i],tmpc,tmpd);
+          cc += dAtoCi(ch[i].P,dA);
+          dd += dAtoDi(ch[i].P,dA);
+        }
+        if(fitdW){
+          float dW = todW(ai[i],bi[i],tmpe,tmpf);
+          ee += dWtoEi(ch[i].A,ch[i].P,dW);
+          ff += dWtoFi(ch[i].A,ch[i].P,dW);
+        }
+        if(fitddA){
+          float ddA = toddA(ai[i],bi[i],tmpe,tmpf);
+          ee += ddAtoEi(ch[i].P,ddA);
+          ff += ddAtoFi(ch[i].P,ddA);
+        }
+        move = tmpa*tmpa + tmpb*tmpb + cc*cc + dd*dd + ee*ee + ff*ff;
+        if(fit_limit>0 && move>fit_limit*fit_limit)flag=1;
+        if(fit_limit<0 && move>1e-14)flag=1;
+      }
+    }
+    if(flag)iter_limit--;
+  }
+
+  for(i=0;i<n;i++){
+    free(cos_table[i]);
+    free(sin_table[i]);
+    free(tcos_table[i]);
+    free(tsin_table[i]);
+    free(ttcos_table[i]);
+    free(ttsin_table[i]);
+  }
+
+  return iter_limit;
+}
+
+/* linear estimation iterator; sets fixed basis functions for each
+   chirp according to the passed in initial estimates. This code
+   follows the paper as exactly as possible (reconstruction is
+   estimated from basis, basis is not chirped regardless of initial
+   estimate, etc). */
+static int linear_iterate(const float *x,
+                          const float *window,
+                          int len,
+                          chirp *ch,
+                          int n,
+
+                          float fit_limit,
+                          int iter_limit,
+                          int fit_gs,
+
+                          int fitW,
+                          int fitdA,
+                          int fitdW,
+                          int fitddA,
+                          int symm_norm,
+                          float E0,
+                          float E1,
+                          float E2){
+
+  float *cos_table[n];
+  float *sin_table[n];
+  float *tcos_table[n];
+  float *tsin_table[n];
+  float *ttcos_table[n];
+  float *ttsin_table[n];
+  float cosE[n], sinE[n];
+  float tcosE[n], tsinE[n];
+  float ttcosE[n], ttsinE[n];
+  float tcosC2[n], tsinC2[n];
+  float ttcosC2[n], ttsinC2[n];
+  float ttcosC3[n], ttsinC3[n];
+
+  float ai[n], bi[n];
+  float ci[n], di[n];
+  float ei[n], fi[n];
+  int i,j;
+  int flag=1;
+  float y[len];
+
+  memset(y,0,sizeof(y));
+
+  for (i=0;i<n;i++){
+    float tmpa=0;
+    float tmpb=0;
+    float tmpc=0;
+    float tmpd=0;
+    float tmpe=0;
+    float tmpf=0;
+
+    float tmpc2=0;
+    float tmpd2=0;
+    float tmpe3=0;
+    float tmpf3=0;
+
+    /* seed the basis accumulators from our passed in estimated parameters */
+    /* Don't add in W; this is already included by the basis */
+    ai[i] = toAi(ch[i].A, ch[i].P);
+    bi[i] = toBi(ch[i].A, ch[i].P);
+    ci[i] = dAtoCi(ch[i].P,ch[i].dA);
+    di[i] = dAtoDi(ch[i].P,ch[i].dA);
+    ei[i] = ddAtoEi(ch[i].P, ch[i].ddA) + dWtoEi(ch[i].A, ch[i].P, ch[i].dW);
+    fi[i] = ddAtoFi(ch[i].P, ch[i].ddA) + dWtoFi(ch[i].A, ch[i].P, ch[i].dW);
+
+    /* Build a table for the basis functions at each frequency */
+    cos_table[i]=calloc(len,sizeof(**cos_table));
+    sin_table[i]=calloc(len,sizeof(**sin_table));
+    tcos_table[i]=calloc(len,sizeof(**tcos_table));
+    tsin_table[i]=calloc(len,sizeof(**tsin_table));
+    ttcos_table[i]=calloc(len,sizeof(**ttcos_table));
+    ttsin_table[i]=calloc(len,sizeof(**ttsin_table));
+
+    for (j=0;j<len;j++){
+      double jj = j-len/2.+.5;
+      double c,s;
+      sincos(ch[i].W*jj,&s,&c);
+
+      /* save basis funcs */
+      /* initial reconstruction */
+      y[j] +=
+        ai[i]*(cos_table[i][j] = c) +
+        bi[i]*(sin_table[i][j] = s) +
+        ci[i]*(tcos_table[i][j] = jj*c) +
+        di[i]*(tsin_table[i][j] = jj*s) +
+        ei[i]*(ttcos_table[i][j] = jj*jj*c) +
+        fi[i]*(ttsin_table[i][j] = jj*jj*s);
+
+      /* The sinusoidal terms */
+      tmpa += cos_table[i][j]*cos_table[i][j]*window[j]*window[j];
+      tmpb += sin_table[i][j]*sin_table[i][j]*window[j]*window[j];
+
+      /* The modulation terms */
+      tmpc += tcos_table[i][j]*tcos_table[i][j]*window[j]*window[j];
+      tmpd += tsin_table[i][j]*tsin_table[i][j]*window[j]*window[j];
+
+      /* The second order modulations */
+      tmpe += ttcos_table[i][j]*ttcos_table[i][j]*window[j]*window[j];
+      tmpf += ttsin_table[i][j]*ttsin_table[i][j]*window[j]*window[j];
+
+      /* gs fit terms */
+      tmpc2 += cos_table[i][j]*tcos_table[i][j]*window[j]*window[j];
+      tmpd2 += sin_table[i][j]*tsin_table[i][j]*window[j]*window[j];
+      tmpe3 += tcos_table[i][j]*ttcos_table[i][j]*window[j]*window[j];
+      tmpf3 += tsin_table[i][j]*ttsin_table[i][j]*window[j]*window[j];
+
+    }
+
+    /* Set basis normalizations */
+    if(symm_norm){
+      cosE[i] = sinE[i] = E0;
+      tcosE[i] = tsinE[i] = E1;
+      ttcosE[i] = ttsinE[i] = E2;
+    }else{
+      cosE[i] = (tmpa>1e-20f ? 1.f/tmpa : 0);
+      sinE[i] = (tmpb>1e-20f ? 1.f/tmpb : 0);
+      tcosE[i] = (tmpc>1e-20f ? 1.f/tmpc : 0);
+      tsinE[i] = (tmpd>1e-20f ? 1.f/tmpd : 0);
+      ttcosE[i] = (tmpe>1e-20f ? 1.f/tmpe : 0);
+      ttsinE[i] = (tmpf>1e-20f ? 1.f/tmpf : 0);
+    }
+
+    /* set gs fit terms */
+    tcosC2[i] = tmpc2;
+    tsinC2[i] = tmpd2;
+    ttcosC2[i] = tmpc;
+    ttsinC2[i] = tmpd;
+    ttcosC3[i] = tmpe3;
+    ttsinC3[i] = tmpf3;
+
+  }
+
+  while(flag && iter_limit>0){
+    flag=0;
+
+    for (i=0;i<n;i++){
+
+      float tmpa=0, tmpb=0;
+      float tmpc=0, tmpd=0;
+      float tmpe=0, tmpf=0;
+
+      /* (Sort of) project the residual on the four basis functions */
+      for (j=0;j<len;j++){
+        float r = (x[j]-y[j])*window[j]*window[j];
+        tmpa += r*cos_table[i][j];
+        tmpb += r*sin_table[i][j];
+        tmpc += r*tcos_table[i][j];
+        tmpd += r*tsin_table[i][j];
+        tmpe += r*ttcos_table[i][j];
+        tmpf += r*ttsin_table[i][j];
+      }
+
+      tmpa*=cosE[i];
+      tmpb*=sinE[i];
+
+      if(fit_gs){
+        tmpc -= tmpa*tcosC2[i];
+        tmpd -= tmpb*tsinC2[i];
+      }
+
+      tmpc*=tcosE[i];
+      tmpd*=tsinE[i];
+
+      if(fit_gs){
+        tmpe -= tmpa*ttcosC2[i] + tmpc*ttcosC3[i];
+        tmpf -= tmpb*ttsinC2[i] + tmpd*ttsinC3[i];
+      }
+
+      tmpe*=ttcosE[i];
+      tmpf*=ttsinE[i];
+
+      if(!fitW && !fitdA)
+        tmpc=tmpd=0;
+      if(!fitdW && !fitddA)
+        tmpe=tmpf=0;
+
+      /* Update the signal approximation for the next iteration */
+      for (j=0;j<len;j++){
+        y[j] +=
+          tmpa*cos_table[i][j]+
+          tmpb*sin_table[i][j]+
+          tmpc*tcos_table[i][j]+
+          tmpd*tsin_table[i][j]+
+          tmpe*ttcos_table[i][j]+
+          tmpf*ttsin_table[i][j];
+      }
+
+
+      /* guard overflow; if we're this far out, assume we're never
+         coming back. drop out now. */
+      if((ai[i]*ai[i] + bi[i]*bi[i])>1e10 ||
+         (ci[i]*ci[i] + di[i]*di[i])>1e10 ||
+         (ei[i]*ei[i] + fi[i]*fi[i])>1e10){
+        iter_limit=0;
+        i=n;
+      }
+
+      ai[i] += tmpa;
+      bi[i] += tmpb;
+      ci[i] += tmpc;
+      di[i] += tmpd;
+      ei[i] += tmpe;
+      fi[i] += tmpf;
+
+      /* base convergence on basis projection movement this
+         iteration */
+      {
+        float move = tmpa*tmpa + tmpb*tmpb +tmpc*tmpc + tmpd*tmpd + tmpe*tmpe + tmpf*tmpf;
+
+        if(fit_limit>0 && move>fit_limit*fit_limit)flag=1;
+        if(fit_limit<0 && move>1e-14)flag=1;
+      }
+
+
+    }
+    if(flag)iter_limit--;
+  }
+
+  for(i=0;i<n;i++){
+    ch[i].A = toA(ai[i],bi[i]);
+    ch[i].P = toP(ai[i],bi[i]);
+    ch[i].W += (fitW ? toW(ai[i],bi[i],ci[i],di[i]) : 0);
+    ch[i].dA = (fitdA ? todA(ai[i],bi[i],ci[i],di[i]) : 0);
+    ch[i].dW += (fitdW ? todW(ai[i],bi[i],ei[i],fi[i]) : 0);
+    ch[i].ddA = (fitddA ? toddA(ai[i],bi[i],ei[i],fi[i]) : 0);
+
+    free(cos_table[i]);
+    free(sin_table[i]);
+    free(tcos_table[i]);
+    free(tsin_table[i]);
+    free(ttcos_table[i]);
+    free(ttsin_table[i]);
+  }
+
+  return iter_limit;
+}
+
+/* Performs an iterative chirp estimation using the passed
+   in chirp paramters as an initial estimate.
+
+   x: input signal (unwindowed)
+   y: output reconstruction (unwindowed)
+   window: window to apply to input and bases during fitting process
+   len: length of x,y,r and window vectors
+   c: chirp initial estimation inputs, chirp fit outputs (sorted in frequency order)
+   n: number of chirps
+   iter_limit: maximum allowable number of iterations
+   fit_limit: desired fit quality limit
+
+   fitW : flag indicating that fitting should include the W parameter.
+   fitdA : flag indicating that fitting should include the dA parameter.
+   fitdW : flag indicating that fitting should include the dW parameter.
+   fitddA : flag indicating that fitting should include the ddA parameter.
+   symm_norm
+   int fit_gs
+   int bound_zero
+
+   Input estimates affect convergence region and speed and fit
+   quality; precise effects depend on the fit estimation algorithm
+   selected.  Note that non-zero initial values for dA, dW or ddA are
+   ignored when fitdA, fitdW, fitddA are unset.  An initial value for
+   W is always used even when W is left unfitted.
+
+   Fitting terminates when no fit parameter changes by more than
+   fit_limit in an iteration or when the fit loop reaches the
+   iteration limit.  The fit parameters as checked are scaled over the
+   length of the block.
+
+*/
+
+int estimate_chirps(const float *x,
+                    const float *window,
+                    int len,
+                    chirp *c,
+                    int n,
+
+                    float fit_limit,
+                    int iter_limit,
+                    int fit_gs,
+
+                    int fitW,
+                    int fitdA,
+                    int fitdW,
+                    int fitddA,
+                    int nonlinear,
+                    float fit_W_alpha,
+                    float fit_dW_alpha,
+                    int symm_norm,
+                    int bound_zero
+                    ){
+
+  int i;
+  float E0=0, E1=0, E2=0;
+
+  if(symm_norm){
+    for(i=0;i<len;i++){
+      float jj = i-len*.5+.5;
+      float w2 = window[i]*window[i];
+      float w2j2 = w2*jj*jj;
+      float w2j4 = w2j2*jj*jj;
+      E0 += w2;
+      E1 += w2j2;
+      E2 += w2j4;
+    }
+    E0=2/E0;
+    E1=2/E1;
+    E2=2/E2;
+  }
+
+  /* zero out initial estimates for params we'll not fit */
+  for(i=0;i<n;i++){
+    if(!fitdA) c[i].dA=0.;
+    if(!fitdW) c[i].dW=0.;
+    if(!fitddA) c[i].ddA=0.;
+  }
+
+  switch(nonlinear){
+  case 0:
+    iter_limit = linear_iterate(x,window,len,c,n,
+                                fit_limit,iter_limit,fit_gs,
+                                fitW,fitdA,fitdW,fitddA,
+                                symm_norm,E0,E1,E2);
+    break;
+  case 1:
+    iter_limit = partial_nonlinear_iterate(x,window,len,c,n,
+                                           fit_limit,iter_limit,fit_gs,
+                                           fitW,fitdA,fitdW,fitddA,
+                                           fit_W_alpha,
+                                           symm_norm,E0,E1,E2);
+    break;
+  case 2:
+    iter_limit = full_nonlinear_iterate(x,window,len,c,n,
+                                        fit_limit,iter_limit,fit_gs,
+                                        fitW,fitdA,fitdW,fitddA,
+                                        fit_W_alpha,fit_dW_alpha,
+                                        symm_norm,E0,E1,E2,
+                                        bound_zero);
+    break;
+  }
+
+  /* Sort by ascending frequency */
+  qsort(c, n, sizeof(*c), cmp_ascending_W);
+
+  return iter_limit;
+}
+
+/* advances/extrapolates the passed in chirps so that the params are
+   centered forward in time by len samples. */
+
+void advance_chirps(chirp *c, int n, int len){
+  int i;
+  for(i=0;i<n;i++){
+    c[i].A += c[i].dA*len;
+    c[i].P += (c[i].W+c[i].dW*len)*len;
+    c[i].W += c[i].dW*len*2;
+  }
+}
+
+/* dead, but leave it around; does a slow, brute force DCFT */
+void slow_dcft_row(float *x, float *y, int k, int n){
+  int i,j;
+  for(i=0;i<n/2+1;i++){ /* frequency loop */
+    float real=0.,imag=0.;
+    for(j=0;j<n;j++){ /* multiply loop */
+      float s,c;
+      float jj = (j-n/2+.5)/n;
+      sincosf(2.*M_PI*(i+k*jj)*jj,&s,&c);
+      real += c*x[j];
+      imag += s*x[j];
+    }
+    y[i] = hypot(real,imag);
+  }
+}
+

Copied: trunk/chirptest/chirp.h (from rev 17944, trunk/ghost/monty/chirp/chirp.h)
===================================================================
--- trunk/chirptest/chirp.h	                        (rev 0)
+++ trunk/chirptest/chirp.h	2011-05-03 20:34:45 UTC (rev 17952)
@@ -0,0 +1,52 @@
+/********************************************************************
+ *                                                                  *
+ * THIS FILE IS PART OF THE OggGhost SOFTWARE CODEC SOURCE CODE.    *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
+ * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
+ * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
+ *                                                                  *
+ * THE OggGhost SOURCE CODE IS (C) COPYRIGHT 2007-2011              *
+ * by the Xiph.Org Foundation http://www.xiph.org/                  *
+ *                                                                  *
+ ********************************************************************
+
+ function: research-grade chirp extraction code
+ last mod: $Id$
+
+ ********************************************************************/
+
+typedef struct {
+  float A;  /* center amplitude (linear) */
+  float W;  /* frequency (radians per sample, not cycles per block) */
+  float P;  /* phase (radians) */
+  float dA; /* amplitude modulation (linear change per sample) */
+  float dW; /* frequency modulation (radians per sample^2) */
+  float ddA; /* amplitude modulation (linear change per sample^2) */
+  int label; /* used for tracking by outside code */
+} chirp;
+
+extern int
+estimate_chirps(const float *x, /* unwindowed input to fit */
+                const float *window, /* window to apply to input/bases */
+                int len,        /* block length */
+                chirp *c,       /* list of chirp estimates/outputs */
+                int n,          /* number of chirps */
+                float fit_limit,/* minimum basis movement to continue iteration */
+                int iter_limit, /* maximum number of iterations */
+                int fit_gs,     /* Use Gauss-Seidel partial updates */
+
+                int fitW,       /* fit the W parameter */
+                int fitdA,      /* fit the dA parameter */
+                int fitdW,      /* fit the dW parameter */
+                int fitddA,     /* fit the ddA parameter */
+                int nonlinear,  /* perform a linear fit (0),
+                                   nonlinear fit recentering W only (1)
+                                   nonlinear fit recentering W and dW (2) */
+                float fit_W_alpha, /* W alpha multiplier for nonlinear fit */
+                float fit_dW_alpha,/* dW alpha multiplier for nonlinear fit */
+                int symm_norm,     /* Use symmetric normalization optimization */
+                int bound_zero);   /* prevent W or dW from fitting to negative or
+                                      greater-then-Nyquist frequencies */
+
+extern void advance_chirps(chirp *c, int n, int len);
+

Copied: trunk/chirptest/chirpgraph.c (from rev 17944, trunk/ghost/monty/chirp/chirpgraph.c)
===================================================================
--- trunk/chirptest/chirpgraph.c	                        (rev 0)
+++ trunk/chirptest/chirpgraph.c	2011-05-03 20:34:45 UTC (rev 17952)
@@ -0,0 +1,483 @@
+/********************************************************************
+ *                                                                  *
+ * THIS FILE IS PART OF THE OggGhost SOFTWARE CODEC SOURCE CODE.    *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
+ * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
+ * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
+ *                                                                  *
+ * THE OggGhost SOURCE CODE IS (C) COPYRIGHT 2007-2011              *
+ * by the Xiph.Org Foundation http://www.xiph.org/                  *
+ *                                                                  *
+ ********************************************************************
+
+ function: graphing code for chirp tests
+ last mod: $Id$
+
+ ********************************************************************/
+
+#define _GNU_SOURCE
+#include <math.h>
+#include "chirp.h"
+#include "chirpgraph.h"
+#include "scales.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <cairo/cairo.h>
+#include <pthread.h>
+
+/********************** colors for graphs *********************************/
+
+void set_iter_color(cairo_t *cC, int ret, float a){
+  if (ret>100){
+    cairo_set_source_rgba(cC,1,1,1,a); /* white */
+  }else if (ret>50){ /* 51(black)-100(white) */
+    cairo_set_source_rgba(cC,(ret-50)*.02,(ret-50)*.02,(ret-50)*.02,a);
+  }else if (ret>30){ /* 31(red) - 50(black) */
+    cairo_set_source_rgba(cC,1-(ret-30)*.05,0,0,a);
+  }else if (ret>10){ /* 11(yellow) - 30(red) */
+    cairo_set_source_rgba(cC,1,1-(ret-10)*.05,0,a);
+  }else if (ret>4){ /*  5 (green) - 10 (yellow) */
+    cairo_set_source_rgba(cC,(ret-5)*.15+.25,1,0,a);
+  }else if (ret==4){ /*  4 brighter cyan */
+    cairo_set_source_rgba(cC,0,.8,.4,a);
+  }else if (ret==3){  /* 3 cyan */
+    cairo_set_source_rgba(cC,0,.6,.6,a);
+  }else if (ret==2){  /* 2 dark cyan */
+    cairo_set_source_rgba(cC,0,.4,.8,a);
+  }else if (ret==1){ /* 1 blue */
+    cairo_set_source_rgba(cC,.1,.2,1,a);
+  }else{ /* dark blue */
+    cairo_set_source_rgba(cC,.2,0,1,a);
+  }
+}
+
+void set_iter_text_color(cairo_t *cC, int ret){
+  if (ret>75){
+    cairo_set_source_rgba(cC,0,0,0,1); /* black on white*/
+  }else if (ret>35){
+    cairo_set_source_rgba(cC,1,1,1,1); /* white on red/black  */
+  }else if (ret>3){
+    cairo_set_source_rgba(cC,0,0,0,1); /* black on cyan */
+  }else{
+    cairo_set_source_rgba(cC,1,1,1,1); /* white on dark */
+  }
+}
+
+void set_error_color(cairo_t *c, float err,float a){
+  if(isnan(err) || fabs(err)>1.){
+    cairo_set_source_rgba(c,1,1,1,a); /* white */
+  }else{
+    err=fabs(err);
+    if(err>.1){
+      cairo_set_source_rgba(c,(err-.1)/.9,(err-.1)/.9,(err-.1)/.9,a); /* white->black */
+    }else if(err>.01){
+      cairo_set_source_rgba(c,1.-((err-.01)/.09),0,0,a); /* black->red */
+    }else if(err>.001){
+      cairo_set_source_rgba(c,1,1.-((err-.001)/.009),0,a); /* red->yellow */
+    }else if(err>.0001){
+      cairo_set_source_rgba(c,((err-.0001)/.0009),1,0,a); /* yellow->green */
+    }else if(err>.00001){
+      cairo_set_source_rgba(c,0,1,1.-((err-.00001)/.00009),a); /* green->cyan */
+    }else if(err>.000001){
+      cairo_set_source_rgba(c,.2-((err-.000001)/.000009)*.2,
+                            ((err-.000001)/.000009)*.8+.2,1,a); /* cyan->blue */
+    }else{
+      cairo_set_source_rgba(c,.1,.1,1,a); /* blue */
+    }
+  }
+}
+
+void set_error_text_color(cairo_t *c, float err){
+  if(isnan(err) || fabs(err)>1.){
+    cairo_set_source_rgba(c,0,0,0,1); /* black on white */
+  }else{
+    err=fabs(err);
+    if(err>.5){
+      cairo_set_source_rgba(c,0,0,0,1); /* black on white */
+    }else if(err>.02){
+      cairo_set_source_rgba(c,1,1,1,1); /* white */
+    }else if(err>.000005){
+      cairo_set_source_rgba(c,0,0,0,1); /* black */
+    }else{
+      cairo_set_source_rgba(c,1,1,1,1); /* white */
+    }
+  }
+}
+
+/********* draw everything in the graph except the graph data itself *******/
+
+static float fontsize=12;
+static int x0s,x1s,xmajor,xmajorf;
+static int y0s,y1s,ymajor,ymajorf;
+
+/* computed configuration globals */
+int leftpad=0;
+static int rightpad=0;
+int toppad=0;
+static int bottompad=0;
+static float legendy=0;
+static float legendh=0;
+static int pic_w=0;
+static int pic_h=0;
+static float titleheight=0.;
+static int x_n=0;
+static int y_n=0;
+
+/* determines padding, etc.  Will not expand the frame to prevent
+   overruns of user-set titles/legends */
+void setup_graphs(int start_x_step,
+                  int end_x_step, /* inclusive; not one past */
+                  int x_major_d,
+                  float x_major_f,
+
+                  int start_y_step,
+                  int end_y_step, /* inclusive; not one past */
+                  int y_major_d,
+                  float y_major_f,
+
+                  int subtitles,
+                  float fs){
+
+  /* determine ~ padding needed */
+  cairo_surface_t *cs=cairo_image_surface_create(CAIRO_FORMAT_RGB24,100,100);
+  cairo_t *ct=cairo_create(cs);
+  cairo_text_extents_t extents;
+  cairo_font_extents_t fextents;
+  int y;
+
+  x_n = end_x_step-start_x_step+1;
+  y_n = end_y_step-start_y_step+1;
+  fontsize=fs;
+  x0s=start_x_step;
+  x1s=end_x_step;
+  xmajor=x_major_d;
+  xmajorf=x_major_f;
+  y0s=start_y_step;
+  y1s=end_y_step;
+  ymajor=y_major_d;
+  ymajorf=y_major_f;
+
+  /* y labels */
+  cairo_set_font_size(ct, fontsize);
+  for(y=y0s+fontsize;y<=y1s;y++){
+    if(y%ymajor==0){
+      char buf[80];
+      snprintf(buf,80,"%.0f",(float)y/ymajor*ymajorf);
+      cairo_text_extents(ct, buf, &extents);
+      if(extents.width + extents.height*3.5>leftpad)leftpad=extents.width + extents.height*3.5;
+    }
+  }
+
+  /* x labels */
+  cairo_font_extents(ct, &fextents);
+  if(fextents.height*3>bottompad)bottompad=fextents.height*3;
+
+  /* center horizontally */
+  if(leftpad<rightpad)leftpad=rightpad;
+  if(rightpad<leftpad)rightpad=leftpad;
+
+  /* top legend */
+  {
+    float sofar=0;
+    cairo_save(ct);
+    cairo_select_font_face (ct, "",
+                            CAIRO_FONT_SLANT_NORMAL,
+                            CAIRO_FONT_WEIGHT_BOLD);
+    cairo_set_font_size(ct,fontsize*1.25);
+    cairo_font_extents(ct, &fextents);
+    sofar+=fextents.height;
+    cairo_restore(ct);
+
+    cairo_font_extents(ct, &fextents);
+    while(subtitles--)
+      sofar+=fextents.height;
+    if(sofar>titleheight)titleheight=sofar;
+    if(toppad<titleheight+fontsize*2)toppad=titleheight+fontsize*2;
+  }
+
+  /* color legend */
+  {
+    float width=0.;
+    float w1,w2,w3;
+
+    cairo_text_extents(ct, "100", &extents);
+    w1=extents.width*2*12;
+    cairo_text_extents(ct, ".000001", &extents);
+    w2=extents.width*1.5*7;
+    cairo_text_extents(ct, ".0001%", &extents);
+    w3=extents.width*1.5*7;
+
+    width+=MAX(w1,MAX(w2,w3));
+
+    legendy = y_n+extents.height*4;
+    legendh = extents.height*1.5;
+
+    if(legendy+legendh*4>bottompad)
+      bottompad=(legendy-y_n)+legendh*2.5;
+  }
+
+  pic_w = x_n + leftpad + rightpad;
+  pic_h = y_n + toppad + bottompad;
+}
+
+/* draws the page surrounding the graph data itself */
+cairo_t *draw_page(char *title,
+                   char *subtitle1,
+                   char *subtitle2,
+                   char *subtitle3,
+                   char *xaxis_label,
+                   char *yaxis_label,
+                   char *legend_label,
+                   int datatype){
+
+  int i;
+  cairo_text_extents_t extents;
+  cairo_font_extents_t fextents;
+  cairo_surface_t *cs=cairo_image_surface_create(CAIRO_FORMAT_ARGB32,pic_w, pic_h);
+  if(!cs || cairo_surface_status(cs)!=CAIRO_STATUS_SUCCESS){
+    fprintf(stderr,"Could not set up Cairo surface.\n\n");
+    exit(1);
+  }
+  cairo_t *c = cairo_create(cs);
+
+  /* clear page to white */
+  cairo_set_source_rgb(c,1,1,1);
+  cairo_rectangle(c,0,0,pic_w,pic_h);
+  cairo_fill(c);
+  cairo_set_font_size(c, fontsize);
+
+  /* set graph area to transparent */
+  cairo_save(c);
+  cairo_set_source_rgba(c,0,0,0,0);
+  cairo_set_operator(c,CAIRO_OPERATOR_SOURCE);
+  cairo_rectangle(c,leftpad,toppad,x_n,y_n);
+  cairo_fill(c);
+  cairo_restore(c);
+
+  /* Y axis numeric labels */
+  cairo_set_source_rgb(c,0,0,0);
+  for(i=y0s+fontsize;i<=y1s;i++){
+    if(i%ymajor==0){
+      int y = toppad+y_n-i+y0s;
+      char buf[80];
+
+      snprintf(buf,80,"%.0f",(float)i/ymajor*ymajorf);
+      cairo_text_extents(c, buf, &extents);
+      cairo_move_to(c,leftpad - fontsize*.5 - extents.width,y+extents.height*.5);
+      cairo_show_text(c,buf);
+    }
+  }
+
+  /* X axis labels */
+  for(i=x0s;i<=x1s;i++){
+    if(i%xmajor==0){
+      char buf[80];
+      int x = leftpad + i - x0s;
+      if(i==0){
+        snprintf(buf,80,"DC");
+      }else{
+        snprintf(buf,80,"%.0f",(float)i/xmajor*xmajorf);
+      }
+      cairo_text_extents(c, buf, &extents);
+      cairo_move_to(c,x - extents.width/2,y_n+toppad+extents.height+fontsize*.5);
+      cairo_show_text(c,buf);
+    }
+  }
+
+  /* Y axis caption */
+  {
+    cairo_matrix_t a;
+    cairo_matrix_t b = {0.,-1., 1.,0., 0.,0.}; // account for border!
+    cairo_matrix_t d;
+    cairo_text_extents(c, yaxis_label, &extents);
+    cairo_move_to(c,extents.height+fontsize*.5,y_n/2+toppad+extents.width*.5);
+
+    cairo_save(c);
+    cairo_get_matrix(c,&a);
+    cairo_matrix_multiply(&d,&a,&b);
+    cairo_set_matrix(c,&d);
+    cairo_show_text(c,yaxis_label);
+    cairo_restore(c);
+  }
+
+  /* X axis caption */
+  {
+    cairo_text_extents(c, xaxis_label, &extents);
+    cairo_move_to(c,pic_w/2-extents.width/2,y_n+toppad+extents.height*2+fontsize*.25);
+    cairo_show_text(c,xaxis_label);
+  }
+
+  /* top title(s) */
+  {
+    float y = (toppad-titleheight);
+    if(title){
+      cairo_save(c);
+      cairo_select_font_face (c, "",
+                              CAIRO_FONT_SLANT_NORMAL,
+                              CAIRO_FONT_WEIGHT_BOLD);
+      cairo_set_font_size(c,fontsize*1.25);
+      cairo_font_extents(c, &fextents);
+      cairo_text_extents(c, title, &extents);
+      cairo_move_to(c,pic_w/2-extents.width/2,y);
+      cairo_show_text(c,title);
+      y+=fextents.height;
+      cairo_restore(c);
+    }
+    cairo_font_extents(c, &fextents);
+    if(subtitle1){
+      cairo_text_extents(c, subtitle1, &extents);
+      cairo_move_to(c,pic_w/2-extents.width/2,y);
+      cairo_show_text(c,subtitle1);
+      y+=fextents.height;
+    }
+    if(subtitle2){
+      cairo_text_extents(c, subtitle2, &extents);
+      cairo_move_to(c,pic_w/2-extents.width/2,y);
+      cairo_show_text(c,subtitle2);
+      y+=fextents.height;
+    }
+    if(subtitle3){
+      cairo_text_extents(c, subtitle3, &extents);
+      cairo_move_to(c,pic_w/2-extents.width/2,y);
+      cairo_show_text(c,subtitle3);
+      y+=fextents.height;
+    }
+  }
+
+  /* color legend */
+  {
+    float cw;
+    cairo_text_extents(c, "100", &extents);
+    cw=extents.width*2*11;
+    cairo_text_extents(c, ".000001", &extents);
+    cw=MAX(cw,extents.width*1.5*7);
+    cairo_text_extents(c, ".0001%", &extents);
+    cw=MAX(cw,extents.width*1.5*7);
+
+    if(datatype==DT_iterations){
+      char buf[80];
+      int i;
+      float w=cw/11,px=leftpad+x_n;
+
+      for(i=0;i<=100;i++){
+        switch(i){
+        case 0:
+        case 1:
+        case 2:
+        case 3:
+        case 4:
+        case 5:
+        case 10:
+        case 20:
+        case 30:
+        case 40:
+        case 50:
+        case 100:
+          px-=w;
+          cairo_rectangle(c,px,legendy+toppad,w,legendh*1.25);
+          cairo_set_source_rgb(c,0,0,0);
+          cairo_stroke_preserve(c);
+          set_iter_color(c, i, 1.);
+          cairo_fill(c);
+
+          snprintf(buf,80,"%d",i);
+          cairo_text_extents(c, buf, &extents);
+          cairo_move_to(c,px+w/2-extents.width*.5,legendy+toppad+legendh*.625+extents.height/2);
+          set_iter_text_color(c, i);
+          cairo_show_text(c, buf);
+
+          break;
+        }
+      }
+      cairo_text_extents(c, legend_label, &extents);
+      cairo_move_to(c,px-extents.width-legendh*.75,toppad+legendy+legendh*.625+extents.height*.5);
+      cairo_show_text(c,legend_label);
+    }else{
+      int per_p = (datatype==DT_percent);
+      int i;
+      float w=cw/7,px=leftpad+x_n;
+
+      for(i=0;i<7;i++){
+        char *buf;
+
+        px-=w;
+        cairo_rectangle(c,px,legendy+toppad,w,legendh*1.25);
+        cairo_set_source_rgb(c,0,0,0);
+        cairo_stroke_preserve(c);
+
+        switch(i){
+        case 0:
+          buf=(per_p?".0001%":".000001");
+          set_error_color(c, 0., 1.);
+          cairo_fill(c);
+          set_error_text_color(c, 0.);
+          break;
+        case 1:
+          buf=(per_p?".001%":".00001");
+          set_error_color(c, .00001, 1.);
+          cairo_fill(c);
+          set_error_text_color(c, .00001);
+          break;
+        case 2:
+          buf=(per_p?".01%":".0001");
+          set_error_color(c, .0001, 1.);
+          cairo_fill(c);
+          set_error_text_color(c, .0001);
+          break;
+        case 3:
+          buf=(per_p?".1%":".001");
+          set_error_color(c, .001, 1.);
+          cairo_fill(c);
+          set_error_text_color(c, .001);
+          break;
+        case 4:
+          buf=(per_p?"1%":".01");
+          set_error_color(c, .01, 1.);
+          cairo_fill(c);
+          set_error_text_color(c, .01);
+          break;
+        case 5:
+          buf=(per_p?"10%":".1");
+          set_error_color(c, .1, 1.);
+          cairo_fill(c);
+          set_error_text_color(c, .1);
+          break;
+        case 6:
+          buf=(per_p?"100%":"1");
+          set_error_color(c, 1., 1.);
+          cairo_fill(c);
+          set_error_text_color(c, 1.);
+          break;
+        }
+
+        cairo_text_extents(c, buf, &extents);
+        cairo_move_to(c,px+w/2-extents.width*.5,legendy+toppad+legendh*.625+extents.height/2);
+        cairo_show_text(c,buf);
+
+      }
+      cairo_text_extents(c, legend_label, &extents);
+      cairo_move_to(c,px-extents.width-legendh*.75,toppad+legendy+legendh*.625+extents.height*.5);
+      cairo_show_text(c,legend_label);
+    }
+  }
+
+  return c;
+}
+
+void to_png(cairo_t *c,char *base, char *name){
+  if(c){
+    cairo_surface_t *cs=cairo_get_target(c);
+    char buf[320];
+    snprintf(buf,320,"%s-%s.png",base,name);
+    cairo_surface_write_to_png(cs,buf);
+  }
+}
+
+void destroy_page(cairo_t *c){
+  if(c){
+    cairo_surface_t *cs=cairo_get_target(c);
+    cairo_destroy(c);
+    cairo_surface_destroy(cs);
+  }
+}

Copied: trunk/chirptest/chirpgraph.h (from rev 17943, trunk/ghost/monty/chirp/chirpgraph.h)
===================================================================
--- trunk/chirptest/chirpgraph.h	                        (rev 0)
+++ trunk/chirptest/chirpgraph.h	2011-05-03 20:34:45 UTC (rev 17952)
@@ -0,0 +1,48 @@
+/********************************************************************
+ *                                                                  *
+ * THIS FILE IS PART OF THE OggGhost SOFTWARE CODEC SOURCE CODE.    *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
+ * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
+ * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
+ *                                                                  *
+ * THE OggGhost SOURCE CODE IS (C) COPYRIGHT 2007-2011              *
+ * by the Xiph.Org Foundation http://www.xiph.org/                  *
+ *                                                                  *
+ ********************************************************************
+
+ function: research-grade chirp extraction code
+ last mod: $Id$
+
+ ********************************************************************/
+
+#include <cairo/cairo.h>
+#define DT_iterations 0
+#define DT_abserror 1
+#define DT_percent 2
+extern int toppad;
+extern int leftpad;
+
+extern void set_error_color(cairo_t *c, float err,float a);
+extern void set_iter_color(cairo_t *cC, int ret, float a);
+extern void to_png(cairo_t *c,char *base, char *name);
+extern void destroy_page(cairo_t *c);
+extern cairo_t *draw_page(char *title,
+                          char *subtitle1,
+                          char *subtitle2,
+                          char *subtitle3,
+                          char *xaxis_label,
+                          char *yaxis_label,
+                          char *legend_label,
+                          int datatype);
+extern void setup_graphs(int start_x_step,
+                         int end_x_step, /* inclusive; not one past */
+                         int x_major_d,
+                         float x_major_f,
+
+                         int start_y_step,
+                         int end_y_step, /* inclusive; not one past */
+                         int y_major_d,
+                         float y_major_f,
+
+                         int subtitles,
+                         float fs);

Copied: trunk/chirptest/chirptest.c (from rev 17948, trunk/ghost/monty/chirp/chirptest.c)
===================================================================
--- trunk/chirptest/chirptest.c	                        (rev 0)
+++ trunk/chirptest/chirptest.c	2011-05-03 20:34:45 UTC (rev 17952)
@@ -0,0 +1,1308 @@
+/********************************************************************
+ *                                                                  *
+ * THIS FILE IS PART OF THE OggGhost SOFTWARE CODEC SOURCE CODE.    *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
+ * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
+ * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
+ *                                                                  *
+ * THE OggGhost SOURCE CODE IS (C) COPYRIGHT 2007-2011              *
+ * by the Xiph.Org Foundation http://www.xiph.org/                  *
+ *                                                                  *
+ ********************************************************************
+
+ function: graphing code for chirp tests
+ last mod: $Id$
+
+ ********************************************************************/
+
+#define _GNU_SOURCE
+#include <math.h>
+#include "chirp.h"
+#include "chirpgraph.h"
+#include "scales.h"
+#include "window.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <cairo/cairo.h>
+#include <pthread.h>
+#include <sys/time.h>
+
+float circular_distance(float A,float B){
+  float ret = A-B;
+  while(ret<-M_PI)ret+=2*M_PI;
+  while(ret>=M_PI)ret-=2*M_PI;
+  return ret;
+}
+
+#define DIM_ESTIMATE_A 1
+#define DIM_ESTIMATE_P 2
+#define DIM_ESTIMATE_W 3
+#define DIM_ESTIMATE_dA 4
+#define DIM_ESTIMATE_dW 5
+#define DIM_ESTIMATE_ddA 6
+#define DIM_ESTIMATE_MASK 0xf
+
+#define DIM_CHIRP_A (1<<4)
+#define DIM_CHIRP_P (2<<4)
+#define DIM_CHIRP_W (3<<4)
+#define DIM_CHIRP_dA (4<<4)
+#define DIM_CHIRP_dW  (5<<4)
+#define DIM_CHIRP_ddA (6<<4)
+#define DIM_CHIRP_MASK 0xf0
+
+/* ranges are inclusive */
+void set_chirp(chirp *c,
+               int xi, int xn, int xdim,
+               int yi, int yn, int ydim,
+               int stepi, int stepn, int rand_p,
+               float A0, float A1,
+               float P0, float P1,
+               float W0, float W1,
+               float dA0, float dA1,
+               float dW0, float dW1,
+               float ddA0, float ddA1){
+
+  float An,Pn,Wn,dAn,dWn,ddAn;
+  float Ai,Pi,Wi,dAi,dWi,ddAi;
+
+  xdim = (xdim | (xdim>>4)) & DIM_ESTIMATE_MASK;
+  ydim = (ydim | (ydim>>4)) & DIM_ESTIMATE_MASK;
+  if(stepn<2)stepn=2;
+
+  An=Pn=Wn=dAn=dWn=ddAn=stepn-1;
+  Ai = (rand_p ? drand48()*An : stepi);
+  Pi = (rand_p ? drand48()*Pn : stepi);
+  Wi = (rand_p ? drand48()*Wn : stepi);
+  dAi = (rand_p ? drand48()*dAn : stepi);
+  dWi = (rand_p ? drand48()*dWn : stepi);
+  ddAi = (rand_p ? drand48()*ddAn : stepi);
+
+  switch(xdim){
+  case DIM_ESTIMATE_A:
+    An = xn-1;
+    Ai = xi;
+    break;
+  case DIM_ESTIMATE_P:
+    Pn = xn-1;
+    Pi = xi;
+    break;
+  case DIM_ESTIMATE_W:
+    Wn = xn-1;
+    Wi = xi;
+    break;
+  case DIM_ESTIMATE_dA:
+    dAn = xn-1;
+    dAi = xi;
+    break;
+  case DIM_ESTIMATE_dW:
+    dWn = xn-1;
+    dWi = xi;
+    break;
+  case DIM_ESTIMATE_ddA:
+    ddAn = xn-1;
+    ddAi = xi;
+    break;
+  }
+
+  switch(ydim){
+  case DIM_ESTIMATE_A:
+    An = yn-1;
+    Ai = yi;
+    break;
+  case DIM_ESTIMATE_P:
+    Pn = yn-1;
+    Pi = yi;
+    break;
+  case DIM_ESTIMATE_W:
+    Wn = yn-1;
+    Wi = yi;
+    break;
+  case DIM_ESTIMATE_dA:
+    dAn = yn-1;
+    dAi = yi;
+    break;
+  case DIM_ESTIMATE_dW:
+    dWn = yn-1;
+    dWi = yi;
+    break;
+  case DIM_ESTIMATE_ddA:
+    ddAn = yn-1;
+    ddAi = yi;
+    break;
+  }
+
+  c->A = A0 + (A1-A0) / An * Ai;
+  c->P = P0 + (P1-P0) / Pn * Pi;
+  c->W = W0 + (W1-W0) / Wn * Wi;
+  c->dA = dA0 + (dA1-dA0) / dAn * dAi;
+  c->dW = dW0 + (dW1-dW0) / dWn * dWi;
+  c->ddA = ddA0 + (ddA1-ddA0) / ddAn * ddAi;
+
+}
+
+/*********************** Plot single estimate vs. chirp **********************/
+
+typedef struct {
+  float *in;
+  float *window;
+  int blocksize;
+  int max_iterations;
+  float fit_tolerance;
+
+  int fit_gauss_seidel;
+  int fit_W;
+  int fit_dA;
+  int fit_dW;
+  int fit_ddA;
+  int fit_nonlinear; /* 0==linear, 1==W-recentered, 2==W,dW-recentered */
+  float fit_W_alpha;
+  float fit_dW_alpha;
+  int fit_symm_norm;
+  int fit_bound_zero;
+
+  chirp *chirp;
+  chirp *estimate;
+  float *rms_error;
+  int *iterations;
+
+} colarg;
+
+
+pthread_mutex_t ymutex = PTHREAD_MUTEX_INITIALIZER;
+int next_y=0;
+int max_y=0;
+
+#define _GNU_SOURCE
+#include <fenv.h>
+
+void *compute_column(void *in){
+  colarg *arg = (colarg *)in;
+  int blocksize=arg->blocksize;
+  int y,i,ret;
+  int except;
+  int localinit = !arg->in;
+  float *chirp = localinit ? malloc(sizeof(*chirp)*blocksize) : arg->in;
+
+  while(1){
+    float rms_acc=0.;
+    float e_acc=0.;
+
+    pthread_mutex_lock(&ymutex);
+    y=next_y;
+    if(y>=max_y){
+      pthread_mutex_unlock(&ymutex);
+      return NULL;
+    }
+    next_y++;
+    pthread_mutex_unlock(&ymutex);
+
+    /* if the input is uninitialized, it's because we're sweeping or
+       randomizing chirp components across the column; generate the
+       input here in the thread */
+    if(localinit){
+      for(i=0;i<blocksize;i++){
+        double jj = i - blocksize/2 + .5;
+        double A = arg->chirp[y].A + (arg->chirp[y].dA + arg->chirp[y].ddA*jj)*jj;
+        double P = arg->chirp[y].P + (arg->chirp[y].W  + arg->chirp[y].dW *jj)*jj;
+        chirp[i] = A*cos(P);
+      }
+    }
+
+    except=feenableexcept(FE_ALL_EXCEPT);
+    fedisableexcept(FE_INEXACT);
+    fedisableexcept(FE_UNDERFLOW);
+
+    ret=estimate_chirps(chirp,arg->window,blocksize,
+                        arg->estimate+y,1,arg->fit_tolerance,arg->max_iterations,
+                        arg->fit_gauss_seidel,
+                        arg->fit_W,
+                        arg->fit_dA,
+                        arg->fit_dW,
+                        arg->fit_ddA,
+                        arg->fit_nonlinear,
+                        arg->fit_W_alpha,
+                        arg->fit_dW_alpha,
+                        arg->fit_symm_norm,
+                        arg->fit_bound_zero);
+
+    for(i=0;i<blocksize;i++){
+      double jj = i - blocksize/2 + .5;
+      double A = arg->estimate[y].A + (arg->estimate[y].dA + arg->estimate[y].ddA*jj)*jj;
+      double P = arg->estimate[y].P + (arg->estimate[y].W  + arg->estimate[y].dW *jj)*jj;
+      float  r = A*cos(P);
+      float ce = chirp[i]*arg->window[i];
+      float ee = (chirp[i]-r)*arg->window[i];
+
+      e_acc += ce*ce;
+      rms_acc += ee*ee;
+    }
+    arg->rms_error[y] = sqrt(rms_acc)/sqrt(e_acc);
+    arg->iterations[y] = arg->max_iterations-ret;
+    feclearexcept(FE_ALL_EXCEPT);
+    feenableexcept(except);
+
+  }
+
+  if(localinit)free(chirp);
+  return NULL;
+}
+
+typedef struct {
+  float fontsize;
+  char *subtitle1;
+  char *subtitle2;
+  char *subtitle3;
+  char *xaxis_label;
+  char *yaxis_label;
+
+  int blocksize;
+  int threads;
+
+  void (*window)(float *,int n);
+  float fit_tolerance;
+
+  int fit_gauss_seidel;
+  int fit_W;
+  int fit_dA;
+  int fit_dW;
+  int fit_ddA;
+  int fit_nonlinear;
+  float fit_W_alpha;
+  float fit_dW_alpha;
+  int fit_symm_norm;
+  int fit_bound_zero;
+
+  int x_dim;
+  int x_steps;
+  float x_major;
+  float x_minor;
+  int y_dim;
+  int y_steps;
+  float y_major;
+  float y_minor;
+  int sweep_steps;
+
+  /* If the randomize flag is unset and min!=max, a param is swept
+     from min to max (inclusive) divided into <sweep_steps>
+     increments.  If the rand flag is set, <sweep_steps> random values
+     in the range min to max instead. */
+  int sweep_or_rand_p;
+
+  float min_est_A;
+  float max_est_A;
+  int rel_est_A;
+
+  float min_est_P;
+  float max_est_P;
+  int rel_est_P;
+
+  float min_est_W;
+  float max_est_W;
+  int rel_est_W;
+
+  float min_est_dA;
+  float max_est_dA;
+  int rel_est_dA;
+
+  float min_est_dW;
+  float max_est_dW;
+  int rel_est_dW;
+
+  float min_est_ddA;
+  float max_est_ddA;
+  int rel_est_ddA;
+
+  float min_chirp_A;
+  float max_chirp_A;
+
+  float min_chirp_P;
+  float max_chirp_P;
+
+  float min_chirp_W;
+  float max_chirp_W;
+
+  float min_chirp_dA;
+  float max_chirp_dA;
+
+  float min_chirp_dW;
+  float max_chirp_dW;
+
+  float min_chirp_ddA;
+  float max_chirp_ddA;
+
+  /* generate which graphs? */
+  int graph_convergence_max;
+  int graph_convergence_delta;
+
+  int graph_Aerror_max;
+  int graph_Aerror_delta;
+
+  int graph_Perror_max;
+  int graph_Perror_delta;
+
+  int graph_Werror_max;
+  int graph_Werror_delta;
+
+  int graph_dAerror_max;
+  int graph_dAerror_delta;
+
+  int graph_dWerror_max;
+  int graph_dWerror_delta;
+
+  int graph_ddAerror_max;
+  int graph_ddAerror_delta;
+
+  int graph_RMSerror_max;
+  int graph_RMSerror_delta;
+
+}  graph_run;
+
+/* performs a W initial estimate error vs chirp W plot.  Ignores the
+   est and chirp arguments for W; these are pulled from the x and y setup */
+
+void w_e(char *filebase,graph_run *arg){
+  int threads=arg->threads;
+  int blocksize = arg->blocksize;
+  float window[blocksize];
+  float in[blocksize];
+  int i,xi,yi;
+
+  int x_n = arg->x_steps;
+  int y_n = arg->y_steps;
+
+  /* graphs:
+
+     convergence
+     Aerror
+     Perror
+     Werror
+     dAerror
+     dWerror
+     ddAerror
+     rms fit error
+
+     generate for:
+     worst case across sweep (0-7)
+     delta across sweep (8-15)
+  */
+
+  cairo_t *cC=NULL;
+  cairo_t *cA=NULL;
+  cairo_t *cP=NULL;
+  cairo_t *cW=NULL;
+  cairo_t *cdA=NULL;
+  cairo_t *cdW=NULL;
+  cairo_t *cddA=NULL;
+  cairo_t *cRMS=NULL;
+
+  cairo_t *cC_d=NULL;
+  cairo_t *cA_d=NULL;
+  cairo_t *cP_d=NULL;
+  cairo_t *cW_d=NULL;
+  cairo_t *cdA_d=NULL;
+  cairo_t *cdW_d=NULL;
+  cairo_t *cddA_d=NULL;
+  cairo_t *cRMS_d=NULL;
+
+  pthread_t threadlist[threads];
+  colarg targ[threads];
+
+  int swept=0;
+  int chirp_swept=0;
+  int est_swept=0;
+
+  float ret_minA[y_n];
+  float ret_minP[y_n];
+  float ret_minW[y_n];
+  float ret_mindA[y_n];
+  float ret_mindW[y_n];
+  float ret_minddA[y_n];
+  float ret_minRMS[y_n];
+  float ret_miniter[y_n];
+
+  float ret_maxA[y_n];
+  float ret_maxP[y_n];
+  float ret_maxW[y_n];
+  float ret_maxdA[y_n];
+  float ret_maxdW[y_n];
+  float ret_maxddA[y_n];
+  float ret_maxRMS[y_n];
+  float ret_maxiter[y_n];
+
+  float minX=0,maxX=0;
+  float minY=0,maxY=0;
+  int x0s,x1s;
+  int y0s,y1s;
+  int xmajori,ymajori;
+  int xminori,yminori;
+
+  struct timeval last;
+  gettimeofday(&last,NULL);
+
+  switch(arg->x_dim){
+  case DIM_ESTIMATE_A:
+    minX = arg->min_est_A;
+    maxX = arg->max_est_A;
+    break;
+  case DIM_ESTIMATE_P:
+    minX = arg->min_est_P;
+    maxX = arg->max_est_P;
+    break;
+  case DIM_ESTIMATE_W:
+    minX = arg->min_est_W;
+    maxX = arg->max_est_W;
+    break;
+  case DIM_ESTIMATE_dA:
+    minX = arg->min_est_dA;
+    maxX = arg->max_est_dA;
+    break;
+  case DIM_ESTIMATE_dW:
+    minX = arg->min_est_dW;
+    maxX = arg->max_est_dW;
+    break;
+  case DIM_ESTIMATE_ddA:
+    minX = arg->min_est_ddA;
+    maxX = arg->max_est_ddA;
+    break;
+  case DIM_CHIRP_A:
+    minX = arg->min_chirp_A;
+    maxX = arg->max_chirp_A;
+    break;
+  case DIM_CHIRP_P:
+    minX = arg->min_chirp_P;
+    maxX = arg->max_chirp_P;
+    break;
+  case DIM_CHIRP_W:
+    minX = arg->min_chirp_W;
+    maxX = arg->max_chirp_W;
+    break;
+  case DIM_CHIRP_dA:
+    minX = arg->min_chirp_dA;
+    maxX = arg->max_chirp_dA;
+    break;
+  case DIM_CHIRP_dW:
+    minX = arg->min_chirp_dW;
+    maxX = arg->max_chirp_dW;
+    break;
+  case DIM_CHIRP_ddA:
+    minX = arg->min_chirp_ddA;
+    maxX = arg->max_chirp_ddA;
+    break;
+  }
+
+  switch(arg->y_dim){
+  case DIM_ESTIMATE_A:
+    minY = arg->min_est_A;
+    maxY = arg->max_est_A;
+    break;
+  case DIM_ESTIMATE_P:
+    minY = arg->min_est_P;
+    maxY = arg->max_est_P;
+    break;
+  case DIM_ESTIMATE_W:
+    minY = arg->min_est_W;
+    maxY = arg->max_est_W;
+    break;
+  case DIM_ESTIMATE_dA:
+    minY = arg->min_est_dA;
+    maxY = arg->max_est_dA;
+    break;
+  case DIM_ESTIMATE_dW:
+    minY = arg->min_est_dW;
+    maxY = arg->max_est_dW;
+    break;
+  case DIM_ESTIMATE_ddA:
+    minY = arg->min_est_ddA;
+    maxY = arg->max_est_ddA;
+    break;
+  case DIM_CHIRP_A:
+    minY = arg->min_chirp_A;
+    maxY = arg->max_chirp_A;
+    break;
+  case DIM_CHIRP_P:
+    minY = arg->min_chirp_P;
+    maxY = arg->max_chirp_P;
+    break;
+  case DIM_CHIRP_W:
+    minY = arg->min_chirp_W;
+    maxY = arg->max_chirp_W;
+    break;
+  case DIM_CHIRP_dA:
+    minY = arg->min_chirp_dA;
+    maxY = arg->max_chirp_dA;
+    break;
+  case DIM_CHIRP_dW:
+    minY = arg->min_chirp_dW;
+    maxY = arg->max_chirp_dW;
+    break;
+  case DIM_CHIRP_ddA:
+    minY = arg->min_chirp_ddA;
+    maxY = arg->max_chirp_ddA;
+    break;
+  }
+
+  x0s = rint((x_n-1)/(maxX-minX)*minX);
+  y0s = rint((y_n-1)/(maxY-minY)*minY);
+  x1s = x0s+x_n-1;
+  y1s = y0s+y_n-1;
+
+  xminori = rint((x_n-1)/(maxX-minX)*arg->x_minor);
+  yminori = rint((y_n-1)/(maxY-minY)*arg->y_minor);
+  xmajori = rint(xminori/arg->x_minor*arg->x_major);
+  ymajori = rint(yminori/arg->y_minor*arg->y_major);
+
+  if(xminori<1 || yminori<1 || xmajori<1 || ymajori<1){
+    fprintf(stderr,"Bad xmajor/xminor/ymajor/yminor value.\n");
+    exit(1);
+  }
+
+  if( rint(xmajori*(maxX-minX)/arg->x_major) != x_n-1){
+    fprintf(stderr,"Xmajor or Xminor results in non-integer pixel increment.\n");
+    exit(1);
+  }
+
+  if( rint(ymajori*(maxY-minY)/arg->y_major) != y_n-1){
+    fprintf(stderr,"Ymajor or Yminor results in non-integer pixel increment.\n");
+    exit(1);
+  }
+
+  /* determine ~ padding needed */
+  setup_graphs(x0s,x1s,xmajori,arg->x_major,
+               y0s,y1s,ymajori,arg->y_major,
+               (arg->subtitle1!=0)+(arg->subtitle2!=0)+(arg->subtitle3!=0),
+               arg->fontsize);
+
+  if(arg->sweep_steps>1){
+    if(!(arg->x_dim==DIM_ESTIMATE_A || arg->y_dim==DIM_ESTIMATE_A) &&
+       arg->min_est_A != arg->max_est_A) est_swept=1;
+    if(!(arg->x_dim==DIM_ESTIMATE_P || arg->y_dim==DIM_ESTIMATE_P) &&
+       arg->min_est_P != arg->max_est_P) est_swept=1;
+    if(!(arg->x_dim==DIM_ESTIMATE_W || arg->y_dim==DIM_ESTIMATE_W) &&
+       arg->min_est_W != arg->max_est_W) est_swept=1;
+    if(!(arg->x_dim==DIM_ESTIMATE_dA || arg->y_dim==DIM_ESTIMATE_dA) &&
+       arg->min_est_dA != arg->max_est_dA) est_swept=1;
+    if(!(arg->x_dim==DIM_ESTIMATE_dW || arg->y_dim==DIM_ESTIMATE_dW) &&
+       arg->min_est_dW != arg->max_est_dW) est_swept=1;
+    if(!(arg->x_dim==DIM_ESTIMATE_ddA || arg->y_dim==DIM_ESTIMATE_ddA) &&
+       arg->min_est_ddA != arg->max_est_ddA) est_swept=1;
+
+    if(!(arg->x_dim==DIM_CHIRP_A || arg->y_dim==DIM_CHIRP_A) &&
+       arg->min_chirp_A != arg->max_chirp_A) chirp_swept=1;
+    if(!(arg->x_dim==DIM_CHIRP_P || arg->y_dim==DIM_CHIRP_P) &&
+       arg->min_chirp_P != arg->max_chirp_P) chirp_swept=1;
+    if(!(arg->x_dim==DIM_CHIRP_W || arg->y_dim==DIM_CHIRP_W) &&
+       arg->min_chirp_W != arg->max_chirp_W) chirp_swept=1;
+    if(!(arg->x_dim==DIM_CHIRP_dA || arg->y_dim==DIM_CHIRP_dA) &&
+       arg->min_chirp_dA != arg->max_chirp_dA) chirp_swept=1;
+    if(!(arg->x_dim==DIM_CHIRP_dW || arg->y_dim==DIM_CHIRP_dW) &&
+       arg->min_chirp_dW != arg->max_chirp_dW) chirp_swept=1;
+    if(!(arg->x_dim==DIM_CHIRP_ddA || arg->y_dim==DIM_CHIRP_ddA) &&
+       arg->min_chirp_ddA != arg->max_chirp_ddA) chirp_swept=1;
+  }
+
+  swept = est_swept | chirp_swept;
+
+  if(arg->y_dim==DIM_CHIRP_A &&
+     arg->min_chirp_A != arg->max_chirp_A) chirp_swept=1;
+  if(arg->y_dim==DIM_CHIRP_P &&
+     arg->min_chirp_P != arg->max_chirp_P) chirp_swept=1;
+  if(arg->y_dim==DIM_CHIRP_W &&
+     arg->min_chirp_W != arg->max_chirp_W) chirp_swept=1;
+  if(arg->y_dim==DIM_CHIRP_dA &&
+     arg->min_chirp_dA != arg->max_chirp_dA) chirp_swept=1;
+  if(arg->y_dim==DIM_CHIRP_dW &&
+     arg->min_chirp_dW != arg->max_chirp_dW) chirp_swept=1;
+  if(arg->y_dim==DIM_CHIRP_ddA &&
+     arg->min_chirp_ddA != arg->max_chirp_ddA) chirp_swept=1;
+
+  if(arg->graph_convergence_max)
+    cC = draw_page(!swept?"Convergence":"Worst Case Convergence",
+                   arg->subtitle1,
+                   arg->subtitle2,
+                   arg->subtitle3,
+                   arg->xaxis_label,
+                   arg->yaxis_label,
+                   "Iterations:",
+                   DT_iterations);
+  if(swept && arg->graph_convergence_delta)
+    cC_d = draw_page("Convergence Delta",
+                     arg->subtitle1,
+                     arg->subtitle2,
+                     arg->subtitle3,
+                     arg->xaxis_label,
+                     arg->yaxis_label,
+                     "Iteration span:",
+                     DT_iterations);
+
+  if(arg->graph_Aerror_max)
+    cA = draw_page(!swept?"A (Amplitude) Error":"Maximum A (Amplitude) Error",
+                   arg->subtitle1,
+                   arg->subtitle2,
+                   arg->subtitle3,
+                   arg->xaxis_label,
+                   arg->yaxis_label,
+                   "Percentage Error:",
+                   DT_percent);
+  if(swept && arg->graph_Aerror_delta)
+    cA_d = draw_page("A (Amplitude) Delta",
+                     arg->subtitle1,
+                     arg->subtitle2,
+                     arg->subtitle3,
+                     arg->xaxis_label,
+                     arg->yaxis_label,
+                     "Delta:",
+                     DT_abserror);
+
+  if(arg->graph_Perror_max)
+    cP = draw_page(!swept?"P (Phase) Error":"Maximum P (Phase) Error",
+                   arg->subtitle1,
+                   arg->subtitle2,
+                   arg->subtitle3,
+                   arg->xaxis_label,
+                   arg->yaxis_label,
+                   "Error (radians):",
+                   DT_abserror);
+  if(swept && arg->graph_Perror_delta)
+    cP_d = draw_page("Phase Delta",
+                     arg->subtitle1,
+                     arg->subtitle2,
+                     arg->subtitle3,
+                     arg->xaxis_label,
+                     arg->yaxis_label,
+                     "Delta (radians):",
+                     DT_abserror);
+
+  if(arg->fit_W){
+    if(arg->graph_Werror_max)
+      cW = draw_page(!swept?"W (Frequency) Error":"Maximum W (Frequency) Error",
+                     arg->subtitle1,
+                     arg->subtitle2,
+                     arg->subtitle3,
+                     arg->xaxis_label,
+                     arg->yaxis_label,
+                     "Error (cycles/block):",
+                     DT_abserror);
+    if(swept && arg->graph_Werror_delta)
+      cW_d = draw_page("Frequency Delta",
+                       arg->subtitle1,
+                       arg->subtitle2,
+                       arg->subtitle3,
+                       arg->xaxis_label,
+                       arg->yaxis_label,
+                       "Delta (cycles/block):",
+                       DT_abserror);
+  }
+
+  if(arg->fit_dA){
+    if(arg->graph_dAerror_max)
+      cdA = draw_page(!swept?"dA (Amplitude Modulation) Error":"Maximum dA (Amplitude Modulation) Error",
+                      arg->subtitle1,
+                      arg->subtitle2,
+                      arg->subtitle3,
+                      arg->xaxis_label,
+                      arg->yaxis_label,
+                      "Error:",
+                      DT_abserror);
+    if(swept && arg->graph_dAerror_delta)
+      cdA_d = draw_page("Amplitude Modulation Delta",
+                        arg->subtitle1,
+                        arg->subtitle2,
+                        arg->subtitle3,
+                        arg->xaxis_label,
+                        arg->yaxis_label,
+                        "Delta:",
+                        DT_abserror);
+  }
+
+  if(arg->fit_dW){
+    if(arg->graph_dWerror_max)
+      cdW = draw_page(!swept?"dW (Chirp Rate) Error":"Maximum dW (Chirp Rate) Error",
+                      arg->subtitle1,
+                      arg->subtitle2,
+                      arg->subtitle3,
+                      arg->xaxis_label,
+                      arg->yaxis_label,
+                      "Error (cycles/block):",
+                      DT_abserror);
+    if(swept && arg->graph_dWerror_delta)
+      cdW_d = draw_page("Chirp Rate Delta",
+                        arg->subtitle1,
+                        arg->subtitle2,
+                        arg->subtitle3,
+                        arg->xaxis_label,
+                        arg->yaxis_label,
+                        "Delta (cycles/block):",
+                        DT_abserror);
+  }
+
+  if(arg->fit_ddA){
+    if(arg->graph_ddAerror_max)
+      cddA = draw_page(!swept?"ddA (Amplitude Modulation Squared) Error":
+                       "Maximum ddA (Amplitude Modulation Squared) Error",
+                       arg->subtitle1,
+                       arg->subtitle2,
+                       arg->subtitle3,
+                       arg->xaxis_label,
+                       arg->yaxis_label,
+                       "Error:",
+                       DT_abserror);
+    if(swept && arg->graph_ddAerror_delta)
+      cddA_d = draw_page("Amplitude Modulation Squared Delta",
+                         arg->subtitle1,
+                         arg->subtitle2,
+                         arg->subtitle3,
+                         arg->xaxis_label,
+                         arg->yaxis_label,
+                         "Delta:",
+                         DT_abserror);
+  }
+
+  if(arg->graph_RMSerror_max)
+    cRMS = draw_page(swept?"Maximum RMS Fit Error":
+                     "RMS Fit Error",
+                     arg->subtitle1,
+                     arg->subtitle2,
+                     arg->subtitle3,
+                     arg->xaxis_label,
+                     arg->yaxis_label,
+                     "Percentage Error:",
+                     DT_percent);
+  if(swept && arg->graph_RMSerror_delta)
+    cRMS_d = draw_page("RMS Error Delta",
+                       arg->subtitle1,
+                       arg->subtitle2,
+                       arg->subtitle3,
+                       arg->xaxis_label,
+                       arg->yaxis_label,
+                       "Percentage Delta:",
+                       DT_percent);
+
+  if(arg->window)
+    arg->window(window,blocksize);
+  else
+    for(i=0;i<blocksize;i++)
+      window[i]=1.;
+
+  /* graph computation */
+  for(xi=0;xi<x_n;xi++){
+    chirp chirps[y_n];
+    chirp estimates[y_n];
+    int iter[y_n];
+    float rms[y_n];
+    int si,sn=(swept && arg->sweep_steps>1 ? arg->sweep_steps : 1);
+
+    fprintf(stderr,"\r%s: column %d/%d...",filebase,xi,x_n-1);
+
+    memset(targ,0,sizeof(targ));
+
+    for(i=0;i<threads;i++){
+      if(!chirp_swept)targ[i].in=in;
+      targ[i].window=window;
+      targ[i].blocksize=blocksize;
+      targ[i].max_iterations=100;
+      targ[i].fit_tolerance=arg->fit_tolerance;
+      targ[i].fit_gauss_seidel=arg->fit_gauss_seidel;
+      targ[i].fit_W=arg->fit_W;
+      targ[i].fit_dA=arg->fit_dA;
+      targ[i].fit_dW=arg->fit_dW;
+      targ[i].fit_ddA=arg->fit_ddA;
+      targ[i].fit_nonlinear=arg->fit_nonlinear;
+      targ[i].fit_W_alpha=arg->fit_W_alpha;
+      targ[i].fit_dW_alpha=arg->fit_dW_alpha;
+      targ[i].fit_symm_norm=arg->fit_symm_norm;
+      targ[i].fit_bound_zero=arg->fit_bound_zero;
+      targ[i].chirp=chirps;
+      //targ[i].in=NULL;
+      targ[i].estimate=estimates;
+      targ[i].rms_error=rms;
+      targ[i].iterations=iter;
+    }
+
+    /* if we're sweeping a parameter, we're going to iterate here for a bit. */
+    for(si=0;si<sn;si++){
+      max_y=y_n;
+      next_y=0;
+
+      /* compute/set chirp and est parameters, potentially compute the
+         chirp waveform */
+      for(yi=0;yi<y_n;yi++){
+
+        set_chirp(chirps+yi,
+                  xi,x_n,
+                  arg->x_dim&DIM_CHIRP_MASK,
+                  y_n-yi-1,y_n,
+                  arg->y_dim&DIM_CHIRP_MASK,
+                  si,sn,
+                  arg->sweep_or_rand_p,
+                  arg->min_chirp_A,
+                  arg->max_chirp_A,
+                  arg->min_chirp_P*2.*M_PI,
+                  arg->max_chirp_P*2.*M_PI,
+                  arg->min_chirp_W*2.*M_PI/blocksize,
+                  arg->max_chirp_W*2.*M_PI/blocksize,
+                  arg->min_chirp_dA/blocksize,
+                  arg->max_chirp_dA/blocksize,
+                  arg->min_chirp_dW*2.*M_PI/blocksize/blocksize,
+                  arg->max_chirp_dW*2.*M_PI/blocksize/blocksize,
+                  arg->min_chirp_ddA/blocksize/blocksize,
+                  arg->max_chirp_ddA/blocksize/blocksize);
+
+        set_chirp(estimates+yi,
+                  xi,x_n,
+                  arg->x_dim&DIM_ESTIMATE_MASK,
+                  y_n-yi-1,y_n,
+                  arg->y_dim&DIM_ESTIMATE_MASK,
+                  si,sn,
+                  arg->sweep_or_rand_p,
+                  arg->min_est_A,
+                  arg->max_est_A,
+                  arg->min_est_P*2.*M_PI,
+                  arg->max_est_P*2.*M_PI,
+                  arg->min_est_W*2.*M_PI/blocksize,
+                  arg->max_est_W*2.*M_PI/blocksize,
+                  arg->min_est_dA/blocksize,
+                  arg->max_est_dA/blocksize,
+                  arg->min_est_dW*2.*M_PI/blocksize/blocksize,
+                  arg->max_est_dW*2.*M_PI/blocksize/blocksize,
+                  arg->min_est_ddA/blocksize/blocksize,
+                  arg->max_est_ddA/blocksize/blocksize);
+
+        if(arg->rel_est_A) estimates[yi].A += chirps[yi].A;
+        if(arg->rel_est_P) estimates[yi].P += chirps[yi].P;
+        if(arg->rel_est_W) estimates[yi].W += chirps[yi].W;
+        if(arg->rel_est_dA) estimates[yi].dA += chirps[yi].dA;
+        if(arg->rel_est_dW) estimates[yi].dW += chirps[yi].dW;
+        if(arg->rel_est_ddA) estimates[yi].ddA += chirps[yi].ddA;
+
+      }
+
+      if(!chirp_swept){
+        for(i=0;i<threads;i++)
+          targ[i].in=in;
+
+        for(i=0;i<blocksize;i++){
+          double jj = i-blocksize/2+.5;
+          double A = chirps[0].A + (chirps[0].dA + chirps[0].ddA*jj)*jj;
+          double P = chirps[0].P + (chirps[0].W  + chirps[0].dW *jj)*jj;
+          in[i] = A*cos(P);
+        }
+      }
+
+      /* compute column */
+      for(i=0;i<threads;i++)
+        pthread_create(threadlist+i,NULL,compute_column,targ+i);
+      for(i=0;i<threads;i++)
+        pthread_join(threadlist[i],NULL);
+
+      /* accumulate results for this pass */
+      if(si==0){
+        for(i=0;i<y_n;i++){
+          ret_minA[i]=ret_maxA[i]=chirps[i].A - estimates[i].A;
+          ret_minP[i]=ret_maxP[i]=circular_distance(chirps[i].P,estimates[i].P);
+          ret_minW[i]=ret_maxW[i]=chirps[i].W - estimates[i].W;
+          ret_mindA[i]=ret_maxdA[i]=chirps[i].dA - estimates[i].dA;
+          ret_mindW[i]=ret_maxdW[i]=chirps[i].dW - estimates[i].dW;
+          ret_minddA[i]=ret_maxddA[i]=chirps[i].ddA - estimates[i].ddA;
+          ret_minRMS[i]=ret_maxRMS[i]=rms[i];
+          ret_miniter[i]=ret_maxiter[i]=iter[i];
+        }
+      }else{
+        for(i=0;i<y_n;i++){
+          float v = chirps[i].A - estimates[i].A;
+          if(ret_minA[i]>v)ret_minA[i]=v;
+          if(ret_maxA[i]<v)ret_maxA[i]=v;
+
+          v = circular_distance(chirps[i].P, estimates[i].P);
+          if(ret_minP[i]>v)ret_minP[i]=v;
+          if(ret_maxP[i]<v)ret_maxP[i]=v;
+
+          v = chirps[i].W - estimates[i].W;
+          if(ret_minW[i]>v)ret_minW[i]=v;
+          if(ret_maxW[i]<v)ret_maxW[i]=v;
+
+          v = chirps[i].dA - estimates[i].dA;
+          if(ret_mindA[i]>v)ret_mindA[i]=v;
+          if(ret_maxdA[i]<v)ret_maxdA[i]=v;
+
+          v = chirps[i].dW - estimates[i].dW;
+          if(ret_mindW[i]>v)ret_mindW[i]=v;
+          if(ret_maxdW[i]<v)ret_maxdW[i]=v;
+
+          v = chirps[i].ddA - estimates[i].ddA;
+          if(ret_minddA[i]>v)ret_minddA[i]=v;
+          if(ret_maxddA[i]<v)ret_maxddA[i]=v;
+
+          if(ret_minRMS[i]>rms[i])ret_minRMS[i]=rms[i];
+          if(ret_maxRMS[i]<rms[i])ret_maxRMS[i]=rms[i];
+
+          if(ret_miniter[i]>iter[i])ret_miniter[i]=iter[i];
+          if(ret_maxiter[i]<iter[i])ret_maxiter[i]=iter[i];
+        }
+      }
+    }
+
+    /* Column sweep complete; plot */
+    for(yi=0;yi<y_n;yi++){
+      int x=x0s+xi;
+      int y=y1s-yi;
+      float a = 1.;
+      if(x%xminori==0 || y%yminori==0) a = .8;
+      if(x%xmajori==0 || y%ymajori==0) a = .3;
+
+      /* Convergence graph */
+      if(cC){
+        set_iter_color(cC,ret_maxiter[yi],a);
+        cairo_rectangle(cC,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cC);
+      }
+
+      /* Convergence delta graph */
+      if(cC_d){
+        set_iter_color(cC_d,ret_maxiter[yi]-ret_miniter[yi],a);
+        cairo_rectangle(cC_d,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cC_d);
+      }
+
+      /* A error graph */
+      if(cA){
+        set_error_color(cA,MAX(fabs(ret_maxA[yi]),fabs(ret_minA[yi])),a);
+        cairo_rectangle(cA,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cA);
+      }
+
+      /* A delta graph */
+      if(cA_d){
+        set_error_color(cA_d,ret_maxA[yi]-ret_minA[yi],a);
+        cairo_rectangle(cA_d,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cA_d);
+      }
+
+      /* P error graph */
+      if(cP){
+        set_error_color(cP,MAX(fabs(ret_maxP[yi]),fabs(ret_minP[yi])),a);
+        cairo_rectangle(cP,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cP);
+      }
+
+      /* P delta graph */
+      if(cP_d){
+        set_error_color(cP_d,circular_distance(ret_maxP[yi],ret_minP[yi]),a);
+        cairo_rectangle(cP_d,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cP_d);
+      }
+
+      if(cW){
+        /* W error graph */
+        set_error_color(cW,MAX(fabs(ret_maxW[yi]),fabs(ret_minW[yi]))/2./M_PI*blocksize,a);
+        cairo_rectangle(cW,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cW);
+      }
+
+        /* W delta graph */
+      if(cW_d){
+        set_error_color(cW_d,(ret_maxW[yi]-ret_minW[yi])/2./M_PI*blocksize,a);
+        cairo_rectangle(cW_d,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cW_d);
+      }
+
+      if(cdA){
+        /* dA error graph */
+        set_error_color(cdA,MAX(fabs(ret_maxdA[yi]),fabs(ret_mindA[yi]))*blocksize,a);
+        cairo_rectangle(cdA,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cdA);
+      }
+        /* dA delta graph */
+      if(cdA_d){
+        set_error_color(cdA_d,(ret_maxdA[yi]-ret_mindA[yi])*blocksize,a);
+        cairo_rectangle(cdA_d,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cdA_d);
+      }
+
+      if(cdW){
+        /* dW error graph */
+        set_error_color(cdW,MAX(fabs(ret_maxdW[yi]),fabs(ret_mindW[yi]))/2./M_PI*blocksize*blocksize,a);
+        cairo_rectangle(cdW,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cdW);
+      }
+
+      /* dW delta graph */
+      if(cdW_d){
+        set_error_color(cdW_d,(ret_maxdW[yi]-ret_mindW[yi])/2./M_PI*blocksize*blocksize,a);
+        cairo_rectangle(cdW_d,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cdW_d);
+      }
+
+      if(cddA){
+        /* ddA error graph */
+        set_error_color(cddA,MAX(fabs(ret_maxddA[yi]),fabs(ret_minddA[yi]))*blocksize*blocksize,a);
+        cairo_rectangle(cddA,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cddA);
+      }
+
+      /* dA delta graph */
+      if(cddA_d){
+        set_error_color(cddA_d,(ret_maxddA[yi]-ret_minddA[yi])*blocksize*blocksize,a);
+        cairo_rectangle(cddA_d,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cddA_d);
+      }
+
+      /* RMS error graph */
+      if(cRMS){
+        set_error_color(cRMS,ret_maxRMS[yi],a);
+        cairo_rectangle(cRMS,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cRMS);
+      }
+
+      /* RMS delta graph */
+      if(cRMS_d){
+        set_error_color(cRMS_d,ret_maxRMS[yi]-ret_minRMS[yi],a);
+        cairo_rectangle(cRMS_d,xi+leftpad,yi+toppad,1,1);
+        cairo_fill(cRMS_d);
+      }
+    }
+
+    /* dump a graph update every 10 seconds */
+    {
+      struct timeval now;
+      gettimeofday(&now,NULL);
+      if(now.tv_sec-last.tv_sec + (now.tv_usec-last.tv_usec)/1000000. > 10.){
+        last=now;
+
+        to_png(cC,filebase,"converge");
+        to_png(cA,filebase,"Aerror");
+        to_png(cP,filebase,"Perror");
+        to_png(cW,filebase,"Werror");
+        to_png(cdA,filebase,"dAerror");
+        to_png(cdW,filebase,"dWerror");
+        to_png(cddA,filebase,"ddAerror");
+        to_png(cRMS,filebase,"RMSerror");
+
+        to_png(cC_d,filebase,"convdelta");
+        to_png(cA_d,filebase,"Adelta");
+        to_png(cP_d,filebase,"Pdelta");
+        to_png(cW_d,filebase,"Wdelta");
+        to_png(cdA_d,filebase,"dAdelta");
+        to_png(cdW_d,filebase,"dWdelta");
+        to_png(cddA_d,filebase,"ddAdelta");
+        to_png(cRMS_d,filebase,"RMSdelta");
+      }
+    }
+  }
+
+  to_png(cC,filebase,"converge");
+  to_png(cA,filebase,"Aerror");
+  to_png(cP,filebase,"Perror");
+  to_png(cW,filebase,"Werror");
+  to_png(cdA,filebase,"dAerror");
+  to_png(cdW,filebase,"dWerror");
+  to_png(cddA,filebase,"ddAerror");
+  to_png(cRMS,filebase,"RMSerror");
+
+  to_png(cC_d,filebase,"convdelta");
+  to_png(cA_d,filebase,"Adelta");
+  to_png(cP_d,filebase,"Pdelta");
+  to_png(cW_d,filebase,"Wdelta");
+  to_png(cdA_d,filebase,"dAdelta");
+  to_png(cdW_d,filebase,"dWdelta");
+  to_png(cddA_d,filebase,"ddAdelta");
+  to_png(cRMS_d,filebase,"RMSdelta");
+
+  fprintf(stderr," done\n");
+}
+
+int main(){
+  graph_run arg={
+    /* fontsize */      18,
+    /* subtitle1 */     "Linear estimation, no ddA fit",
+    /* subtitle2 */     "chirp: A=1.0, dA=0., swept phase | estimate A=P=dA=dW=0, estimate W=chirp W",
+    /* subtitle3 */     "sine window",
+    /* xaxis label */   "W (cycles/block)",
+    /* yaxis label */   "dW (cycles/block)",
+
+    /* blocksize */     128,
+    /* threads */       32,
+
+    /* window */        window_functions.sine,
+    /* fit_tol */       .000001,
+    /* gauss_seidel */  1,
+    /* fit_W */         1,
+    /* fit_dA */        1,
+    /* fit_dW */        1,
+    /* fit_ddA */       0,
+    /* nonlinear */     0,
+    /* W_alpha */       1.,
+    /* dW_alpha */      1.,
+    /* symm_norm */     0,
+    /* bound_zero */    0,
+
+    /* x dimension */   DIM_CHIRP_W,
+    /* x steps */       1001,
+    /* x major */       1.,
+    /* x minor */       .25,
+    /* y dimension */   DIM_CHIRP_dW,
+    /* y steps */       601,
+    /* y major */       1.,
+    /* y minor */       .25,
+    /* sweep_steps */   32,
+    /* randomize_p */   0,
+
+    /* est A range */     0.,  0.,  0, /* relative flag */
+    /* est P range */     0.,  0.,  0, /* relative flag */
+    /* est W range */     0.,  0.,  1, /* relative flag */
+    /* est dA range */    0.,  0.,  0, /* relative flag */
+    /* est dW range */    0.,  0.,  0, /* relative flag */
+    /* est ddA range */   0.,  0.,  0, /* relative flag */
+
+    /* ch A range */    1.,1.,
+    /* ch P range */    0.,1.-1./32.,
+    /* ch W range */    0.,10.,
+    /* ch dA range */   0.,0.,
+    /* ch dW range */   -2.5,2.5,
+    /* ch ddA range */  0.,0.,
+
+    /* converge max */    1,
+    /* converge del */    0,
+    /* max A error */     1,
+    /* A error delta */   0,
+    /* max P error */     1,
+    /* P error delta */   0,
+    /* max W error */     1,
+    /* W error delta */   0,
+    /* max dA error */    1,
+    /* dA error delta */  0,
+    /* max dW error */    1,
+    /* dW error delta */  0,
+    /* max ddA error */   1,
+    /* ddA error delta */ 0,
+    /* max RMS error */   1,
+    /* RMS error delta */ 0,
+
+  };
+
+  /* Graphs for dW vs W ****************************************/
+
+  //w_e("linear-dW-vs-W",&arg);
+  arg.fit_nonlinear=1;
+  arg.subtitle1="Partial nonlinear estimation, no ddA fit";
+  //w_e("partial-nonlinear-dW-vs-W",&arg);
+  arg.subtitle1="Full nonlinear estimation, no ddA fit";
+  arg.fit_nonlinear=2;
+  //w_e("full-nonlinear-dW-vs-W",&arg);
+
+  /* Graphs for W estimate distance vs W ************************/
+
+  arg.subtitle1="Linear estimation, no ddA fit";
+  arg.fit_nonlinear=0;
+  arg.yaxis_label="initial distance from W (cycles/block)";
+  arg.y_dim = DIM_ESTIMATE_W;
+  arg.min_est_W = -2.5;
+  arg.max_est_W =  2.5;
+  arg.min_chirp_dW=0.;
+  arg.max_chirp_dW=0.;
+
+  //w_e("linear-estW-vs-W",&arg);
+  arg.subtitle1="Partial nonlinear estimation, no ddA fit";
+  arg.subtitle2="chirp: A=1.0, dA=dW=0., swept phase | estimate A=P=dA=dW=0";
+  arg.fit_nonlinear=1;
+  //w_e("partial-nonlinear-estW-vs-W",&arg);
+  arg.subtitle1="Full nonlinear estimation, no ddA fit";
+  arg.fit_nonlinear=2;
+  //w_e("full-nonlinear-estW-vs-W",&arg);
+  arg.fit_nonlinear=0;
+
+  /* graphs for different windows */
+
+  arg.min_est_W = -2.5;
+  arg.max_est_W =  2.5;
+  arg.max_chirp_W =  10;
+  arg.fit_nonlinear = 0;
+  arg.subtitle1="Linear estimation, no ddA fit";
+  arg.window = window_functions.rectangle;
+  arg.subtitle3 = "rectangular window";
+  //w_e("linear-estW-vs-W-rectangular",&arg);
+
+  arg.window = window_functions.sine;
+  arg.subtitle3 = "sine window";
+  //w_e("linear-estW-vs-W-sine",&arg);
+
+  arg.window = window_functions.hanning;
+  arg.subtitle3 = "hanning window";
+  //w_e("linear-estW-vs-W-hanning",&arg);
+
+  arg.window = window_functions.tgauss_deep;
+  arg.subtitle3 = "unimodal triangular/gaussian window";
+  //w_e("linear-estW-vs-W-unimodal",&arg);
+
+  arg.window = window_functions.maxwell1;
+  arg.subtitle3 = "maxwell (optimized) window";
+  //w_e("linear-estW-vs-W-maxwell",&arg);
+
+  arg.min_est_W = -15;
+  arg.max_est_W =  15;
+  arg.max_chirp_W =  25;
+  arg.fit_nonlinear = 2;
+  arg.subtitle1="Fully nonlinear estimation, no ddA fit";
+  arg.window = window_functions.rectangle;
+  arg.subtitle3 = "rectangular window";
+  //w_e("full-nonlinear-estW-vs-W-rectangular",&arg);
+
+  arg.window = window_functions.sine;
+  arg.subtitle3 = "sine window";
+  //w_e("full-nonlinear-estW-vs-W-sine",&arg);
+
+  arg.window = window_functions.hanning;
+  arg.subtitle3 = "hanning window";
+  //w_e("full-nonlinear-estW-vs-W-hanning",&arg);
+
+  arg.window = window_functions.tgauss_deep;
+  arg.subtitle3 = "unimodal triangular/gaussian window";
+  //w_e("full-nonlinear-estW-vs-W-unimodal",&arg);
+
+  arg.window = window_functions.maxwell1;
+  arg.subtitle3 = "maxwell (optimized) window";
+  w_e("full-nonlinear-estW-vs-W-maxwell",&arg);
+
+  /* 1, 1.5, 2nd order */
+  arg.min_est_W = -3;
+  arg.max_est_W =  3;
+  arg.max_chirp_W =  10;
+  arg.window = window_functions.hanning;
+  arg.subtitle3 = "hanning window";
+  arg.fit_dW = 0;
+  arg.fit_dA = 0;
+
+  arg.fit_nonlinear = 2;
+  arg.subtitle1="Fully nonlinear estimation, half-order fit (W only)";
+  //w_e("nonlinear-estW-vs-W-.5order",&arg);
+
+  arg.fit_dA=1;
+
+  arg.fit_nonlinear = 0;
+  arg.subtitle1="Linear estimation, first-order fit";
+  //w_e("linear-estW-vs-W-1order",&arg);
+
+  arg.fit_nonlinear = 2;
+  arg.subtitle1="Fully nonlinear estimation, first-order fit";
+  //w_e("nonlinear-estW-vs-W-1order",&arg);
+
+  arg.fit_dW=1;
+
+  arg.fit_nonlinear = 2;
+  arg.subtitle1="Fully nonlinear estimation, 1.5th-order fit (dW only)";
+  //w_e("nonlinear-estW-vs-W-1.5order",&arg);
+
+  arg.fit_ddA=1;
+
+  arg.fit_nonlinear = 0;
+  arg.subtitle1="Linear estimation, second-order fit";
+  //w_e("linear-estW-vs-W-2order",&arg);
+
+  arg.fit_nonlinear = 2;
+  arg.subtitle1="Fully nonlinear estimation, second-order fit";
+  //w_e("nonlinear-estW-vs-W-2order",&arg);
+
+
+  return 0;
+}
+

Copied: trunk/chirptest/configure.ac (from rev 13251, trunk/ghost/configure.ac)
===================================================================
--- trunk/chirptest/configure.ac	                        (rev 0)
+++ trunk/chirptest/configure.ac	2011-05-03 20:34:45 UTC (rev 17952)
@@ -0,0 +1,77 @@
+dnl Process this file with autoconf to produce a configure script. -*-m4-*-
+
+cflags_save="$CFLAGS"
+AC_PREREQ(2.57)
+AC_INIT(chirptest.c)
+AM_INIT_AUTOMAKE(chirp, 20110503, [ghost-dev at xiph.org])
+AC_CONFIG_FILES([Makefile])
+m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])], [AC_SUBST([AM_DEFAULT_VERBOSITY], [1])])
+
+# Checks for programs.
+AC_PROG_CC
+
+# Checks for libraries.
+AC_CHECK_LIB([m], [sincos])
+AC_CHECK_LIB([pthread], [pthread_create])
+AC_CHECK_LIB([cairo], [cairo_create])
+
+# Checks for header files.
+AC_HEADER_STDC
+AC_CHECK_HEADERS([stdlib.h string.h stdio.h])
+
+# Checks for typedefs, structures, and compiler characteristics.
+AC_C_CONST
+
+if test -z "$GCC"; then
+        case $host in 
+        *-*-irix*)
+                DEBUG="-g -signed" 
+                CFLAGS="-O2 -w -signed"
+                PROFILE="-p -g3 -O2 -signed"
+                ;;
+        sparc-sun-solaris*)
+                DEBUG="-v -g"
+                CFLAGS="-xO4 -fast -w -fsimple -native -xcg92"
+                PROFILE="-v -xpg -g -xO4 -fast -native -fsimple -xcg92 -Dsuncc"
+                ;;
+        *)
+                DEBUG="-g"
+                CFLAGS="-O"
+                PROFILE="-g -p" 
+                ;;
+        esac
+else
+        case $host in 
+        *-*-linux*)
+                DEBUG="-g -Wall -fsigned-char"
+                CFLAGS="-O2 -fsigned-char -ffast-math"
+                PROFILE="-Wall -W -pg -g -O2 -fsigned-char -ffast-math"
+                ;;
+        sparc-sun-*)
+                DEBUG="-g -Wall -fsigned-char"
+                CFLAGS="-O2 -fsigned-char"
+                PROFILE="-pg -g -O2 -fsigned-char" 
+                ;;
+        *-*-darwin*)
+                DEBUG="-fno-common -g -Wall -fsigned-char"
+                CFLAGS="-fno-common -O2 -Wall -fsigned-char -ffast-math"
+                PROFILE="-fno-common -O2 -Wall -pg -g -fsigned-char -ffast-math"
+                ;;
+        *)
+                DEBUG="-g -Wall -fsigned-char"
+                CFLAGS="-O2 -fsigned-char -ffast-math"
+                PROFILE="-O2 -g -pg -fsigned-char -ffast-math" 
+                ;;
+        esac
+fi
+
+CFLAGS="$CFLAGS -DVERSION='\"$VERSION\"' $COMMON_FLAGS"
+DEBUG="$DEBUG -DVERSION='\\\"$VERSION\\\"' $COMMON_FLAGS"
+PROFILE="$PROFILE -DVERSION='\\\"$VERSION\\\"' $COMMON_FLAGS"
+LIBS="$LIBS"
+
+AC_SUBST(DEBUG)
+AC_SUBST(PROFILE)
+AC_OUTPUT
+
+echo "Type \"make\" to compile";

Copied: trunk/chirptest/scales.h (from rev 17943, trunk/ghost/monty/chirp/scales.h)
===================================================================
--- trunk/chirptest/scales.h	                        (rev 0)
+++ trunk/chirptest/scales.h	2011-05-03 20:34:45 UTC (rev 17952)
@@ -0,0 +1,38 @@
+/********************************************************************
+ *                                                                  *
+ * THIS FILE IS PART OF THE OggGhost SOFTWARE CODEC SOURCE CODE.    *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
+ * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
+ * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
+ *                                                                  *
+ * THE OggGhost SOURCE CODE IS (C) COPYRIGHT 2007                   *
+ * by the Xiph.Org Foundation http://www.xiph.org/                  *
+ *                                                                  *
+ ********************************************************************
+
+ function: scale inlines
+ last mod: $Id$
+
+ ********************************************************************/
+
+#include <math.h>
+
+#define MIN(a,b)  ((a)<(b) ? (a):(b))
+#define MAX(a,b)  ((a)>(b) ? (a):(b))
+
+static inline float todB(float x){
+  return log(x*x+1e-24f)*4.34294480f;
+}
+
+static inline float fromdB(float x){
+  return exp((x)*.11512925f);
+}
+
+static inline float toBark(float x){
+  return 13.1f*atan(.00074f*x) + 2.24f*atan(x*x*1.85e-8f) + 1e-4f*x;
+}
+
+static inline float fromBark(float x){
+  return 102.f*x - 2.f*pow(x,2.f) + .4f*pow(x,3.f) + pow(1.46f,x) - 1.f;
+}
+

Copied: trunk/chirptest/window.c (from rev 17943, trunk/ghost/monty/chirp/window.c)
===================================================================
--- trunk/chirptest/window.c	                        (rev 0)
+++ trunk/chirptest/window.c	2011-05-03 20:34:45 UTC (rev 17952)
@@ -0,0 +1,141 @@
+/********************************************************************
+ *                                                                  *
+ * THIS FILE IS PART OF THE OggGhost SOFTWARE CODEC SOURCE CODE.    *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
+ * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
+ * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
+ *                                                                  *
+ * THE OggGhost SOURCE CODE IS (C) COPYRIGHT 2007-2011              *
+ * by the Xiph.Org Foundation http://www.xiph.org/                  *
+ *                                                                  *
+ ********************************************************************
+
+ function: window functions for research code
+ last mod: $Id$
+
+ ********************************************************************/
+
+#define _GNU_SOURCE
+#include <math.h>
+#include "window.h"
+
+/* A few windows */
+
+/* Not-a-window */
+static void rectangle(float *x, int n){
+  int i;
+  for(i=0;i<n;i++)
+    x[i]=1.;
+}
+
+/* sine window */
+static void sine(float *x, int n){
+  int i;
+  float scale = M_PI/n;
+
+  for(i=0;i<n;i++){
+    float i5 = i+.5;
+    x[i] = sin(scale*i5);
+  }
+}
+
+/* Minimum 4-term Blackman Harris; highest sidelobe = -92dB */
+#define A0 .35875f
+#define A1 .48829f
+#define A2 .14128f
+#define A3 .01168f
+
+static void blackmann_harris(float *x, int n){
+  int i;
+  float scale = 2*M_PI/n;
+
+  for(i=0;i<n;i++){
+    float i5 = i+.5;
+    float w = A0 - A1*cos(scale*i5) + A2*cos(scale*i5*2) - A3*cos(scale*i5*3);
+    x[i] = w;
+  }
+}
+
+/* Good 'ol Hanning window (narrow mainlobe, fast falloff, high sidelobes) */
+static void hanning(float *x, int n){
+  int i;
+  float scale = 2*M_PI/n;
+
+  for(i=0;i<n;i++){
+    float i5 = i+.5;
+    x[i] = (.5-.5*cos(scale*i5));
+  }
+}
+
+/* Triangular * gaussian (sidelobeless window) */
+#define TGA 1.e-4
+#define TGB 21.6
+static void tgauss_deep(float *x, int n){
+  int i;
+  for(i=0;i<n;i++){
+    float f = (i-n/2.)/(n/2.);
+    x[i] = exp(-TGB*pow(f,2)) * pow(1.-fabs(f),TGA);
+  }
+}
+
+static void vorbis(float *d, int n){
+  int i;
+  for(i=0;i<n;i++)
+    d[i] = sin(0.5 * M_PI * sin((i+.5)/n * M_PI)*sin((i+.5)/n * M_PI));
+}
+
+static float beta(int n, float alpha){
+  return cosh (acosh(pow(10,alpha))/(n-1));
+}
+
+static double T(double n, double x){
+  if(fabs(x)<=1){
+    return cos(n*acos(x));
+  }else{
+    return cosh(n*acosh(x));
+  }
+}
+
+/* Dolph-Chebyshev window (a=6., all sidelobes < -120dB) */
+static void dolphcheb(float *d, int n){
+  int i,k;
+  float a = 6.;
+  int M=n/2;
+  int N=M*2;
+  double b = beta(N,a);
+
+  for(i=0;i<n;i++){
+    double sum=0;
+    for(k=0;k<M;k++)
+      sum += (k&1?-1:1)*T(N,b*cos(M_PI*k/N)) * cos (2*i*k*M_PI/N);
+    sum /= T(N,b);
+    sum-=.5;
+    d[i]=sum;
+  }
+}
+
+/* sidelobeless window machine optimized for fast convergence */
+static void maxwell1(float *x, int n){
+  int i;
+  float scale = 2*M_PI/n;
+  for(i=0;i<n;i++){
+    float i5 = i+.5;
+    x[i] = pow( 119.72981
+                - 119.24098*cos(scale*i5)
+                + 0.10283622*cos(2*scale*i5)
+                + 0.044013144*cos(3*scale*i5)
+                + 0.97203713*cos(4*scale*i5),
+                1.9730763)/50000.;
+  }
+}
+
+window_bundle window_functions = {
+  rectangle,
+  sine,
+  hanning,
+  vorbis,
+  blackmann_harris,
+  tgauss_deep,
+  dolphcheb,
+  maxwell1,
+};

Copied: trunk/chirptest/window.h (from rev 17943, trunk/ghost/monty/chirp/window.h)
===================================================================
--- trunk/chirptest/window.h	                        (rev 0)
+++ trunk/chirptest/window.h	2011-05-03 20:34:45 UTC (rev 17952)
@@ -0,0 +1,29 @@
+/********************************************************************
+ *                                                                  *
+ * THIS FILE IS PART OF THE OggGhost SOFTWARE CODEC SOURCE CODE.    *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS     *
+ * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
+ * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING.       *
+ *                                                                  *
+ * THE OggGhost SOURCE CODE IS (C) COPYRIGHT 2007-2011              *
+ * by the Xiph.Org Foundation http://www.xiph.org/                  *
+ *                                                                  *
+ ********************************************************************
+
+ function: window functions for research code
+ last mod: $Id$
+
+ ********************************************************************/
+
+typedef struct {
+  void (*rectangle)(float *,int);
+  void (*sine)(float *,int);
+  void (*hanning)(float *,int);
+  void (*vorbis)(float *,int);
+  void (*blackman_harris)(float *,int);
+  void (*tgauss_deep)(float *,int);
+  void (*dolphcheb)(float *,int);
+  void (*maxwell1)(float *,int);
+} window_bundle;
+
+extern window_bundle window_functions;



More information about the commits mailing list