[xiph-commits] r17895 - trunk/ghost/monty/chirp

xiphmont at svn.xiph.org xiphmont at svn.xiph.org
Thu Mar 17 23:06:42 PDT 2011


Author: xiphmont
Date: 2011-03-17 23:06:42 -0700 (Thu, 17 Mar 2011)
New Revision: 17895

Added:
   trunk/ghost/monty/chirp/bessel.c
   trunk/ghost/monty/chirp/chirp.c
   trunk/ghost/monty/chirp/chirp.h
   trunk/ghost/monty/chirp/harness.c
   trunk/ghost/monty/chirp/scales.h
   trunk/ghost/monty/chirp/smallft.c
   trunk/ghost/monty/chirp/smallft.h
   trunk/ghost/monty/chirp/window.c
Log:
continue moving around chirp estimation code



Copied: trunk/ghost/monty/chirp/bessel.c (from rev 17893, trunk/ghost/monty/sinusoid/bessel.c)
===================================================================
--- trunk/ghost/monty/chirp/bessel.c	                        (rev 0)
+++ trunk/ghost/monty/chirp/bessel.c	2011-03-18 06:06:42 UTC (rev 17895)
@@ -0,0 +1,68 @@
+/********************************************************************
+ *                                                                  *
+ * 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: second order bessel filters
+ last mod: $Id$
+
+ ********************************************************************/
+
+#define _GNU_SOURCE
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+typedef struct {
+  float c;
+  float igain;
+} bessel1;
+
+typedef struct {
+  float c[2];
+  float igain;
+} bessel2;
+
+bessel2 mkbessel2(float alpha){
+  bessel2 ret={{0}};
+  float w = 2.f*tan(M_PIl * alpha);  
+  float r = w * -1.10160133059;
+  float i = w * .636009824757;
+  float mag = 1.f/(4.f + (r-4.f)*r + i*i);
+  float br = (4.f - r*r - i*i)*mag;
+  float bi = 4.f*i*mag;
+
+  ret.c[1] = -br*br - bi*bi;
+  ret.c[0] = 2.f*br;
+  ret.igain=(1.f-ret.c[0]-ret.c[1])*.25f;
+  return ret;
+}
+
+bessel1 mkbessel1(float alpha){
+  bessel1 ret={0};
+  float r = -2.f*tan(M_PIl * alpha);
+  float mag = 1.f/(4.f + (r-4.f)*r);
+  ret.c = (4.f - r*r)*mag;
+  ret.igain = (1.f-ret.c)*.5f;
+  return ret;
+}
+
+//static inline float compbessel2(float x2, float x1, float x0, float y2, float y1, bessel2 *b){
+//  return (x0+x1*2.+x2)*b->igain + y1*b->c1+y2*b->c2;
+//}
+
+static inline float bessel1_floor(float x1, float x0, float y1, bessel1 *a, bessel1 *b){
+  float ay = (x0+x1)*a->igain + y1*a->c;
+  //float by = (x0+x1)*b->igain + y1*b->c;
+  float by = x0;
+  return (ay<by ? ay : by);
+}
+

