[xiph-commits] r13235 - in trunk/ghost/monty: . lift sinusoid
xiphmont at svn.xiph.org
xiphmont at svn.xiph.org
Fri Jul 6 07:21:35 PDT 2007
Author: xiphmont
Date: 2007-07-06 07:21:34 -0700 (Fri, 06 Jul 2007)
New Revision: 13235
Added:
trunk/ghost/monty/sinusoid/
trunk/ghost/monty/sinusoid/scales.h
trunk/ghost/monty/sinusoid/sinusoids.c
trunk/ghost/monty/sinusoid/sinusoids.h
trunk/ghost/monty/sinusoid/smallft.c
trunk/ghost/monty/sinusoid/smallft.h
trunk/ghost/monty/sinusoid/test_sinusoids.c
Modified:
trunk/ghost/monty/lift/liftsearch.c
trunk/ghost/monty/lift/smallft.c
Log:
"Thinking aloud" code, nothing more.
Modified: trunk/ghost/monty/lift/liftsearch.c
===================================================================
--- trunk/ghost/monty/lift/liftsearch.c 2007-07-06 13:44:36 UTC (rev 13234)
+++ trunk/ghost/monty/lift/liftsearch.c 2007-07-06 14:21:34 UTC (rev 13235)
@@ -6,13 +6,13 @@
#include "smallft.h"
int fixed_bits = 12;
+int order = 14;
int num_coefficients = 14;
-int p_delay = 2;
-int q_delay = 3;
+int p_delay, q_delay, r_delay, s_delay;
double TB = .1;
-#define MAX_PZ 32
-#define FS 256
+#define MAX_PZ 64
+#define FS 512
drft_lookup fft;
typedef struct local_minimum{
@@ -48,31 +48,32 @@
for(j=0;j<num_coefficients;j++) mc[j] = c[j];
minimum_list_size++;
+ fprintf(stderr,"\rfinal cost: %f \n",cost);
+
return minimum_list_size-1;
}
-double fixed_to_double(int f){
- return f*pow(2.,(double)-fixed_bits);
-}
-
typedef struct {
int n;
- double a[MAX_PZ]; // normalized filter; a[0] is actually a_1
- double b[MAX_PZ+1];
+ double scale[MAX_PZ];
+ double iir[MAX_PZ];
+ double fir[MAX_PZ+1];
double pre[MAX_PZ*2];
} iir_t;
+
+//Direct form
double Az2(iir_t *f, double v){
double iir = v;
int i;
int n = f->n;
for(i=0;i<n;i++)
- iir -= f->pre[i*2+1] * f->a[i];
+ iir -= f->pre[i*2+1] * f->iir[i];
- double fir = f->b[0] * iir;
+ double fir = f->fir[0] * iir;
for(i=0;i<n;i++)
- fir += f->pre[i*2+1] * f->b[i+1];
+ fir += f->pre[i*2+1] * f->fir[i+1];
for(i = n*2-1; i>0; i--)
f->pre[i] = f->pre[i-1];
@@ -81,6 +82,22 @@
return fir;
}
+// Cascade form
+double Cz2(iir_t *f, double v){
+ int i;
+ int n = f->n;
+
+ v*=f->scale[n-1];
+ for(i=0;i<n;i++){
+ double iir = v - f->pre[i*2] * f->iir[i];
+ v = iir + f->pre[i*2] * f->fir[i];
+
+ f->pre[i*2] = f->pre[i*2+1];
+ f->pre[i*2+1] = iir;
+ }
+ return v;
+}
+
void mag_e(double *d){
int i;
@@ -93,30 +110,57 @@
}
-int lift_run(int *c, double *lp, double *hp, int dump){
- iir_t P;
- iir_t Q;
+void lift_run(double *c, double *lp, double *hp, int dump){
+ iir_t P,Q,R,S;
int i,j;
- int offset = (abs(p_delay) + abs(q_delay))*2+3;
- int ret = 0;
+ int offset = (abs(p_delay) + abs(q_delay))*2+1;
memset(&P,0,sizeof(P));
memset(&Q,0,sizeof(Q));
+ //memset(&R,0,sizeof(R));
+ //memset(&S,0,sizeof(S));
+
lp[offset] = 1.;
hp[offset+1] = 1.;
+ P.n = Q.n = order;
- P.n = Q.n = ((num_coefficients/2) -1)/2;
+#if 1
- for(i=0,j=0; i<P.n; i++, j++)
- P.a[i] = fixed_to_double(c[j]);
- for(i=0; i<P.n+1; i++, j++)
- P.b[i] = fixed_to_double(c[j]);
+ j=0;
+ for(i=0; i<P.n; i++)
+ Q.iir[i] = P.iir[i] = c[j++];
+ for(i=0; i<P.n; i++)
+ Q.fir[i] = P.fir[i] = c[--j];
+ Q.fir[i] = P.fir[i] = 1.;
- for(i=0; i<Q.n; i++, j++)
- Q.a[i] = fixed_to_double(c[j]);
- for(i=0; i<Q.n+1; i++, j++)
- Q.b[i] = fixed_to_double(c[j]);
+#else
+ j=0;
+ for(i=0; i<P.n; i++)
+ P.iir[i] = c[j++];
+ for(i=0; i<P.n+1; i++)
+ P.fir[i] = c[j++];
+ for(i=0; i<Q.n; i++)
+ Q.iir[i] = c[j++];
+ for(i=0; i<Q.n+1; i++)
+ Q.fir[i] = c[j++];
+ //for(i=0; i<R.n; i++)
+ // R.iir[i] = c[j++];
+ //for(i=0; i<R.n+1; i++)
+ // R.fir[i] = c[j++];
+ //for(i=0; i<S.n; i++)
+ // S.iir[i] = c[j++];
+ //for(i=0; i<S.n+1; i++)
+ // S.fir[i] = c[j++];
+#endif
+
+ /*i=0; j=0;
+ if(r_delay>0) i += r_delay*2;
+ if(r_delay<0) j += -r_delay*2;
+
+ for(; i<FS && j<FS; i++, j++)
+ hp[j] -= Az2(&R, lp[i]);*/
+
i=0; j=0;
if(p_delay>0) i += p_delay*2;
if(p_delay<0) j += -p_delay*2;
@@ -131,63 +175,59 @@
for(; i<FS && j<FS; i++, j++)
lp[j] += .5*Az2(&Q, hp[i]);
- // check for "stable enough"
- for(i=0;i<16;i++)
- if(todB(lp[FS-i-1])>-120 ||
- todB(hp[FS-i-1])>-120){
- ret=1;
- if(dump)
- fprintf(stderr,"warning, insufficiently stable filter!\n");
- break;
+
+ /*i=0; j=0;
+ if(s_delay>0) i += s_delay*2;
+ if(s_delay<0) j += -s_delay*2;
+
+ for(; i<FS && j<FS; i++, j++)
+ lp[j] += .5*Az2(&S, hp[i]);*/
+
+ if(dump){
+ for(i=0; i<64; i++){
+ if(todB(hp[FS-i-1])>-100 ||
+ todB(lp[FS-i-1])>-100){
+ fprintf(stderr,"Warning: Insufficiently stable filter!\n");
+ break;
+ }
+ if(todB(hp[FS-i-1])>-120 ||
+ todB(lp[FS-i-1])>-120){
+ fprintf(stderr,"Warning: Gibbs ringing overflow!\n");
+ break;
+ }
}
+ }
+
// energy
mag_e(hp);
mag_e(lp);
if(dump){
- fprintf(stdout,"Filter result: \n");
- fprintf(stdout," P: (");
- for(i=0;i<P.n+1;i++)
- fprintf(stdout,"%s%f",(i==0?"":", "),P.b[i]);
- fprintf(stdout,") / (1");
- for(i=0;i<P.n;i++)
- fprintf(stdout,", %f",P.a[i]);
- fprintf(stdout,")\n");
-
- fprintf(stdout," Q: (");
- for(i=0;i<Q.n+1;i++)
- fprintf(stdout,"%s%f",(i==0?"":", "),Q.b[i]);
- fprintf(stdout,") / (1");
+ fprintf(stderr," a: ( ");
for(i=0;i<Q.n;i++)
- fprintf(stdout,", %f",Q.a[i]);
- fprintf(stdout,")\n");
+ fprintf(stderr,"%f ",Q.iir[i]);
+ fprintf(stderr,")\n");
+
- fprintf(stdout," response: {\n");
for(i=0;i<=FS/2;i++)
fprintf(stdout,"%f %f\n",(double)i/FS, todB(lp[i])/2);
fprintf(stdout,"\n");
for(i=0;i<=FS/2;i++)
fprintf(stdout,"%f %f\n",(double)i/FS, todB(hp[i])/2);
- fprintf(stdout,"}\n\n");
+ fprintf(stdout,"\n");
+
}
}
-double lift_cost(int *c){
+double lift_cost(double *c){
double *hp = calloc (FS,sizeof(*hp));
double *lp = calloc (FS,sizeof(*lp));
int i;
- iir_t P,Q;
-
- memset(&P,0,sizeof(P));
- memset(&Q,0,sizeof(Q));
-
- P.n = Q.n = ((num_coefficients/2) -1)/2;
-
lift_run(c, lp, hp, 0);
// cost
@@ -198,9 +238,17 @@
// energy in stopband
for(i=TBhi; i<=FS/2; i++)
cost += lp[i];
- for(i=0; i<=TBlo; i++)
+ for(i=0; i<TBlo; i++)
cost += hp[i]*.25;
+ // overshoot energy
+ /*for(i=0; i<=FS/2; i++){
+ if(lp[i]>2.)
+ cost += lp[i]-2.;
+ if(hp[i]>8.)
+ cost += (hp[i]*.25-2.);
+ }*/
+
// deviation from 0dB in passband
//for(i=0; i<=TBlo; i++)
// cost += fabs(lp[i]-1.)/4096.;
@@ -218,180 +266,12 @@
free(hp);
free(lp);
- return todB(sqrt(cost/(FS-TBhi+TBlo)));
+ return cost*cost;//(FS/2-TBhi + TBlo);
}
-void drft_sub(double *p, double *n, double *out){
- int i;
- for(i=0;i<FS;i++)
- out[i] = p[i] - n[i];
-}
-
-void drft_add(double *p, double *n, double *out){
- int i;
- for(i=0;i<FS;i++)
- out[i] = p[i] + n[i];
-}
-
-void drft_ediv(double *n, double *d, double *out){
- int i;
- for(i=0;i<=FS/2;i++)
- out[i] = n[i] / d[i];
-}
-
-void drft_mul(double *a, double *b, double *out){
- int i;
- out[0] = a[0] * b[0];
- out[FS-1] = a[FS-1] * b[FS-1];
-
- for(i=2;i<FS-1;i+=2){
- double ra = a[i-1];
- double rb = b[i-1];
- double ia = a[i];
- double ib = b[i];
-
- out[i-1] = ra*rb - ia*ib;
- out[i] = ra*ib + ia*rb;
- }
-}
-
-void drft_energy(double *d){
- int i;
-
- d[0] = d[0]*d[0];
- for(i=2;i<FS;i+=2)
- d[i>>1] =d[i-1]*d[i-1] + d[i]*d[i];
- d[FS/2] = d[FS-1]*d[FS-1];
-}
-
-
-// Directly compute response magnitude of upper/lower branches from
-// filter coefficients. Much more complicated, but also doesn't care about
-// relative filter stability (which may be desirable or undesirable
-// depending) and can use much shorter eval vectors.
-double lift_cost2(int *c){
- double *pn = calloc (FS,sizeof(*pn));
- double *pd = calloc (FS,sizeof(*pd));
- double *qn = calloc (FS,sizeof(*qn));
- double *qd = calloc (FS,sizeof(*qd));
- double *zM = calloc (FS,sizeof(*zM));
- double *one = calloc (FS,sizeof(*one));
- int i,j,k;
-
- int n = ((num_coefficients/2) -1)/2;
-
- /*P.n = Q.n = n;
- for(i=0,j=0; i<P.n; i++, j++)
- P.a[i] = fixed_to_double(c[j]);
- for(i=0; i<P.n+1; i++, j++)
- P.b[i] = fixed_to_double(c[j]);
-
- for(i=0; i<Q.n; i++, j++)
- Q.a[i] = fixed_to_double(c[j]);
- for(i=0; i<Q.n+1; i++, j++)
- Q.b[i] = fixed_to_double(c[j]);
- */
-
- /* set up response computation of highpass section */
- /* D(z) */
- i=0; j=0;
- pd[j] += 1; j+=2;
- for(i=0;i<n;i++,j+=2)
- pd[j] += fixed_to_double(c[i]); //P.a
-
- /* N(z)z^M where M>=0 */
- i=0; j=0;
- if(p_delay<0) j -= p_delay*2;
- for(i=0;i<n+1;i++,j+=2)
- pn[j] += fixed_to_double(c[i+n]); // P.b
-
- /* z^M (shift case for M<0) */
- if(p_delay>0)
- zM[1+p_delay*2] = 1.;
- else
- zM[1] = 1.;
-
- /* 1 */
- one[0] = 1.;
-
- drft_forward(&fft,pn);
- drft_forward(&fft,pd);
- drft_forward(&fft,zM);
- drft_forward(&fft,one);
-
- drft_mul(zM,pd,zM);
-
- drft_mul(pd,one,pd);
- drft_mul(pn,one,pn);
- drft_sub(zM,pn,pn);
-
- /* now the lowpass */
-
- /* D(z) */
- i=0; j=0;
- qd[j] += 1; j+=2;
- for(i=0;i<n;i++,j+=2)
- qd[j] += fixed_to_double(c[i+2*n+1]); //Q.a
-
- /* N(z)z^M where M>=0 */
- i=0; j=0;
- if(q_delay+p_delay<0) j -= (q_delay+p_delay)*2;
- for(i=0;i<n+1;i++,j+=2)
- qn[j] += .5*fixed_to_double(c[i+3*n+1]); // Q.b
-
- /* z^M (shift case for M<0) */
- memset(zM,0,FS*sizeof(*zM));
- if(p_delay+q_delay>0)
- zM[p_delay*2+q_delay*2] = 1.;
- else
- zM[0] = 1.;
-
- drft_forward(&fft,qn);
- drft_forward(&fft,qd);
- drft_forward(&fft,zM);
-
- drft_mul(pn,qn,qn);
- drft_mul(pd,qd,qd);
-
- drft_mul(zM,qd,zM);
-
- drft_mul(qd,one,qd);
- drft_mul(qn,one,qn);
- drft_add(zM,qn,qn);
-
-
-
- drft_energy(pn);
- drft_energy(pd);
- drft_energy(qn);
- drft_energy(qd);
- drft_ediv(pn,pd,pn);
- drft_ediv(qn,qd,qn);
-
- fprintf(stdout," test2 hp response: {\n");
- for(i=0;i<=FS/2;i++)
- fprintf(stdout,"%f %f\n",(double)i/FS, todB(pn[i])/2);
- fprintf(stdout,"}\n\n");
-
- fprintf(stdout," test2 lp response: {\n");
- for(i=0;i<=FS/2;i++)
- fprintf(stdout,"%f %f\n",(double)i/FS, todB(qn[i])/2);
- fprintf(stdout,"}\n\n");
-
- free(pn);
- free(pd);
- free(qn);
- free(qd);
- free(zM);
- free(one);
-
-}
-
-
-void lift_dump(int *c){
+void lift_dump(double *c){
double *hp = calloc (FS,sizeof(*hp));
double *lp = calloc (FS,sizeof(*lp));
- int i;
lift_run(c, lp, hp, 1);
@@ -399,28 +279,28 @@
free(lp);
}
-int walk_to_minimum(int *c){
+int walk_to_minimum_A(double *c, double ep){
int last_change = num_coefficients-1;
int cur_dim = 0;
double cost = lift_cost(c);
while(1){
- c[cur_dim]++;
+ double val = c[cur_dim];
+ c[cur_dim] = val+ep;
double up_cost = lift_cost(c);
- c[cur_dim]-=2;
+ c[cur_dim] = val-ep;
double down_cost = lift_cost(c);
- c[cur_dim]++;
+ c[cur_dim] = val;
if(up_cost<cost && up_cost<down_cost){
- c[cur_dim]++;
+ c[cur_dim]+=ep;
last_change = cur_dim;
cost = up_cost;
}else if (down_cost<cost){
- c[cur_dim]--;
+ c[cur_dim]-=ep;
last_change = cur_dim;
cost = down_cost;
-
}else{
if(cur_dim == last_change)break;
}
@@ -429,11 +309,108 @@
if(cur_dim >= num_coefficients)
cur_dim = 0;
- fprintf(stderr,"\rwalking... current cost: %g ",cost);
+ fprintf(stderr,"\rwalking A... current cost: %g ",todB(cost)/2);
}
- fprintf(stderr,"\rfinal cost: %f \n",cost);
+ fprintf(stderr,"\n");
+ return 0;
+}
- log_minimum(c,cost);
+int walk_to_minimum_CSD(double *c, double ep){
+ double prev[num_coefficients];
+ double r[num_coefficients];
+ double rr[num_coefficients];
+ int have_prev = 0;
+ int flag = 1;
+ double ep2 = 1./(1<<24);
+ double mul=1.;
+ while(flag){
+ int i;
+ flag = 0;
+
+ // compute gradient
+ double d[num_coefficients];
+ double mag = 0;
+ for(i=0;i<num_coefficients;i++){
+ double val = c[i];
+ c[i] = val + ep2;
+ d[i] = lift_cost(c);
+ c[i] = val - ep2;
+ d[i] -= lift_cost(c);
+ c[i] = val;
+ d[i] *= .5;
+ mag += d[i]*d[i];
+ }
+
+ if(mag == 0) return 0;
+ mag = sqrt(mag);
+
+ // normalize gradient to |ep|
+ for(i=0;i<num_coefficients;i++)
+ d[i] *= ep/mag;
+
+ // walk this line
+ memcpy(r,c,sizeof(r));
+ double cost = lift_cost(r);
+ while(1){
+ for(i=0;i<num_coefficients;i++)
+ r[i] -= d[i]*mul;
+ double test_cost = lift_cost(r);
+ if(test_cost < cost){
+ cost = test_cost;
+ flag = 1;
+ mul *= 1.5;
+ if(mul>32)mul=32;
+ }else{
+ for(i=0;i<num_coefficients;i++)
+ r[i] += d[i]*mul;
+ if(mul<1)break;
+ mul *=.5;
+ }
+ }
+
+ /*if(have_prev){
+ mag = 0;
+ for(i=0;i<num_coefficients;i++){
+ d[i] = prev[i] - (c[i]+r[i])*.5;
+ mag += d[i]*d[i];
+ }
+ for(i=0;i<num_coefficients;i++)
+ d[i] /= mag;
+
+ // walk conditioned line
+ memcpy(rr,prev,sizeof(rr));
+ double r_cost = lift_cost(rr);
+ while(1){
+ for(i=0;i<num_coefficients;i++)
+ rr[i] -= d[i]*mul;
+ double test_cost = lift_cost(rr);
+ if(test_cost < r_cost){
+ r_cost = test_cost;
+ mul *= 2;
+ if(mul>32)mul=32;
+ }else{
+ for(i=0;i<num_coefficients;i++)
+ rr[i] += d[i]*mul;
+ mul *=.5;
+ if(mul<1./(1<<6))break;
+ }
+ }
+
+ if(r_cost < cost){
+ flag = 1;
+ memcpy(r,rr,sizeof(rr));
+ cost = r_cost;
+ }
+ }*/
+
+ memcpy(prev,c,sizeof(prev));
+ memcpy(c,r,sizeof(r));
+ have_prev = 1;
+
+ fprintf(stderr,"\rwalking CSD... current cost: %g ",todB(cost)/2);
+ }
+ fprintf(stderr,"\n");
+ return 0;
}
double n_choose_m(int n, int m){
@@ -461,48 +438,99 @@
return ret;
}
+void set_maxflat(double *fc){
+ int i, j;
+ for(i=0,j=0;i<order;i++,j++)
+ fc[j] = def_a(order, i, order-1)+(drand48()-.5)*.01;
+
+ //for(i=0;i<order;i++,j++)
+ // fc[j] = def_a(order, order-i-1, order-1)+(drand48()-.5)*.01;
+ //fc[j++]=1.;
+ //for(i=0;i<order;i++,j++)
+ // fc[j] = def_a(order, i, order-1)+(drand48()-.5)*.01;
+ //for(i=0;i<order;i++,j++)
+ // fc[j] = def_a(order, order-i-1, order-1)+(drand48()-.5)*.01;
+ //fc[j++]=1.;
+}
+
int main(int argc, char *argv[]){
- int *c;
- double *fc;
- int i,j;
drft_init(&fft,FS);
- int order = 8;
- num_coefficients = (order*2+1)*2;
- TB = .1;
- p_delay = 7;
- q_delay = 8;
+ order = 16;
+ num_coefficients = order;//(order*2+1)*2;
+ TB = .20;
+ p_delay = 15;
+ q_delay = 16;
+ r_delay = 0;
+ s_delay = 0;
+
+ //double *fc = calloc(num_coefficients, sizeof(*fc));
- c = alloca(num_coefficients * sizeof(*c));
- fc = alloca(num_coefficients * sizeof(*fc));
+ //set_maxflat(fc);
+ //num_coefficients = (order*2+1)*2;
- fixed_bits = 8;
- for(i=0,j=0;i<order;i++,j++)
- fc[j] = def_a(order, i, order-1);
- for(i=0;i<order;i++,j++)
- fc[j] = def_a(order, order-i-1, order-1);
- fc[j++]=1.;
- for(i=0;i<order;i++,j++)
- fc[j] = def_a(order, i, order-1);
- for(i=0;i<order;i++,j++)
- fc[j] = def_a(order, order-i-1, order-1);
- fc[j++]=1.;
-
- for(i=0;i<num_coefficients;i++)
- c[i] = rint(fc[i] * pow(2, fixed_bits));
+ //walk_to_minimum_A(fc,1./(1<<10));
+ //walk_to_minimum_A(fc,1./(1<<12));
+ //walk_to_minimum_A(fc,1./(1<<14));
+ //walk_to_minimum_A(fc,1./(1<<16));
+ //walk_to_minimum_A(fc,1./(1<<18));
+
+ //lift_dump(fc);
+ //fflush(stdout);
- walk_to_minimum(c);
+ double fc[16]={ 0.494814, -0.117399, 0.053387, -0.029044, 0.016904, -0.010041,
+ 0.005928, -0.003415, 0.001890, -0.000988, 0.000478, -0.000208,
+ 0.000078, -0.000022, 0.000003, 0.000001 };
- for(j=9;j<=16;j++){
- fixed_bits = j;
- for(i=0;i<num_coefficients;i++)
- c[i] <<=1;
- walk_to_minimum(c);
- }
+ //set_maxflat(fc);
+ //walk_to_minimum_CSD(fc,1./(1<<16));
+ //lift_dump(fc);
+ //fflush(stdout);
+ //walk_to_minimum_CSD(fc,1./(1<<17));
+ //lift_dump(fc);
+ //fflush(stdout);
+ //walk_to_minimum_CSD(fc,1./(1<<18));
+ //lift_dump(fc);
+ //fflush(stdout);
+ //walk_to_minimum_CSD(fc,1./(1<<20));
+ //walk_to_minimum_CSD(fc,1./(1<<21));
+ //lift_dump(fc);
+ //fflush(stdout);
+ //walk_to_minimum_CSD(fc,1./(1<<22));
+ //lift_dump(fc);
+ //fflush(stdout);
+ //walk_to_minimum_CSD(fc,1./(1<<24));
+ //lift_dump(fc);
+ //fflush(stdout);
+ walk_to_minimum_CSD(fc,1./(1<<26));
+ //lift_dump(fc);
+ //fflush(stdout);
+ walk_to_minimum_CSD(fc,1./(1<<28));
+ //lift_dump(fc);
+ //fflush(stdout);
+ walk_to_minimum_CSD(fc,1./(1<<30));
+ //lift_dump(fc);
+ //fflush(stdout);
- lift_dump(c);
- lift_cost2(c);
+ lift_dump(fc);
+
return 0;
}
+
+#if 0
+Good ones so far:
+
+16/.21
+ a: ( 0.494814 -0.117399 0.053387 -0.029044 0.016904 -0.010041 0.005928 -0.003415 0.001890 -0.000988 0.000478 -0.000208 0.000078 -0.000022 0.000003 0.000001 )
+ a: ( 0.494801 -0.117383 0.053372 -0.029035 0.016902 -0.010047 0.005942 -0.003436 0.001913 -0.001012 0.000500 -0.000225 0.000090 -0.000030 0.000007 -0.000001 )
+16/.19
+ a: ( 0.495075 -0.117772 0.053811 -0.029476 0.017309 -0.010393 0.006214 -0.003631 0.002040 -0.001083 0.000531 -0.000232 0.000085 -0.000021 0.000000 0.000003 )
+16/.16
+ a: ( 0.496384 -0.119662 0.056015 -0.031800 0.019595 -0.012518 0.008089 -0.005205 0.003295 -0.002031 0.001205 -0.000679 0.000357 -0.000170 0.000070 -0.000022 )
+16/.16
+ a: ( 0.496302 -0.119542 0.055873 -0.031647 0.019441 -0.012371 0.007955 -0.005087 0.003196 -0.001951 0.001144 -0.000635 0.000327 -0.000152 0.000060 -0.000017 )
+
+
+#endif
Modified: trunk/ghost/monty/lift/smallft.c
===================================================================
--- trunk/ghost/monty/lift/smallft.c 2007-07-06 13:44:36 UTC (rev 13234)
+++ trunk/ghost/monty/lift/smallft.c 2007-07-06 14:21:34 UTC (rev 13235)
@@ -77,7 +77,7 @@
nfm1=nf-1;
l1=1;
- if(nfm1==0)return;
+ if(nfm1==0)return -1;
for (k1=0;k1<nfm1;k1++){
ip=ifac[k1+2];
Added: trunk/ghost/monty/sinusoid/scales.h
===================================================================
--- trunk/ghost/monty/sinusoid/scales.h (rev 0)
+++ trunk/ghost/monty/sinusoid/scales.h 2007-07-06 14:21:34 UTC (rev 13235)
@@ -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-80f)*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;
+}
+
Added: trunk/ghost/monty/sinusoid/sinusoids.c
===================================================================
--- trunk/ghost/monty/sinusoid/sinusoids.c (rev 0)
+++ trunk/ghost/monty/sinusoid/sinusoids.c 2007-07-06 14:21:34 UTC (rev 13235)
@@ -0,0 +1,289 @@
+/********************************************************************
+ * *
+ * 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 sode
+ last mod: $Id$
+
+ ********************************************************************/
+
+#include <math.h>
+#include "sinusoids.h"
+#include "scales.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+/* f: input; log magnitude (decibel scale) of input spectrum */
+
+void level_mean(float *f,float *out,int n, float lowindow, float hiwindow, int min, int rate){
+ float nevent[n];
+ float devent[n];
+
+ float nadd = 0.f;
+ float dadd = 0.f;
+ float nacc = 0.f;
+ float dacc = 0.f;
+
+ float binwidth = rate*.5f/n;
+ float ibinwidth = 1.f/binwidth;
+
+ int head=-1;
+ int tail=-1;
+ int i;
+
+ memset(nevent,0,sizeof(nevent));
+ memset(devent,0,sizeof(devent));
+
+ /* init case for hi-side window */
+ int lo;
+ int hi = rint(fromBark(hiwindow)*ibinwidth);
+ if(hi<min)hi=min;
+
+ for(head=0;head<hi;head++){
+ float de = 1./hi;
+ float nu = f[head]*de;
+
+ nevent[head] -= nu;
+ devent[head] -= de;
+
+ nadd += nu;
+ dadd += de;
+
+ nacc +=nadd;
+ dacc +=dadd;
+ }
+
+ for(i=0;i<n;i++){
+ float bark = toBark(i*binwidth);
+ hi = rint(fromBark(bark+hiwindow)*ibinwidth)-i;
+ if(bark>lowindow)
+ lo = i-rint(fromBark(bark-lowindow)*ibinwidth);
+ else
+ lo = rint(fromBark(lowindow)*ibinwidth);
+ if(hi<min)hi=min;
+ if(lo<min)lo=min;
+
+ /* high side window*/
+ for(;head<i+hi;head++){
+ if(head<n){
+ float de = 1./hi;
+ float nu = f[head]*de;
+
+ nevent[i+hi] -= nu;
+ devent[i+hi] -= de;
+
+ nadd += nu;
+ dadd += de;
+
+ }
+ }
+
+ /* low side window */
+ {
+ float de = 1./lo;
+ float nu = f[i]*de;
+
+ if(i+lo<n){
+ nevent[i+lo] += nu;
+ devent[i+lo] += de;
+ }
+ if(i<n){
+ nevent[i] -= nu;
+ devent[i] -= de;
+ }
+ }
+
+ nadd += nevent[i];
+ dadd += devent[i];
+
+ nacc += nadd;
+ dacc += dadd;
+
+ out[i] = nacc/dacc;
+ }
+}
+
+
+
+void extract_sinusoids(float *x, float *w, float *window, float *ai, float *bi, float *y, int N, int len)
+{
+ float cos_table[N][len];
+ float sin_table[N][len];
+ float cosE[N], sinE[N];
+ int i,j, iter;
+ for (i=0;i<N;i++)
+ {
+ float tmp1=0, tmp2=0;
+ for (j=0;j<len;j++)
+ {
+ cos_table[i][j] = cos(w[i]*j)*window[j];
+ sin_table[i][j] = sin(w[i]*j)*window[j];
+ tmp1 += cos_table[i][j]*cos_table[i][j];
+ tmp2 += sin_table[i][j]*sin_table[i][j];
+ }
+ cosE[i] = tmp1;
+ sinE[i] = tmp2;
+ }
+ for (j=0;j<len;j++)
+ y[j] = 0;
+ for (i=0;i<N;i++)
+ ai[i] = bi[i] = 0;
+
+ for (iter=0;iter<5;iter++)
+ {
+ for (i=0;i<N;i++)
+ {
+ float tmp1=0, tmp2=0;
+ for (j=0;j<len;j++)
+ {
+ tmp1 += (x[j]-y[j])*cos_table[i][j];
+ tmp2 += (x[j]-y[j])*sin_table[i][j];
+ }
+ tmp1 /= cosE[i];
+ //Just in case it's a DC! Must fix that anyway
+ tmp2 /= (.0001+sinE[i]);
+ for (j=0;j<len;j++)
+ {
+ y[j] += tmp1*cos_table[i][j] + tmp2*sin_table[i][j];
+ }
+ ai[i] += tmp1;
+ bi[i] += tmp2;
+ if (iter==4)
+ {
+ printf ("%f %f ", w[i], (float)sqrt(ai[i]*ai[i] + bi[i]*bi[i]));
+ }
+ }
+ }
+ printf ("\n");
+}
+
+/* Models the signal as modulated sinusoids
+ * x is the signal
+ * w are the frequencies
+ * window is the analysis window (you can use a rectangular window)
+ * ai are the cos(x) coefficients
+ * bi are the sin(x) coefficients
+ * ci are the x*cos(x) coefficients
+ * di are the x*sin(x) coefficients
+ * y is the approximated signal by summing all the params
+ * N is the number of sinusoids
+ * len is the frame size
+*/
+void extract_modulated_sinusoids(float *x, float *w, float *window, float *ai, float *bi, float *ci, float *di, float *y, int N, int len)
+{
+ float cos_table[N][len];
+ float sin_table[N][len];
+ float tcos_table[N][len];
+ float tsin_table[N][len];
+ float cosE[N], sinE[N];
+ float costE[N], sintE[N];
+ int i,j, iter;
+ /* Build a table for the four basis functions at each frequency: cos(x), sin(x), x*cos(x), x*sin(x)*/
+ for (i=0;i<N;i++)
+ {
+ float tmp1=0, tmp2=0;
+ float tmp3=0, tmp4=0;
+ for (j=0;j<len;j++)
+ {
+ float jj = j-len/2+.5;
+ cos_table[i][j] = cos(w[i]*jj)*window[j];
+ sin_table[i][j] = sin(w[i]*jj)*window[j];
+ tcos_table[i][j] = ((jj))*cos_table[i][j];
+ tsin_table[i][j] = ((jj))*sin_table[i][j];
+ /* The sinusoidal terms */
+ tmp1 += cos_table[i][j]*cos_table[i][j];
+ tmp2 += sin_table[i][j]*sin_table[i][j];
+ /* The modulation terms */
+ tmp3 += tcos_table[i][j]*tcos_table[i][j];
+ tmp4 += tsin_table[i][j]*tsin_table[i][j];
+ }
+ cosE[i] = sqrt(tmp1);
+ sinE[i] = sqrt(tmp2);
+ costE[i] = sqrt(tmp3);
+ sintE[i] = sqrt(tmp4);
+ for (j=0;j<len;j++)
+ {
+ cos_table[i][j] /= cosE[i];
+ sin_table[i][j] /= sinE[i];
+ tcos_table[i][j] /= costE[i];
+ tsin_table[i][j] /= sintE[i];
+ }
+ }
+ /* y is the initial approximation of the signal */
+ for (j=0;j<len;j++)
+ y[j] = 0;
+ for (i=0;i<N;i++)
+ ai[i] = bi[i] = ci[i] = di[i] = 0;
+ int tata=0;
+ /* This is an iterative solution -- much quicker than inverting a matrix */
+ for (iter=0;iter<5;iter++)
+ {
+ for (i=0;i<N;i++)
+ {
+ float tmp1=0, tmp2=0;
+ float tmp3=0, tmp4=0;
+ /* (Sort of) project the residual on the four basis functions */
+ for (j=0;j<len;j++)
+ {
+ tmp1 += (x[j]-y[j])*cos_table[i][j];
+ tmp2 += (x[j]-y[j])*sin_table[i][j];
+ tmp3 += (x[j]-y[j])*tcos_table[i][j];
+ tmp4 += (x[j]-y[j])*tsin_table[i][j];
+ }
+
+ //tmp3=tmp4 = 0;
+
+ /* Update the signal approximation for the next iteration */
+ for (j=0;j<len;j++)
+ {
+ y[j] += tmp1*cos_table[i][j] + tmp2*sin_table[i][j] + tmp3*tcos_table[i][j] + tmp4*tsin_table[i][j];
+ }
+ ai[i] += tmp1;
+ bi[i] += tmp2;
+ ci[i] += tmp3;
+ di[i] += tmp4;
+ if (iter==4)
+ {
+ if (w[i] > .49 && w[i] < .53 && !tata)
+ {
+ //printf ("%f %f %f %f %f\n", w[i], ai[i], bi[i], ci[i], di[i]);
+ tata = 1;
+ }
+ //printf ("%f %f ", w[i], (float)sqrt(ai[i]*ai[i] + bi[i]*bi[i]));
+ }
+ }
+ }
+ for (i=0;i<N;i++)
+ {
+ ai[i] /= cosE[i];
+ bi[i] /= sinE[i];
+ ci[i] /= costE[i];
+ di[i] /= sintE[i];
+ }
+#if 0
+ if (N)
+ for (i=0;i<1;i++)
+ {
+ float A, phi, dA, dw;
+ A = sqrt(ai[i]*ai[i] + bi[i]*bi[i]);
+ phi = atan2(bi[i], ai[i]);
+ //phi = ai[i]*ai[i] + bi[i]*bi[i];
+ dA = (ci[i]*ai[i] + bi[i]*di[i])/(.1+A);
+ dw = (ci[i]*bi[i] - di[i]*ai[i])/(.1+A*A);
+ printf ("%f %f %f %f %f %f %f %f %f\n", w[i], ai[i], bi[i], ci[i], di[i], A, phi, dA, dw);
+ }
+#endif
+ //if(!tata)
+ //printf ("0 0 0 0 0\n");
+
+ //printf ("\n");
+}
Added: trunk/ghost/monty/sinusoid/sinusoids.h
===================================================================
--- trunk/ghost/monty/sinusoid/sinusoids.h (rev 0)
+++ trunk/ghost/monty/sinusoid/sinusoids.h 2007-07-06 14:21:34 UTC (rev 13235)
@@ -0,0 +1,18 @@
+/********************************************************************
+ * *
+ * 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 sode
+ last mod: $Id$
+
+ ********************************************************************/
+
+extern void level_mean(float *f,float *out,int n, float lowindow, float hiwindow, int min, int rate);
Added: trunk/ghost/monty/sinusoid/smallft.c
===================================================================
--- trunk/ghost/monty/sinusoid/smallft.c (rev 0)
+++ trunk/ghost/monty/sinusoid/smallft.c 2007-07-06 14:21:34 UTC (rev 13235)
@@ -0,0 +1,335 @@
+/********************************************************************
+ * *
+ * 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];
+}
+
+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);
+}
+
+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));
+ }
+}
Added: trunk/ghost/monty/sinusoid/smallft.h
===================================================================
--- trunk/ghost/monty/sinusoid/smallft.h (rev 0)
+++ trunk/ghost/monty/sinusoid/smallft.h 2007-07-06 14:21:34 UTC (rev 13235)
@@ -0,0 +1,44 @@
+/********************************************************************
+ * *
+ * 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 int drft_init(drft_lookup *l,int n);
+extern void drft_clear(drft_lookup *l);
+
+#endif
Added: trunk/ghost/monty/sinusoid/test_sinusoids.c
===================================================================
--- trunk/ghost/monty/sinusoid/test_sinusoids.c (rev 0)
+++ trunk/ghost/monty/sinusoid/test_sinusoids.c 2007-07-06 14:21:34 UTC (rev 13235)
@@ -0,0 +1,114 @@
+#define _GNU_SOURCE
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "sinusoids.h"
+#include "smallft.h"
+#include "scales.h"
+
+#define BLOCK_SIZE 1024
+
+static float float_in[BLOCK_SIZE*2];
+static short short_in[BLOCK_SIZE];
+
+static void dump_vec(float *data, int n, char *base, int i){
+ char *filename;
+ FILE *f;
+
+ asprintf(&filename,"%s_%d.m",base,i);
+ f=fopen(filename,"w");
+
+ for(i=0;i<n;i++)
+ fprintf(f,"%f %f\n",
+ (float)i/n,data[i]);
+
+ fclose(f);
+
+}
+
+#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;
+ }
+}
+
+void mag_dB(float *d, int n){
+ int i;
+ d[0] = todB(d[0]*d[0])*.5;
+ for(i=2;i<n;i+=2)
+ d[i>>1] = todB(d[i-1]*d[i-1] + d[i]*d[i])*.5;
+ d[n/2] = todB(d[n-1]*d[n-1])*.5;
+}
+
+int main(int argc, char **argv){
+ FILE *fin, *fout;
+ int i,frame=0;
+ drft_lookup fft;
+
+ drft_init(&fft, BLOCK_SIZE*2);
+
+ if (argc != 3){
+ fprintf (stderr, "usage: testghost input_file output_file\nWhere the input and output are raw mono files sampled at 44.1 kHz or 48 kHz\n");
+ exit(1);
+ }
+
+ fin = fopen(argv[1], "r");
+ fout = fopen(argv[2], "w");
+
+ /*
+ for(i=0;i<BLOCK_SIZE*2;i++)
+ float_in[i]=1.;
+ blackmann_harris(float_in, float_in, BLOCK_SIZE*2);
+ dump_vec(float_in,BLOCK_SIZE*2,"window",0);
+ memset(float_in,0,sizeof(float_in));
+ */
+
+ while (1){
+ int i;
+ memmove(float_in,float_in+BLOCK_SIZE,sizeof(*float_in)*BLOCK_SIZE);
+ fread(short_in, sizeof(short), BLOCK_SIZE, fin);
+
+ if (feof(fin))
+ break;
+ for (i=0;i<BLOCK_SIZE;i++)
+ float_in[i+BLOCK_SIZE] = short_in[i] * .000030517578125;
+
+ {
+ /* generate a log spectrum */
+ float fft_buf[BLOCK_SIZE*2];
+ float weight[BLOCK_SIZE+1];
+
+ blackmann_harris(fft_buf, float_in, BLOCK_SIZE*2);
+ //dump_vec(float_in,BLOCK_SIZE*2,"data",frame);
+ dump_vec(fft_buf,BLOCK_SIZE*2,"win_data",frame);
+ drft_forward(&fft, fft_buf);
+ for(i=0;i<BLOCK_SIZE*2;i++)fft_buf[i] *= 2./BLOCK_SIZE;
+ //dump_vec(fft_buf,BLOCK_SIZE*2,"fft",frame);
+
+ mag_dB(fft_buf,BLOCK_SIZE*2);
+ dump_vec(fft_buf,BLOCK_SIZE+1,"logmag",frame);
+
+ level_mean(fft_buf,weight,BLOCK_SIZE+1, 1.f,1.f, 10, 44100);
+ dump_vec(weight,BLOCK_SIZE+1,"weight",frame);
+
+ }
+
+
+ //fwrite(short_in, sizeof(short), BLOCK_SIZE, fout);
+ frame++;
+ }
+
+ return 0;
+}
+
More information about the commits
mailing list