[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