Copied: trunk/ghost/monty/chirp/chirp.c (from rev 17893, trunk/ghost/monty/sinusoid/chirp.c)
===================================================================
--- trunk/ghost/monty/chirp/chirp.c	                        (rev 0)
+++ trunk/ghost/monty/chirp/chirp.c	2011-03-18 06:06:42 UTC (rev 17895)
@@ -0,0 +1,278 @@
+/********************************************************************
+ *                                                                  *
+ * 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$
+
+ ********************************************************************/
+
+#define _GNU_SOURCE
+#include <math.h>
+#include "sinusoids.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);
+}
+
+/* performs a second-order nonlinear chirp fit. Passed in chirps are
+   used as initial estimation inputs to the iteration. Applies the
+   square of passed in window to the input/basis functions.
+
+   x: input signal (unwindowed)
+   y: output reconstruction (unwindowed)
+   window: window to apply to input/basis during fitting process
+   len: length of x and window vectors
+   c: chirp estimation inputs, chirp fit outputs (sorted in frequency order)
+   n: number of chirps
+   iter_limit: maximum allowable number of iterations
+   fit_limit: desired fit limit
+
+   Chirp frequency must be within ~ 1 FFT bin of fit frequency (see
+   paper).  Amplitude, amplitude' and phase are fit wouthout a need
+   for initial estimate.  Frequency' fit will converge from any point,
+   though it will fit to the closest lobe (wheter mainlobe or
+   sidelobe).  All parameters fit within 3 iterations with the
+   exception of frequency' which shows approximately linear
+   convergence outside of ~ .5 frequency bins.
+
+   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(float *x, float *y, float *window, int len,
+                    chirp *c, int n, int iter_limit, float fit_limit){
+
+  int i,j,flag=1;
+  float cwork[len];
+  float swork[len];
+  float residue[len];
+
+  /* Build a starting reconstruction based on the passied in initial
+     chirp estimates */
+  memset(y,0,sizeof(*y)*len);
+  for(i=0;i<n;i++){
+    for(j=0;j<len;j++){
+      float jj = j-len*.5+.5;
+      float a = c[i].A + c[i].dA*jj;
+      float w = (c[i].W + c[i].dW*jj)*jj;
+      y[j] += a*cos(w+c[i].P);
+    }
+  }
+
+  /* outer fit iteration */
+  while(flag && iter_limit--){
+    flag=0;
+
+    /* precompute the portion of the projection/fit estimate shared by
+       the zero, first and second order fits.  Subtractsthe current
+       best fits from the input signal and widows the result. */
+    for(j=0;j<len;j++)
+      residue[j]=(x[j]-y[j])*window[j];
+    memset(y,0,sizeof(*y)*len);
+
+    /* Sort chirps by descending amplitude */
+    qsort(c, n, sizeof(*c), cmp_descending_A);
+
+    /* reestimate the next chirp fit */
+    for(i=0;i<n;i++){
+      float lA2,lA,lW,ldA,lP,ldW;
+      float aP=0, bP=0;
+      float cP=0, dP=0;
+      float eP=0, fP=0;
+      float aE=0, bE=0;
+      float cE=0, dE=0;
+      float eE=0, fE=0;
+      float aC = cos(c[i].P);
+      float aS = sin(c[i].P);
+
+      for(j=0;j<len;j++){
+	float jj = j-len*.5+.5;
+	float jj2 = jj*jj;
+	float co = cos((c[i].W + c[i].dW*jj)*jj)*window[j];
+	float si = sin((c[i].W + c[i].dW*jj)*jj)*window[j];
+        float c2 = co*co;
+        float s2 = si*si;
+
+        /* add the current estimate back to the residue vector */
+        float yy = residue[j] += (aC*co-aS*si) * (c[i].A + c[i].dA*jj);
+
+        /* Cache the basis function premultiplied by jj for the two
+           later fitting loops */
+        cwork[j] = c2*jj;
+        swork[j] = s2*jj;
+
+        /* zero order projection */
+	aP += co*yy;
+	bP += si*yy;
+        /* the part of the first order fit we can perform before
+           having an updated amplitude/phase.  This saves needing to
+           cache or recomupute yy later.  */
+        cP += co*yy*jj;
+        dP += si*yy*jj;
+        /* the part of the second order fit we can perform before
+           having an updated amplitude/phase/frequency/ampliude'.
+           This saves needing to cache or recomupute yy later. */
+        eP += co*yy*jj2;
+        fP += si*yy*jj2;
+
+        /* integrate the total energy under the curve of each basis
+           function */
+	aE += c2;
+	bE += s2;
+	cE += c2*jj2;
+	dE += s2*jj2;
+        eE += c2*jj2*jj2;
+        fE += s2*jj2*jj2;
+      }
+      /* aP/bP in terms of total basis energy */
+      aP = (aE>0.f ? aP/aE : 0.f);
+      bP = (bE>0.f ? bP/bE : 0.f);
+
+      /* the frequency/amplitude' fits converge faster now that we
+         have an updated amplitude/phase estimate.  Finish the portion
+         of the first order projection we couldn't compute above. */
+      for(j=0;j<len;j++){
+        cP -= aP*cwork[j];
+        dP -= bP*swork[j];
+      }
+      /* cP/dP in terms of total basis energy */
+      cP = (cE>0.f ? cP/cE : 0.f);
+      dP = (dE>0.f ? dP/dE : 0.f);
+
+      /* the frequency' fit converges faster now that we have an
+         updated estimates for the other parameters.  Finish the
+         portion of the second order projection we couldn't compute
+         above. */
+      for(j=0;j<len;j++){
+        float jj = j-len*.5+.5;
+        eP -= (aP+cP*jj)*cwork[j]*jj;
+        fP -= (bP+dP*jj)*swork[j]*jj;
+      }
+      /* eP/fP in terms of total basis energy */
+      eP = (eE>0.f ? eP/eE : 0.f);
+      fP = (fE>0.f ? fP/fE : 0.f);
+
+      /* computer new chirp fit estimate from basis projections */
+      lA2 = aP*aP + bP*bP;
+      lA = sqrt(lA2);
+      lP = -atan2(bP, aP);
+      lW = (cP*bP - dP*aP)/lA2;
+      ldA = (cP*aP + dP*bP)/lA;
+      ldW = 1.75*(eP*bP - fP*aP)/lA2;
+
+      /* have we converged to within a requested limit */
+      if(c[i].A==0 || fabs((lA-c[i].A)/c[i].A)>fit_limit) flag=1;
+      if(fabs(lP - c[i].P)>fit_limit) flag=1;
+      if(fabs(lW*len)>fit_limit) flag=1;
+      if(c[i].A==0 || fabs((ldA-c[i].dA)/c[i].A)>fit_limit) flag=1;
+      if(fabs(ldW*len*len)>fit_limit*2.*M_PI) flag=1;
+
+      /* save new fit estimate */
+      c[i].A = lA;
+      c[i].P = lP;
+      c[i].W += lW;
+      c[i].dA = ldA;
+      c[i].dW += ldW;
+
+      /* update the reconstruction/residue vectors with new fit */
+      for(j=0;j<len;j++){
+        float jj = j-len*.5+.5;
+        float a = c[i].A + c[i].dA*jj;
+        float w = (c[i].W + c[i].dW*jj)*jj;
+        float v = a*cos(w+c[i].P);
+        residue[j] -= v*window[j];
+        y[j] += v;
+      }
+    }
+  }
+
+  /* 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].W += c[i].dW*len;
+    c[i].P += c[i].W*len;
+  }
+}
+
+/* OMG!  An example! */
+
+#if 0
+
+void hanning(float *x, int n){
+  float scale = 2*M_PI/n;
+  int i;
+  for(i=0;i<n;i++){
+    float i5 = i+.5;
+    x[i] = .5-.5*cos(scale*i5);
+  }
+}
+
+int main(){
+  int BLOCKSIZE=1024;
+  int i;
+  float w[1024];
+  float x[1024];
+  float y[1024];
+
+  chirp c[]={
+    {.3, 100.*2.*M_PI/BLOCKSIZE, .2, 0, 0},
+    {1., 95.*2.*M_PI/BLOCKSIZE, .2, 0, 2.5*2.*M_PI/BLOCKSIZE/BLOCKSIZE}
+  };
+
+  hanning(w,BLOCKSIZE);
+
+  for(i=0;i<BLOCKSIZE;i++){
+    float jj = i-BLOCKSIZE*.5+.5;
+    x[i]=(.3+.1*jj/BLOCKSIZE)*cos(jj*2.*M_PI/BLOCKSIZE * 100.2 +1.2);
+    x[i]+=cos(jj*2.*M_PI/BLOCKSIZE * (96. + 3./BLOCKSIZE*jj) +1.2);
+  }
+
+  int ret=estimate_chirps(x,y,w,BLOCKSIZE,c,2,55,.01);
+
+  for(i=0;i<sizeof(c)/sizeof(*c);i++)
+    fprintf(stderr,"W:%f\nP:%f\nA:%f\ndW:%f\ndA:%f\nreturned=%d\n\n",
+            c[i].W*BLOCKSIZE/2./M_PI,
+            c[i].P,
+            c[i].A,
+            c[i].dW*BLOCKSIZE*BLOCKSIZE/2./M_PI,
+            c[i].dA*BLOCKSIZE,ret);
+
+  return 0;
+}
+#endif

Copied: trunk/ghost/monty/chirp/chirp.h (from rev 17893, trunk/ghost/monty/sinusoid/chirp.h)
===================================================================
--- trunk/ghost/monty/chirp/chirp.h	                        (rev 0)
+++ trunk/ghost/monty/chirp/chirp.h	2011-03-18 06:06:42 UTC (rev 17895)
@@ -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: 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) */
+} chirp;
+
+extern int estimate_chirps(float *x, float *y, float *window, int len,
+                           chirp *c, int n, int iter_limit, float fit_limit);
+extern void advance_chirps(chirp *c, int n, int len);
+

