[xiph-commits] r17893 - trunk/ghost/monty/sinusoid
xiphmont at svn.xiph.org
xiphmont at svn.xiph.org
Thu Mar 17 23:03:20 PDT 2011
Author: xiphmont
Date: 2011-03-17 23:03:19 -0700 (Thu, 17 Mar 2011)
New Revision: 17893
Added:
trunk/ghost/monty/sinusoid/bessel.c
trunk/ghost/monty/sinusoid/chirp.c
trunk/ghost/monty/sinusoid/chirp.h
trunk/ghost/monty/sinusoid/harness.c
trunk/ghost/monty/sinusoid/window.c
Removed:
trunk/ghost/monty/sinusoid/sinusoids.h
Modified:
trunk/ghost/monty/sinusoid/smallft.c
trunk/ghost/monty/sinusoid/smallft.h
Log:
Clean up a little, get the most recent chirp estimation code in SVN
Added: trunk/ghost/monty/sinusoid/bessel.c
===================================================================
--- trunk/ghost/monty/sinusoid/bessel.c (rev 0)
+++ trunk/ghost/monty/sinusoid/bessel.c 2011-03-18 06:03:19 UTC (rev 17893)
@@ -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);
+}
+
Added: trunk/ghost/monty/sinusoid/chirp.c
===================================================================
--- trunk/ghost/monty/sinusoid/chirp.c (rev 0)
+++ trunk/ghost/monty/sinusoid/chirp.c 2011-03-18 06:03:19 UTC (rev 17893)
@@ -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/sinusoid/chirp.h (from rev 13511, trunk/ghost/monty/sinusoid/sinusoids.h)
===================================================================
--- trunk/ghost/monty/sinusoid/chirp.h (rev 0)
+++ trunk/ghost/monty/sinusoid/chirp.h 2011-03-18 06:03:19 UTC (rev 17893)
@@ -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);
+
Added: trunk/ghost/monty/sinusoid/harness.c
===================================================================
--- trunk/ghost/monty/sinusoid/harness.c (rev 0)
+++ trunk/ghost/monty/sinusoid/harness.c 2011-03-18 06:03:19 UTC (rev 17893)
@@ -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;
+}
+
Deleted: trunk/ghost/monty/sinusoid/sinusoids.h
===================================================================
--- trunk/ghost/monty/sinusoid/sinusoids.h 2011-03-17 17:01:48 UTC (rev 17892)
+++ trunk/ghost/monty/sinusoid/sinusoids.h 2011-03-18 06:03:19 UTC (rev 17893)
@@ -1,37 +0,0 @@
-/********************************************************************
- * *
- * 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: research-grade sinusoidal extraction code
- last mod: $Id$
-
- ********************************************************************/
-
-extern void window_weight(float *logf, float *out, int n, float flatbias,
- int lowindow, int hiwindow, int min, int rate);
-
-extern int extract_sinusoids(float *x,
- float *A, float *W, float *P,
- float *dA, float *dW,
- float *y, int N, int len, int count);
-
-extern void extract_modulated_sinusoids(float *x, float *w, float *Aout, float *Wout, float *Pout,
- float *delAout, float *delWout, float *ddAout, float *y, int N, int len);
-
-extern void extract_modulated_sinusoidsB(float *x, float *w,
- float *Aout, float *Wout, float *Pout,
- float *dAout, float *dWout, float *ddAout,
- float *y, int N, int len);
-
-extern void extract_modulated_sinusoids_nonlinear(float *x, float *window,
- float *A, float *W, float *P,
- float *dA, float *dW,
- float *y, int N, int len);
Modified: trunk/ghost/monty/sinusoid/smallft.c
===================================================================
--- trunk/ghost/monty/sinusoid/smallft.c 2011-03-17 17:01:48 UTC (rev 17892)
+++ trunk/ghost/monty/sinusoid/smallft.c 2011-03-18 06:03:19 UTC (rev 17893)
@@ -313,12 +313,200 @@
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));
Modified: trunk/ghost/monty/sinusoid/smallft.h
===================================================================
--- trunk/ghost/monty/sinusoid/smallft.h 2011-03-17 17:01:48 UTC (rev 17892)
+++ trunk/ghost/monty/sinusoid/smallft.h 2011-03-18 06:03:19 UTC (rev 17893)
@@ -38,6 +38,7 @@
} 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);
Added: trunk/ghost/monty/sinusoid/window.c
===================================================================
--- trunk/ghost/monty/sinusoid/window.c (rev 0)
+++ trunk/ghost/monty/sinusoid/window.c 2011-03-18 06:03:19 UTC (rev 17893)
@@ -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