[xiph-commits] r12956 - trunk/ghost/monty/lift
xiphmont at svn.xiph.org
xiphmont at svn.xiph.org
Wed May 16 10:32:30 PDT 2007
Author: xiphmont
Date: 2007-05-16 10:32:29 -0700 (Wed, 16 May 2007)
New Revision: 12956
Modified:
trunk/ghost/monty/lift/liftsearch.c
Log:
Random additional working cost eval code, not really used for anything
yet
Modified: trunk/ghost/monty/lift/liftsearch.c
===================================================================
--- trunk/ghost/monty/lift/liftsearch.c 2007-05-16 17:22:31 UTC (rev 12955)
+++ trunk/ghost/monty/lift/liftsearch.c 2007-05-16 17:32:29 UTC (rev 12956)
@@ -93,11 +93,12 @@
}
-void lift_run(int *c, double *lp, double *hp, int dump){
+int lift_run(int *c, double *lp, double *hp, int dump){
iir_t P;
iir_t Q;
int i,j;
- int offset = (abs(p_delay) + abs(q_delay))*2+1;
+ int offset = (abs(p_delay) + abs(q_delay))*2+3;
+ int ret = 0;
memset(&P,0,sizeof(P));
memset(&Q,0,sizeof(Q));
@@ -130,6 +131,16 @@
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;
+ }
+
// energy
mag_e(hp);
mag_e(lp);
@@ -170,6 +181,13 @@
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
@@ -181,7 +199,7 @@
for(i=TBhi; i<=FS/2; i++)
cost += lp[i];
for(i=0; i<=TBlo; i++)
- cost += hp[i];
+ cost += hp[i]*.25;
// deviation from 0dB in passband
//for(i=0; i<=TBlo; i++)
@@ -200,9 +218,176 @@
free(hp);
free(lp);
- return cost;
+ return todB(sqrt(cost/(FS-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){
double *hp = calloc (FS,sizeof(*hp));
double *lp = calloc (FS,sizeof(*lp));
@@ -215,7 +400,7 @@
}
int walk_to_minimum(int *c){
- int last_change = 0;
+ int last_change = num_coefficients-1;
int cur_dim = 0;
double cost = lift_cost(c);
@@ -231,36 +416,11 @@
last_change = cur_dim;
cost = up_cost;
- // keep walking until we hit the minimum in this dim
- /*while(1){
- c[cur_dim]++;
- up_cost = lift_cost(c);
- if(up_cost<cost){
- cost = up_cost;
- }else{
- c[cur_dim]--;
- break;
- }
- }*/
-
}else if (down_cost<cost){
c[cur_dim]--;
last_change = cur_dim;
cost = down_cost;
- // keep walking until we hit the minimum in this dim
- /*
- while(1){
- c[cur_dim]--;
- down_cost = lift_cost(c);
- if(down_cost<cost){
- cost = down_cost;
- }else{
- c[cur_dim]++;
- break;
- }
- }*/
-
}else{
if(cur_dim == last_change)break;
}
@@ -271,63 +431,11 @@
fprintf(stderr,"\rwalking... current cost: %g ",cost);
}
- fprintf(stderr,"\r ",cost);
+ fprintf(stderr,"\rfinal cost: %f \n",cost);
log_minimum(c,cost);
}
-
-int walk_to_minimum2(int *c){
- double cost = lift_cost(c);
- int i;
-
- while(1){
- double up_cost[num_coefficients];
- double down_cost[num_coefficients];
- double min_cost = cost;
- int min_cost_i = -1;
- int min_cost_d = -1;
-
- for(i=0;i<num_coefficients;i++){
- c[i]++;
- up_cost[i] = lift_cost(c);
- c[i]-=2;
- down_cost[i] = lift_cost(c);
- c[i]++;
- }
-
- for(i=0;i<num_coefficients;i++){
- if(up_cost[i] < cost){
- cost = up_cost[i];
- min_cost_i = i;
- min_cost_d = 1;
- }
- if(down_cost[i] < cost){
- cost = down_cost[i];
- min_cost_i = i;
- min_cost_d = -1;
- }
- }
-
- if(min_cost_i == -1) break;
- c[min_cost_i] += min_cost_d;
-
- fprintf(stderr,"\rwalking... current cost: %f ",cost);
- }
- fprintf(stderr,"\r ",cost);
-
- log_minimum(c,cost);
-}
-
-
-void bisect_Ncube_recurse(int *lowrange, int *highrange, int min_num){
-
-
-
-
-
-}
-
double n_choose_m(int n, int m){
double work[n+1];
int i,j;
@@ -363,22 +471,22 @@
int order = 8;
num_coefficients = (order*2+1)*2;
TB = .1;
- p_delay = order-1;
- q_delay = order;
+ p_delay = 7;
+ q_delay = 8;
c = alloca(num_coefficients * sizeof(*c));
fc = alloca(num_coefficients * sizeof(*fc));
fixed_bits = 8;
for(i=0,j=0;i<order;i++,j++)
- fc[j] = def_a(order, i, p_delay);
+ fc[j] = def_a(order, i, order-1);
for(i=0;i<order;i++,j++)
- fc[j] = def_a(order, order-i-1, p_delay);
+ 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, p_delay);
+ fc[j] = def_a(order, i, order-1);
for(i=0;i<order;i++,j++)
- fc[j] = def_a(order, order-i-1, p_delay);
+ fc[j] = def_a(order, order-i-1, order-1);
fc[j++]=1.;
for(i=0;i<num_coefficients;i++)
@@ -386,7 +494,7 @@
walk_to_minimum(c);
- for(j=9;j<=20;j++){
+ for(j=9;j<=16;j++){
fixed_bits = j;
for(i=0;i<num_coefficients;i++)
c[i] <<=1;
@@ -394,6 +502,7 @@
}
lift_dump(c);
+ lift_cost2(c);
return 0;
}
More information about the commits
mailing list