Copied: trunk/ghost/monty/chirp/harness.c (from rev 17893, trunk/ghost/monty/sinusoid/harness.c)
===================================================================
--- trunk/ghost/monty/chirp/harness.c	                        (rev 0)
+++ trunk/ghost/monty/chirp/harness.c	2011-03-18 06:06:42 UTC (rev 17895)
@@ -0,0 +1,740 @@
+#define _GNU_SOURCE
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <cairo/cairo.h>
+#include <squishyio/squishyio.h>
+#include "sinusoids.h"
+#include "smallft.h"
+#include "scales.h"
+
+static int BLOCKSIZE=1024;
+
+static float MINdB=-100;
+static float MAXdB=0;
+
+#define AREAS 2
+static float AREA_SIZE[AREAS]={1.,1.};
+static int   AREA_HZOOM[AREAS]={1,1};
+static int   AREA_VZOOM[AREAS]={2,2};
+static int   AREA_FRAMEOUT_TRIGGER=1;
+
+static int   SPECTROGRAM_AREA=0;
+
+static int   SPECTRUM_AREA=-1;
+static int   SPECTRUM_GRIDFONTSIZE=8;
+
+static int   ENVELOPE_AREA=-1;
+
+static int   pic_w,pic_h;
+static drft_lookup fft;
+
+
+#define A0 .35875f
+#define A1 .48829f
+#define A2 .14128f
+#define A3 .01168f
+
+void spread_bark(float *in, float *out, int n, int rate,
+                 float loslope, float width, float hislope);
+
+
+void blackmann_harris(float *out, float *in, 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);
+    out[i] = in[i]*w;
+  }
+}
+
+static void hanning(float *out, float *in, int n){
+  int i;
+  float scale = 2*M_PI/n;
+
+  for(i=0;i<n;i++){
+    float i5 = i+.5;
+    out[i] = in[i]*(.5-.5*cos(scale*i5));
+  }
+}
+
+void mag_dB(float *log,float *d, int n){
+  int i;
+  log[0] = todB(d[0]*d[0])*.5;
+  for(i=2;i<n;i+=2)
+    log[i>>1] = todB(d[i-1]*d[i-1] + d[i]*d[i])*.5;
+  log[n/2] = todB(d[n-1]*d[n-1])*.5;
+}
+
+static int frameno=0;
+void write_yuv(pcm_t *pcm, cairo_surface_t *s){
+  int w=cairo_image_surface_get_width(s);
+  int h=cairo_image_surface_get_height(s);
+  unsigned char *d=cairo_image_surface_get_data(s);
+  unsigned char *p=malloc(w*h*3);
+  int i;
+  unsigned char *p1,*p2,*p3,*dd;
+
+  printf("FRAME S0 L%d P%.3f\n",w*h*3,
+         (float)frameno*(BLOCKSIZE/AREA_VZOOM[AREA_FRAMEOUT_TRIGGER])/
+         pcm->rate);
+  p1=p;
+  p2=p+w*h;
+  p3=p2+w*h;
+  dd=d;
+  for(i=0;i<w*h;i++){
+    *(p1++)=(65481*dd[2]+128553*dd[1]
+             +24966*dd[0]+4207500)/255000;
+    *(p2++)=(-8372*4*dd[2]-16436*4*dd[1]+
+             +24808*4*dd[0]+29032005)/225930;
+    *(p3++)=(39256*4*dd[2]-32872*4*dd[1]
+             -6384*4*dd[0]+45940035)/357510;
+    dd+=4;
+  }
+  fwrite(p,1,w*h*3,stdout);
+
+  free(p);
+  frameno++;
+}
+
+static off_t pcmno=0;
+void write_pcm(pcm_t *pcm,float **d, int pos, int samples){
+  int i,j;
+  unsigned char *p=malloc(samples*pcm->ch*3);
+  unsigned char *pp=p;
+  printf("FRAME S1 L%d P%.3f\n",samples*pcm->ch*3,(float)pcmno/pcm->rate);
+
+  for(i=pos;i<samples+pos;i++)
+    for(j=0;j<pcm->ch;j++){
+      int val = rint(d[j][i]*8388608.f);
+      if(val<-8388608)val=-8388608;
+      if(val>8388607)val=8388607;
+      p[0]=(val&0xff);
+      p[1]=((val>>8)&0xff);
+      p[2]=((val>>16)&0xff);
+      p+=3;
+    }
+  fwrite(pp,1,samples*pcm->ch*3,stdout);
+  free(pp);
+  pcmno+=samples;
+}
+
+int area_y(int area){
+  int i;
+  float n=0,d=0;
+  for(i=0;i<area;i++)
+    n+=AREA_SIZE[i];
+  for(;i<AREAS;i++)
+    d+=AREA_SIZE[i];
+  d+=n;
+  return (int)rint(n/d*pic_h);
+}
+
+int area_h(int area){
+  return area_y(area+1)-area_y(area);
+}
+
+void clip_area(cairo_t *c, int area){
+  cairo_rectangle(c,0,area_y(area),pic_w,area_h(area));
+  cairo_clip(c);
+}
+
+
+void clear_area(cairo_t *c, int area){
+  cairo_set_source_rgb(c,0,0,0);
+  cairo_rectangle(c,0,area_y(area),pic_w,area_h(area));
+  cairo_fill(c);
+}
+
+
+void scroll_area(cairo_t *c, int area){
+  /* shift image up by a pixel */
+  cairo_set_source_surface(c,cairo_get_target(c),0,-1);
+  cairo_paint(c);
+  cairo_set_source_rgb(c,0,0,0);
+  cairo_rectangle(c,0,area_y(area)+area_h(area)-1,pic_w,1);
+  cairo_fill(c);
+}
+
+void pick_a_color(int channel, float *r, float *g, float *b){
+
+  *r=*g=*b=0.;
+  switch(channel){
+  case 5:
+    *r=1.;
+    break;
+  case 1:
+    *g=1.;
+    break;
+  case 2:
+    *b=1.;
+    break;
+  case 3:
+    *r=.7;
+    *g=.7;
+    break;
+  case 4:
+    *g=.7;
+    *b=.7;
+    break;
+  case 0:
+    *r=.7;
+    *b=.7;
+    break;
+  case 6:
+    *r=.45;
+    *b=.45;
+    *g=.45;
+    break;
+  case 7:
+    *r=.8;
+    *g=.3;
+    break;
+  }
+}
+
+void spectrogram(pcm_t *pcm, cairo_t *c,off_t pos){
+  int i,j;
+  float r,g,b;
+  float work[BLOCKSIZE];
+  float work2[BLOCKSIZE];
+  float pc[BLOCKSIZE];
+  int y = area_y(SPECTROGRAM_AREA)+area_h(SPECTROGRAM_AREA)-1;
+  int hzoom = AREA_HZOOM[SPECTROGRAM_AREA];
+  cairo_save(c);
+  clip_area(c,SPECTROGRAM_AREA);
+  scroll_area(c,SPECTROGRAM_AREA);
+
+  cairo_set_operator(c,CAIRO_OPERATOR_ADD);
+
+  for(i=0;i<pcm->ch;i++){
+    if(pos+BLOCKSIZE<=pcm->samples)
+      memcpy(work, &pcm->data[i][pos],BLOCKSIZE*sizeof(*work));
+    else{
+      memset(work, 0, sizeof(work));
+      memcpy(work, &pcm->data[i][pos],(pcm->samples-pos)*sizeof(*work));
+    }
+
+    /* window and drft */
+    hanning(work, work,BLOCKSIZE);
+    drft_forward(&fft, work);
+    for(j=0;j<BLOCKSIZE;j++)
+      work[j]*=4./BLOCKSIZE;
+
+    mag_dB(work,work,BLOCKSIZE);
+    spread_bark(work, work2, BLOCKSIZE/2+1, pcm->rate, -20, 1., -20);
+    for(j=0;j<BLOCKSIZE/2;j++)
+    work[j] = work[j]-work2[j]-50.;
+    for(j=0;j<BLOCKSIZE/2;j++)
+    pc[j] = pc[j]-work2[j]-50.;
+
+    pick_a_color(i,&r,&g,&b);
+    for(j=0;j<BLOCKSIZE/2;j++){
+      float alpha = work[j]*(1./-MINdB)+1.f;
+      if(alpha<0.)alpha=0.;
+      if(alpha>1.)alpha=1.;
+      cairo_set_source_rgba(c,r,g,b,alpha);
+      cairo_rectangle(c,j*hzoom-hzoom/2,y,hzoom,1);
+      cairo_fill(c);
+    }
+  }
+
+  cairo_restore(c);
+}
+
+void grid_lines(cairo_t *c,int area){
+  int y1 = area_y(area);
+  int h = area_h(area);
+  int y2 = y1+h;
+  int j;
+
+  cairo_save(c);
+  clip_area(c,area);
+  cairo_set_line_width(c,1.);
+  cairo_set_source_rgb(c,.2,.2,.2);
+  for(j=0;j<pic_w;j+=AREA_HZOOM[area]){
+    cairo_move_to(c,j+.5,y1);
+    cairo_line_to(c,j+.5,y2);
+  }
+  for(j=MINdB;j<=MAXdB;j+=20){
+    float y = j*(1./MINdB);
+    cairo_move_to(c,.5,y1+h*y+.5);
+    cairo_line_to(c,pic_w-.5,y1+h*y+.5);
+  }
+  cairo_stroke(c);
+  for(j=MINdB;j<=MAXdB;j+=20){
+    float y = j*(1./MINdB);
+    cairo_text_extents_t extents;
+    char b[80];
+    sprintf(b,"%ddB",j);
+    cairo_text_extents(c, b, &extents);
+    cairo_move_to(c,.5,y1+h*y+extents.height+2.5);
+
+    cairo_set_source_rgb(c, 0,0,0);
+    cairo_text_path (c, b);
+    cairo_set_line_width(c,3.);
+    cairo_set_line_join(c,CAIRO_LINE_JOIN_ROUND);
+    cairo_stroke(c);
+
+    cairo_set_source_rgb(c,.2,.2,.2);
+    cairo_move_to(c,.5,y1+h*y+extents.height+2.5);
+    cairo_show_text(c, b);
+  }
+  cairo_restore(c);
+}
+
+void frequency_vector_dB(float *vec,int ch,int n,
+                         int area, cairo_t *c,int fillp){
+  float r,g,b;
+  int y1 = area_y(area);
+  int h = area_h(area);
+  int hzoom = AREA_HZOOM[area];
+  int j;
+
+  for(j=0;j<n;j++){
+    float y = vec[j]*(1./MINdB);
+    if(y<0.)y=0.;
+    if(y>1.)y=1.;
+    if(j==0){
+      cairo_move_to(c,j*hzoom+.5,y1+h*y+.5);
+    }else{
+      cairo_line_to(c,j*hzoom+.5,y1+h*y+.5);
+    }
+  }
+  pick_a_color(ch,&r,&g,&b);
+  if(fillp){
+    cairo_line_to(c,pic_w-.5,y1+h-.5);
+    cairo_line_to(c,.5,y1+h-.5);
+    cairo_close_path(c);
+    cairo_set_source_rgba(c,r,g,b,.2);
+    cairo_fill(c);
+    cairo_new_path(c);
+  }else{
+    cairo_set_source_rgb(c,r,g,b);
+    cairo_stroke(c);
+  }
+
+}
+
+void spectrum(pcm_t *pcm,cairo_t *c,off_t pos){
+  float work[BLOCKSIZE];
+  int i,j;
+  cairo_save(c);
+  clip_area(c,SPECTRUM_AREA);
+  clear_area(c,SPECTRUM_AREA);
+  grid_lines(c,SPECTRUM_AREA);
+  cairo_set_line_width(c,.7);
+
+  for(i=0;i<pcm->ch;i++){
+    if(pos+BLOCKSIZE<=pcm->samples)
+      memcpy(work, &pcm->data[i][pos],BLOCKSIZE*sizeof(*work));
+    else{
+      memset(work, 0, sizeof(work));
+      memcpy(work, &pcm->data[i][pos],(pcm->samples-pos)*sizeof(*work));
+    }
+
+    /* window and drft */
+    hanning(work, work,BLOCKSIZE);
+    drft_forward(&fft, work);
+    for(j=0;j<BLOCKSIZE;j++)
+      work[j]*=4./BLOCKSIZE;
+
+    mag_dB(work,work,BLOCKSIZE);
+    frequency_vector_dB(work,i,BLOCKSIZE/2+1,SPECTRUM_AREA,c,1);
+  }
+  cairo_restore(c);
+}
+
+void spread_bark_naive(float *in, float *out, int n, int rate,
+                       float loslope, float width, float hislope){
+
+  static float *spread_integral=NULL;
+  static float *bark=NULL;
+  static int spread_integral_n=0;
+
+  int i,j;
+  memset(out,0,sizeof(*out)*n);
+
+  if(spread_integral_n!=n){
+    float lobarkwidth=toBark(.5/n*rate)-toBark(0);
+    float hibarkwidth=toBark(.5*rate)-toBark(.5*(n-1)/n*rate);
+
+    if(spread_integral)free(spread_integral);
+    if(bark)free(bark);
+    spread_integral = calloc(n,sizeof(*spread_integral));
+    bark = calloc(n,sizeof(*bark));
+    spread_integral_n = n;
+
+    for(i=0;i<n;i++){
+      float total=1.;
+      float att=0.;
+      float centerbark = toBark(.5*i/n*rate);
+      j=i;
+      bark[i]=centerbark;
+
+      /* loop for computation that assumes we've not fallen off bottom
+         edge of the spectrum */
+      while(--j>=0){
+        float b = toBark(.5*j/n*rate);
+        if(b>centerbark-width*.5){
+          /* still inside the F3 point */
+          att = (centerbark-b)/(width*.5)*-3.;
+        }else{
+          /* outside the F3 point */
+          att = ((centerbark-width*.5)-b)*loslope-3.;
+        }
+        total += fromdB(att);
+      }
+
+      /* avoid area contraction at edges of spectrum */
+      while(att>-120.){
+        att += (lobarkwidth*loslope);
+        total += fromdB(att);
+      }
+
+      j=i;
+      /* loop for computation that assumes we've not fallen off top
+         edge of the spectrum */
+      while(++j<n){
+        float b = toBark(.5*j/n*rate);
+        if(b<centerbark+width*.5){
+          /* still inside the F3 point */
+          att = (b-centerbark)/(width*.5)*-3.;
+        }else{
+          /* outside the F3 point */
+          att = (b-(width*.5+centerbark))*hislope-3.;
+        }
+        total += fromdB(att);
+      }
+
+      /* avoid area contraction at edges of spectrum */
+      while(att>-120.){
+        att += (hibarkwidth*hislope);
+        total += fromdB(att);
+      }
+
+      /* numerical integration complete to well beyond limits of reason */
+      spread_integral[i]=total;
+    }
+  }
+
+  for(i=0;i<n;i++){
+    float centerbark = bark[i];
+    float e = fromdB(in[i])*fromdB(in[i])/spread_integral[i];
+    j=i;
+
+    out[j] += e;
+    while(--j>=0){
+      float b = bark[j];
+      float att;
+      if(b>centerbark-width*.5){
+        /* still inside the F3 point */
+        att = (centerbark-b)/(width*.5)*-3.;
+      }else{
+        /* outside the F3 point */
+        att = ((centerbark-width*.5)-b)*loslope-3.;
+      }
+      out[j] += fromdB(att)*e;
+    }
+
+    j=i;
+    while(++j<n){
+      float b = bark[j];
+      float att;
+      if(b<centerbark+width*.5){
+        /* still inside the F3 point */
+        att = (b-centerbark)/(width*.5)*-3.;
+      }else{
+        /* outside the F3 point */
+        att = (b-(width*.5+centerbark))*hislope-3.;
+      }
+      out[j] += fromdB(att)*e;
+    }
+  }
+
+  for(i=0;i<n;i++)
+    out[i] = todB(sqrt(out[i]));
+
+}
+
+void spread_bark(float *in, float *out, int n, int rate,
+                 float loslope, float width, float hislope){
+
+  static float *spread_att=NULL;
+  static float *barkw=NULL;
+  static int spread_n=0;
+  static int *transition_lo=NULL;
+  static int *transition_hi=NULL;
+  static float *transition_lo_att=NULL;
+  static float *transition_hi_att=NULL;
+
+  float work[n];
+
+  int i,j;
+
+  if(spread_n!=n){
+    float b=toBark(0);
+
+    if(spread_att){
+      free(spread_att);
+      free(barkw);
+      free(transition_lo);
+      free(transition_hi);
+      free(transition_lo_att);
+      free(transition_hi_att);
+    }
+    spread_att = calloc(n,sizeof(*spread_att));
+    barkw = calloc(n,sizeof(*barkw));
+    transition_lo = calloc(n,sizeof(*transition_lo));
+    transition_hi = calloc(n,sizeof(*transition_hi));
+    transition_lo_att = calloc(n,sizeof(*transition_lo_att));
+    transition_hi_att = calloc(n,sizeof(*transition_hi_att));
+    spread_n = n;
+
+    /* used to determine rolloff per linear step */
+    for(i=0;i<n;i++){
+      float nb = toBark(.5*(i+1)/n*rate);
+      barkw[i]=nb-b;
+      b=nb;
+    }
+
+    /* compute the energy integral */
+    for(i=0;i<n;i++){
+      float total=1.;
+      float att=0.;
+      float centerbark = toBark(.5*i/n*rate);
+      j=i;
+
+      /* loop for computation that assumes we've not fallen off bottom
+         edge of the spectrum */
+      while(--j>=0 && att>-120.){
+        float b = toBark(.5*j/n*rate);
+        if(b>centerbark-width*.5){
+          /* still inside the F3 point */
+          att = (centerbark-b)/(width*.5)*-3.;
+        }else{
+          /* outside the F3 point */
+          att = ((centerbark-width*.5)-b)*loslope-3.;
+        }
+        total += fromdB(att);
+      }
+
+      /* avoid area contraction at edges of spectrum */
+      while(att>-120.){
+        att += (barkw[0]*loslope);
+        total += fromdB(att);
+      }
+
+      j=i;
+      att=0;
+      /* loop for computation that assumes we've not fallen off top
+         edge of the spectrum */
+      while(++j<n && att>-120.){
+        float b = toBark(.5*j/n*rate);
+        if(b<centerbark+width*.5){
+          /* still inside the F3 point */
+          att = (b-centerbark)/(width*.5)*-3.;
+        }else{
+          /* outside the F3 point */
+          att = (b-(width*.5+centerbark))*hislope-3.;
+        }
+        total += fromdB(att);
+      }
+
+      /* avoid area contraction at edges of spectrum */
+      while(att>-120.){
+        att += (barkw[n-1]*hislope);
+        total += fromdB(att);
+      }
+
+      /* numerical integration complete to well beyond limits of reason */
+      spread_att[i]=1./total;
+    }
+
+    /* precompute transition bins and attenuation at transition
+       point. Do it empirically in a form parallel to actual
+       computation for good roundoff match (and easy
+       implementation) */
+
+    for(i=0;i<n;i++){
+      float att = 1.;
+      float bw=barkw[i];
+      for(j=i-1;j>=0 && bw<width*.5;j--){
+        bw += barkw[j];
+        att *= fromdB(-6.*barkw[j+1]/width);
+      }
+      if(j>=0){
+        transition_lo[i]=j+1;
+        transition_lo_att[i]=att;
+      }else{
+        transition_lo[i]=0;
+        transition_lo_att[i]=0.;
+      }
+
+      att = 1.;
+      bw=0.;
+      for(j=i+1;j<n && bw<width*.5;j++){
+        bw += barkw[j];
+        att *= fromdB(-6.*barkw[j-1]/width);
+      }
+      if(j<n){
+        transition_hi[i]=j-1;
+        transition_hi_att[i]=att;
+      }else{
+        transition_hi[i]=n-1;
+        transition_hi_att[i]=0.;
+      }
+    }
+  }
+
+  /* seed filter transitions */
+  memset(out,0,sizeof(*out)*n);
+  memset(work,0,sizeof(work));
+  for(i=0;i<n;i++){
+    float e0 = fromdB(in[i])*fromdB(in[i])*spread_att[i];
+    out[transition_lo[i]]+=e0*transition_lo_att[i];
+  }
+
+  float dB3=0.;
+  float dBlo=0.;
+  for(i=n-1;i>=0;i--){
+    float e0 = fromdB(in[i])*fromdB(in[i])*spread_att[i];
+    dB3  += e0;
+    dB3  -= out[i];
+    dB3   = (dB3<0?0:dB3);
+    dBlo += out[i];
+    out[i] = dB3+dBlo;
+    dB3  *= fromdB(-6.*barkw[i]/width);
+    dBlo *= fromdB(loslope*barkw[i]);
+  }
+
+  for(i=0;i<n;i++){
+    float e0 = fromdB(in[i])*fromdB(in[i])*spread_att[i];
+    work[transition_hi[i]]+=e0*transition_hi_att[i];
+  }
+
+  dB3=0.;
+  dBlo=0.;
+  for(i=0;i<n;i++){
+    float e0 = fromdB(in[i])*fromdB(in[i])*spread_att[i];
+    out[i] += dB3+dBlo;
+
+    dB3  += e0;
+    dB3  -= work[i];
+    dB3   = (dB3<0?0:dB3);
+    dBlo += work[i];
+    dB3  *= fromdB(-6.*barkw[i]/width);
+    dBlo *= fromdB(hislope*barkw[i]);
+  }
+
+
+
+
+  for(i=0;i<n;i++)
+    out[i]=todB(sqrt(out[i]));
+
+}
+
+void envelope(pcm_t *pcm,cairo_t *c,off_t pos){
+  float work[BLOCKSIZE];
+  float work2[BLOCKSIZE];
+  int i,j;
+  cairo_save(c);
+  clip_area(c,SPECTRUM_AREA);
+  cairo_set_line_width(c,.7);
+
+  for(i=0;i<pcm->ch;i++){
+    if(pos+BLOCKSIZE<=pcm->samples)
+      memcpy(work, &pcm->data[i][pos],BLOCKSIZE*sizeof(*work));
+    else{
+      memset(work, 0, sizeof(work));
+      memcpy(work, &pcm->data[i][pos],(pcm->samples-pos)*sizeof(*work));
+    }
+
+    /* window and drft */
+    hanning(work, work,BLOCKSIZE);
+    drft_forward(&fft, work);
+    for(j=0;j<BLOCKSIZE;j++)
+      work[j]*=4./BLOCKSIZE;
+
+    mag_dB(work,work,BLOCKSIZE);
+    spread_bark(work, work2, BLOCKSIZE/2+1, pcm->rate,
+                -20., 1., -20.);
+    frequency_vector_dB(work2,i+pcm->ch,BLOCKSIZE/2+1,ENVELOPE_AREA,c,0);
+
+  }
+  cairo_restore(c);
+}
+
+int main(int argc, char **argv){
+  pcm_t *pcm;
+  cairo_t *c;
+
+  int areas_active[AREAS]={1,1};
+  off_t areas_pos[AREAS]={0,0};
+
+  pic_w = BLOCKSIZE/2;
+  pic_h = pic_w*9/16;
+
+
+  drft_init(&fft, BLOCKSIZE);
+
+  /* squishyload a file. */
+  pcm = squishyio_load_file(argv[1]);
+
+  /* Make a cairo drawable */
+  cairo_surface_t *pane=cairo_image_surface_create(CAIRO_FORMAT_RGB24,pic_w, pic_h);
+  if(!pane || cairo_surface_status(pane)!=CAIRO_STATUS_SUCCESS){
+    fprintf(stderr,"Could not set up Cairo surface.\n\n");
+    exit(1);
+  }
+
+  c=cairo_create(pane);
+  cairo_set_source_rgb(c,0,0,0);
+  cairo_rectangle(c,0,0,pic_w,pic_h);
+  cairo_fill(c);
+  cairo_set_font_size(c, SPECTRUM_GRIDFONTSIZE);
+
+  printf("YUV4OGG Sy\n");
+  printf("VIDEO W%d H%d F%d:%d Ip A1:1 C444\n",
+         pic_w,pic_h,pcm->rate,BLOCKSIZE/AREA_VZOOM[AREA_FRAMEOUT_TRIGGER]);
+  printf("AUDIO R%d C%d\n",pcm->rate,pcm->ch);
+
+  while(1){
+    int i;
+
+    /* which area is next? */
+    off_t min=0;
+    off_t min_area=-1;
+    for(i=0;i<AREAS;i++)
+      if(areas_active[i] &&
+         (min_area==-1 || min>areas_pos[i])){
+        min=areas_pos[i];
+        min_area=i;
+      }
+
+    if(min>=pcm->samples) break;
+
+    if(min_area == SPECTROGRAM_AREA)
+      spectrogram(pcm,c,areas_pos[min_area]);
+
+    if(min_area == SPECTRUM_AREA)
+      spectrum(pcm,c,areas_pos[min_area]);
+
+    if(min_area == ENVELOPE_AREA)
+      envelope(pcm,c,areas_pos[min_area]);
+
+
+    /* lastly */
+    if(min_area == AREA_FRAMEOUT_TRIGGER){
+      write_yuv(pcm,pane);
+      write_pcm(pcm,pcm->data,areas_pos[min_area],BLOCKSIZE/AREA_VZOOM[min_area]);
+    }
+    areas_pos[min_area]+=BLOCKSIZE/AREA_VZOOM[min_area];
+  }
+  return 0;
+}
+

