[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