Copied: trunk/ghost/monty/chirp/scales.h (from rev 17893, trunk/ghost/monty/sinusoid/scales.h)
===================================================================
--- trunk/ghost/monty/chirp/scales.h	                        (rev 0)
+++ trunk/ghost/monty/chirp/scales.h	2011-03-18 06:06:42 UTC (rev 17895)
@@ -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/ghost/monty/chirp/smallft.c (from rev 17893, trunk/ghost/monty/sinusoid/smallft.c)
===================================================================
--- trunk/ghost/monty/chirp/smallft.c	                        (rev 0)
+++ trunk/ghost/monty/chirp/smallft.c	2011-03-18 06:06:42 UTC (rev 17895)
@@ -0,0 +1,523 @@
+/********************************************************************
+ *                                                                  *
+ * 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 1994-2007              *
+ * the Xiph.Org FOundation http://www.xiph.org/                     *
+ *                                                                  *
+ ********************************************************************
+
+ function: unnormalized powers-of-two forward fft transform
+ last mod: $Id: smallft.c 7573 2004-08-16 01:26:52Z conrad $
+
+ ********************************************************************/
+
+/* FFT implementation from OggSquish, minus cosine transforms,
+ * minus all but radix 2/4 case.  In Vorbis we only need this
+ * cut-down version.
+ *
+ * To do more than just power-of-two sized vectors, see the full
+ * version I wrote for NetLib.
+ *
+ * Note that the packing is a little strange; rather than the FFT r/i
+ * packing following R_0, R_n, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1,
+ * it follows R_0, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1, I_n like the
+ * FORTRAN version
+ */
+
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "smallft.h"
+
+static int drfti1(int n, float *wa, int *ifac){
+  static int ntryh[2] = { 4,2 };
+  static float tpi = 6.28318530717958648f;
+  float arg,argh,argld,fi;
+  int ntry=0,i,j=-1;
+  int k1, l1, l2, ib;
+  int ld, ii, ip, is, nq, nr;
+  int ido, ipm, nfm1;
+  int nl=n;
+  int nf=0;
+
+ L101:
+  j++;
+  if (j < 2)
+    ntry=ntryh[j];
+  else // powers of two only in this one, sorry
+    return -1;
+
+ L104:
+  nq=nl/ntry;
+  nr=nl-ntry*nq;
+  if (nr!=0) goto L101;
+
+  nf++;
+  ifac[nf+1]=ntry;
+  nl=nq;
+  if(ntry!=2)goto L107;
+  if(nf==1)goto L107;
+
+  for (i=1;i<nf;i++){
+    ib=nf-i+1;
+    ifac[ib+1]=ifac[ib];
+  }
+  ifac[2] = 2;
+
+ L107:
+  if(nl!=1)goto L104;
+  ifac[0]=n;
+  ifac[1]=nf;
+  argh=tpi/n;
+  is=0;
+  nfm1=nf-1;
+  l1=1;
+
+  if(nfm1==0)return -1;
+
+  for (k1=0;k1<nfm1;k1++){
+    ip=ifac[k1+2];
+    ld=0;
+    l2=l1*ip;
+    ido=n/l2;
+    ipm=ip-1;
+
+    for (j=0;j<ipm;j++){
+      ld+=l1;
+      i=is;
+      argld=(float)ld*argh;
+      fi=0.;
+      for (ii=2;ii<ido;ii+=2){
+	fi+=1.;
+	arg=fi*argld;
+	wa[i++]=cos(arg);
+	wa[i++]=sin(arg);
+      }
+      is+=ido;
+    }
+    l1=l2;
+  }
+  return 0;
+}
+
+static int fdrffti(int n, float *wsave, int *ifac){
+  if (n == 1) return 0;
+  return drfti1(n, wsave, ifac);
+}
+
+static void dradf2(int ido,int l1,float *cc,float *ch,float *wa1){
+  int i,k;
+  float ti2,tr2;
+  int t0,t1,t2,t3,t4,t5,t6;
+
+  t1=0;
+  t0=(t2=l1*ido);
+  t3=ido<<1;
+  for(k=0;k<l1;k++){
+    ch[t1<<1]=cc[t1]+cc[t2];
+    ch[(t1<<1)+t3-1]=cc[t1]-cc[t2];
+    t1+=ido;
+    t2+=ido;
+  }
+    
+  if(ido<2)return;
+  if(ido==2)goto L105;
+
+  t1=0;
+  t2=t0;
+  for(k=0;k<l1;k++){
+    t3=t2;
+    t4=(t1<<1)+(ido<<1);
+    t5=t1;
+    t6=t1+t1;
+    for(i=2;i<ido;i+=2){
+      t3+=2;
+      t4-=2;
+      t5+=2;
+      t6+=2;
+      tr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];
+      ti2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];
+      ch[t6]=cc[t5]+ti2;
+      ch[t4]=ti2-cc[t5];
+      ch[t6-1]=cc[t5-1]+tr2;
+      ch[t4-1]=cc[t5-1]-tr2;
+    }
+    t1+=ido;
+    t2+=ido;
+  }
+
+  if(ido%2==1)return;
+
+ L105:
+  t3=(t2=(t1=ido)-1);
+  t2+=t0;
+  for(k=0;k<l1;k++){
+    ch[t1]=-cc[t2];
+    ch[t1-1]=cc[t3];
+    t1+=ido<<1;
+    t2+=ido;
+    t3+=ido;
+  }
+}
+
+static void dradf4(int ido,int l1,float *cc,float *ch,float *wa1,
+	    float *wa2,float *wa3){
+  static float hsqt2 = .70710678118654752f;
+  int i,k,t0,t1,t2,t3,t4,t5,t6;
+  float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
+  t0=l1*ido;
+  
+  t1=t0;
+  t4=t1<<1;
+  t2=t1+(t1<<1);
+  t3=0;
+
+  for(k=0;k<l1;k++){
+    tr1=cc[t1]+cc[t2];
+    tr2=cc[t3]+cc[t4];
+
+    ch[t5=t3<<2]=tr1+tr2;
+    ch[(ido<<2)+t5-1]=tr2-tr1;
+    ch[(t5+=(ido<<1))-1]=cc[t3]-cc[t4];
+    ch[t5]=cc[t2]-cc[t1];
+
+    t1+=ido;
+    t2+=ido;
+    t3+=ido;
+    t4+=ido;
+  }
+
+  if(ido<2)return;
+  if(ido==2)goto L105;
+
+
+  t1=0;
+  for(k=0;k<l1;k++){
+    t2=t1;
+    t4=t1<<2;
+    t5=(t6=ido<<1)+t4;
+    for(i=2;i<ido;i+=2){
+      t3=(t2+=2);
+      t4+=2;
+      t5-=2;
+
+      t3+=t0;
+      cr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];
+      ci2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];
+      t3+=t0;
+      cr3=wa2[i-2]*cc[t3-1]+wa2[i-1]*cc[t3];
+      ci3=wa2[i-2]*cc[t3]-wa2[i-1]*cc[t3-1];
+      t3+=t0;
+      cr4=wa3[i-2]*cc[t3-1]+wa3[i-1]*cc[t3];
+      ci4=wa3[i-2]*cc[t3]-wa3[i-1]*cc[t3-1];
+
+      tr1=cr2+cr4;
+      tr4=cr4-cr2;
+      ti1=ci2+ci4;
+      ti4=ci2-ci4;
+
+      ti2=cc[t2]+ci3;
+      ti3=cc[t2]-ci3;
+      tr2=cc[t2-1]+cr3;
+      tr3=cc[t2-1]-cr3;
+
+      ch[t4-1]=tr1+tr2;
+      ch[t4]=ti1+ti2;
+
+      ch[t5-1]=tr3-ti4;
+      ch[t5]=tr4-ti3;
+
+      ch[t4+t6-1]=ti4+tr3;
+      ch[t4+t6]=tr4+ti3;
+
+      ch[t5+t6-1]=tr2-tr1;
+      ch[t5+t6]=ti1-ti2;
+    }
+    t1+=ido;
+  }
+  if(ido&1)return;
+
+ L105:
+  
+  t2=(t1=t0+ido-1)+(t0<<1);
+  t3=ido<<2;
+  t4=ido;
+  t5=ido<<1;
+  t6=ido;
+
+  for(k=0;k<l1;k++){
+    ti1=-hsqt2*(cc[t1]+cc[t2]);
+    tr1=hsqt2*(cc[t1]-cc[t2]);
+
+    ch[t4-1]=tr1+cc[t6-1];
+    ch[t4+t5-1]=cc[t6-1]-tr1;
+
+    ch[t4]=ti1-cc[t1+t0];
+    ch[t4+t5]=ti1+cc[t1+t0];
+
+    t1+=ido;
+    t2+=ido;
+    t4+=t3;
+    t6+=ido;
+  }
+}
+
+static void drftf1(int n,float *c,float *ch,float *wa,int *ifac){
+  int i,k1,l1,l2;
+  int na,kh,nf;
+  int ip,iw,ido,idl1,ix2,ix3;
+
+  nf=ifac[1];
+  na=1;
+  l2=n;
+  iw=n;
+
+  for(k1=0;k1<nf;k1++){
+    kh=nf-k1;
+    ip=ifac[kh+1];
+    l1=l2/ip;
+    ido=n/l2;
+    idl1=ido*l1;
+    iw-=(ip-1)*ido;
+    na=1-na;
+
+    if(ip!=4)goto L102;
+
+    ix2=iw+ido;
+    ix3=ix2+ido;
+    if(na!=0)
+      dradf4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);
+    else
+      dradf4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);
+    goto L110;
+
+ L102:
+    if(na!=0)goto L103;
+
+    dradf2(ido,l1,c,ch,wa+iw-1);
+    goto L110;
+
+  L103:
+    dradf2(ido,l1,ch,c,wa+iw-1);
+
+  L110:
+    l2=l1;
+  }
+
+  if(na==1)return;
+
+  for(i=0;i<n;i++)c[i]=ch[i];
+}
+
+static void dradb2(int ido,int l1,float *cc,float *ch,float *wa1){
+  int i,k,t0,t1,t2,t3,t4,t5,t6;
+  float ti2,tr2;
+
+  t0=l1*ido;
+
+  t1=0;
+  t2=0;
+  t3=(ido<<1)-1;
+  for(k=0;k<l1;k++){
+    ch[t1]=cc[t2]+cc[t3+t2];
+    ch[t1+t0]=cc[t2]-cc[t3+t2];
+    t2=(t1+=ido)<<1;
+  }
+
+  if(ido<2)return;
+  if(ido==2)goto L105;
+
+  t1=0;
+  t2=0;
+  for(k=0;k<l1;k++){
+    t3=t1;
+    t5=(t4=t2)+(ido<<1);
+    t6=t0+t1;
+    for(i=2;i<ido;i+=2){
+      t3+=2;
+      t4+=2;
+      t5-=2;
+      t6+=2;
+      ch[t3-1]=cc[t4-1]+cc[t5-1];
+      tr2=cc[t4-1]-cc[t5-1];
+      ch[t3]=cc[t4]-cc[t5];
+      ti2=cc[t4]+cc[t5];
+      ch[t6-1]=wa1[i-2]*tr2-wa1[i-1]*ti2;
+      ch[t6]=wa1[i-2]*ti2+wa1[i-1]*tr2;
+    }
+    t2=(t1+=ido)<<1;
+  }
+
+  if(ido%2==1)return;
+
+L105:
+  t1=ido-1;
+  t2=ido-1;
+  for(k=0;k<l1;k++){
+    ch[t1]=cc[t2]+cc[t2];
+    ch[t1+t0]=-(cc[t2+1]+cc[t2+1]);
+    t1+=ido;
+    t2+=ido<<1;
+  }
+}
+
+static void dradb4(int ido,int l1,float *cc,float *ch,float *wa1,
+                          float *wa2,float *wa3){
+  static float sqrt2=1.414213562373095f;
+  int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8;
+  float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
+  t0=l1*ido;
+
+  t1=0;
+  t2=ido<<2;
+  t3=0;
+  t6=ido<<1;
+  for(k=0;k<l1;k++){
+    t4=t3+t6;
+    t5=t1;
+    tr3=cc[t4-1]+cc[t4-1];
+    tr4=cc[t4]+cc[t4];
+    tr1=cc[t3]-cc[(t4+=t6)-1];
+    tr2=cc[t3]+cc[t4-1];
+    ch[t5]=tr2+tr3;
+    ch[t5+=t0]=tr1-tr4;
+    ch[t5+=t0]=tr2-tr3;
+    ch[t5+=t0]=tr1+tr4;
+    t1+=ido;
+    t3+=t2;
+  }
+
+  if(ido<2)return;
+  if(ido==2)goto L105;
+
+  t1=0;
+  for(k=0;k<l1;k++){
+    t5=(t4=(t3=(t2=t1<<2)+t6))+t6;
+    t7=t1;
+    for(i=2;i<ido;i+=2){
+      t2+=2;
+      t3+=2;
+      t4-=2;
+      t5-=2;
+      t7+=2;
+      ti1=cc[t2]+cc[t5];
+      ti2=cc[t2]-cc[t5];
+      ti3=cc[t3]-cc[t4];
+      tr4=cc[t3]+cc[t4];
+      tr1=cc[t2-1]-cc[t5-1];
+      tr2=cc[t2-1]+cc[t5-1];
+      ti4=cc[t3-1]-cc[t4-1];
+      tr3=cc[t3-1]+cc[t4-1];
+      ch[t7-1]=tr2+tr3;
+      cr3=tr2-tr3;
+      ch[t7]=ti2+ti3;
+      ci3=ti2-ti3;
+      cr2=tr1-tr4;
+      cr4=tr1+tr4;
+      ci2=ti1+ti4;
+      ci4=ti1-ti4;
+
+      ch[(t8=t7+t0)-1]=wa1[i-2]*cr2-wa1[i-1]*ci2;
+      ch[t8]=wa1[i-2]*ci2+wa1[i-1]*cr2;
+      ch[(t8+=t0)-1]=wa2[i-2]*cr3-wa2[i-1]*ci3;
+      ch[t8]=wa2[i-2]*ci3+wa2[i-1]*cr3;
+      ch[(t8+=t0)-1]=wa3[i-2]*cr4-wa3[i-1]*ci4;
+      ch[t8]=wa3[i-2]*ci4+wa3[i-1]*cr4;
+    }
+    t1+=ido;
+  }
+
+  if(ido%2 == 1)return;
+
+ L105:
+
+  t1=ido;
+  t2=ido<<2;
+  t3=ido-1;
+  t4=ido+(ido<<1);
+  for(k=0;k<l1;k++){
+    t5=t3;
+    ti1=cc[t1]+cc[t4];
+    ti2=cc[t4]-cc[t1];
+    tr1=cc[t1-1]-cc[t4-1];
+    tr2=cc[t1-1]+cc[t4-1];
+    ch[t5]=tr2+tr2;
+    ch[t5+=t0]=sqrt2*(tr1-ti1);
+    ch[t5+=t0]=ti2+ti2;
+    ch[t5+=t0]=-sqrt2*(tr1+ti1);
+
+    t3+=ido;
+    t1+=t2;
+    t4+=t2;
+  }
+}
+
+static void drftb1(int n, float *c, float *ch, float *wa, int *ifac){
+  int i,k1,l1,l2;
+  int na;
+  int nf,ip,iw,ix2,ix3,ido,idl1;
+
+  nf=ifac[1];
+  na=0;
+  l1=1;
+  iw=1;
+
+  for(k1=0;k1<nf;k1++){
+    ip=ifac[k1 + 2];
+    l2=ip*l1;
+    ido=n/l2;
+    idl1=ido*l1;
+
+    if(ip!=4){
+      if(na!=0)
+        dradb2(ido,l1,ch,c,wa+iw-1);
+      else
+        dradb2(ido,l1,c,ch,wa+iw-1);
+    } else {
+      ix2=iw+ido;
+      ix3=ix2+ido;
+      if(na!=0)
+        dradb4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);
+      else
+        dradb4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);
+    }
+    na=1-na;
+    l1=l2;
+    iw+=(ip-1)*ido;
+  }
+
+  if(na==0)return;
+
+  for(i=0;i<n;i++)c[i]=ch[i];
+}
+
+void drft_forward(drft_lookup *l,float *data){
+  float work[l->n];
+  if(l->n==1)return;
+  drftf1(l->n,data,work,l->trigcache,l->splitcache);
+}
+
+void drft_backward(drft_lookup *l,float *data){
+  float work[l->n];
+  if (l->n==1)return;
+  drftb1(l->n,data,work,l->trigcache,l->splitcache);
+}
+
+int drft_init(drft_lookup *l,int n){
+  l->n=n;
+  l->trigcache=calloc(2*n,sizeof(*l->trigcache));
+  l->splitcache=calloc(32,sizeof(*l->splitcache));
+  return fdrffti(n, l->trigcache, l->splitcache);
+}
+
+void drft_clear(drft_lookup *l){
+  if(l){
+    if(l->trigcache)free(l->trigcache);
+    if(l->splitcache)free(l->splitcache);
+    memset(l,0,sizeof(*l));
+  }
+}

Copied: trunk/ghost/monty/chirp/smallft.h (from rev 17893, trunk/ghost/monty/sinusoid/smallft.h)
===================================================================
--- trunk/ghost/monty/chirp/smallft.h	                        (rev 0)
+++ trunk/ghost/monty/chirp/smallft.h	2011-03-18 06:06:42 UTC (rev 17895)
@@ -0,0 +1,45 @@
+/********************************************************************
+ *                                                                  *
+ * 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 1994-2007              *
+ * the Xiph.Org FOundation http://www.xiph.org/                     *
+ *                                                                  *
+ ********************************************************************
+
+ function: unnormalized powers-of-two forward fft transform
+ last mod: $Id: smallft.c 7573 2004-08-16 01:26:52Z conrad $
+
+ ********************************************************************/
+
+/* FFT implementation from OggSquish, minus cosine transforms,
+ * minus all but radix 2/4 case.  In Vorbis we only need this
+ * cut-down version.
+ *
+ * To do more than just power-of-two sized vectors, see the full
+ * version I wrote for NetLib.
+ *
+ * Note that the packing is a little strange; rather than the FFT r/i
+ * packing following R_0, R_n, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1,
+ * it follows R_0, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1, I_n like the
+ * FORTRAN version
+ */
+
+#ifndef _V_SMFT_H_
+#define _V_SMFT_H_
+
+typedef struct {
+  int n;
+  float *trigcache;
+  int *splitcache;
+} drft_lookup;
+
+extern void drft_forward(drft_lookup *l,float *data);
+extern void drft_backward(drft_lookup *l,float *data);
+extern int drft_init(drft_lookup *l,int n);
+extern void drft_clear(drft_lookup *l);
+
+#endif

Copied: trunk/ghost/monty/chirp/window.c (from rev 17893, trunk/ghost/monty/sinusoid/window.c)
===================================================================
--- trunk/ghost/monty/chirp/window.c	                        (rev 0)
+++ trunk/ghost/monty/chirp/window.c	2011-03-18 06:06:42 UTC (rev 17895)
@@ -0,0 +1,38 @@
+#define _GNU_SOURCE
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <cairo/cairo.h>
+#include <squishyio/squishyio.h>
+#include "sinusoids.h"
+#include "smallft.h"
+#include "scales.h"
+
+/* A few windows */
+
+#define A0 .35875f
+#define A1 .48829f
+#define A2 .14128f
+#define A3 .01168f
+
+void blackmann_harris(float *out, float *in, 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);
+    out[i] = in[i]*w;
+  }
+}
+
+static void hanning(float *out, float *in, int n){
+  int i;
+  float scale = 2*M_PI/n;
+
+  for(i=0;i<n;i++){
+    float i5 = i+.5;
+    out[i] = in[i]*(.5-.5*cos(scale*i5));
+  }
+}



More information about the commits mailing list