[xiph-commits] r13511 - trunk/ghost/monty/sinusoid
xiphmont at svn.xiph.org
xiphmont at svn.xiph.org
Fri Aug 10 07:27:59 PDT 2007
Author: xiphmont
Date: 2007-08-10 07:27:59 -0700 (Fri, 10 Aug 2007)
New Revision: 13511
Modified:
trunk/ghost/monty/sinusoid/sinusoids.c
trunk/ghost/monty/sinusoid/sinusoids.h
trunk/ghost/monty/sinusoid/test_sinusoids.c
Log:
Latest monty version of the nonlinear code. Done with it for now, back
to concentrating on picking seeds well.
Modified: trunk/ghost/monty/sinusoid/sinusoids.c
===================================================================
--- trunk/ghost/monty/sinusoid/sinusoids.c 2007-08-10 14:21:57 UTC (rev 13510)
+++ trunk/ghost/monty/sinusoid/sinusoids.c 2007-08-10 14:27:59 UTC (rev 13511)
@@ -15,6 +15,7 @@
********************************************************************/
+#define _GNU_SOURCE
#include <math.h>
#include "sinusoids.h"
#include "scales.h"
@@ -180,7 +181,7 @@
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 (iter=0;iter<50;iter++)
{
for (i=0;i<N;i++)
{
@@ -281,144 +282,155 @@
}
}
-static void compute_sinusoid(float A, float W, float P,
- float delA, float delW,
- float *y, int len){
- int i,j;
- float ilen = 1./len;
+int extract_sinusoids(float *x,
+ float *A, float *W, float *P,
+ float *dA, float *dW,
+ float *y, int N, int len, int count){
- for(j=0;j<len;j++){
- float jj = (j-len*.5+.5)*ilen;
- float a = A + delA*jj;
- float w = W + delW*jj;
- y[j] += a*cos(2.*M_PI*jj*w+P);
- }
-}
-
-int extract_sinusoids_iteration(float *x, float *window,
- float *A, float *W, float *P,
- float *dA, float *dW,
- float *y, int N, int len){
-
- float scale = 2.*M_PI/len ;
+ float scale = 2.*M_PI/len;
float iscale = len/(2.*M_PI);
- int i,j,flag=0;
+ float ilen = 1./len;
+ int i,j,flag=1;
float cwork[len];
float swork[len];
+ float window[len];
- for(i=0;i<N;i++){
+ float aprev[N],bprev[N];
+ float cprev[N],dprev[N];
+ float eprev[N],fprev[N];
+
+ hanning(window,len);
+ memset(y,0,sizeof(*y)*len);
+ memset(aprev,0,sizeof(aprev));
+ memset(bprev,0,sizeof(bprev));
+ memset(cprev,0,sizeof(cprev));
+ memset(dprev,0,sizeof(dprev));
+ memset(eprev,0,sizeof(eprev));
+ memset(fprev,0,sizeof(fprev));
+
+ while(count-- && flag){
+ flag=0;
- float aP=0, bP=0, cP=0;
- float dP=0, eP=0, fP=0;
- float aE=0, bE=0, cE=0;
- float dE=0, eE=0, fE=0;
- float ldW = (dW?dW[i]:0.f);
-
- if(A[i]!=0. || dA[i]!=0.)
- compute_sinusoid(-A[i],W[i],P[i],-dA[i],dW[i],y,len);
+ for(i=0;i<N;i++){
+
+ 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(P[i]);
+ float aS=sin(P[i]);
+
+ for(j=0;j<len;j++){
+ float jj = j-len*.5+.5;
+ float jj2 = jj*jj;
+ float jj4 = jj2*jj2;
+ float a = A[i] + dA[i]*jj*ilen;
+ float c = cos((W[i] + dW[i]*jj*ilen)*scale*jj);
+ float s = sin((W[i] + dW[i]*jj*ilen)*scale*jj);
+ float cw = c*window[j];
+ float sw = s*window[j];
+ float ccw = cw*c;
+ float ssw = sw*s;
+
+ cwork[j] = c;
+ swork[j] = s;
+ y[j] -= aC*a*c - aS*a*s;
+
+ aE += ccw;
+ bE += ssw;
+ cE += ccw*jj2;
+ dE += ssw*jj2;
+ eE += ccw*jj4;
+ fE += ssw*jj4;
- for(j=0;j<len;j++){
- float jj = j-len/2.+.5;
- float jj2 = jj*jj;
- float jj4 = jj2*jj2;
- float c = cos((W[i]*scale + ldW*scale*jj)*jj);
- float s = sin((W[i]*scale + ldW*scale*jj)*jj);
- float cw = c*window[j];
- float sw = s*window[j];
- float ccw = cw*c;
- float ssw = sw*s;
-
- aE += ccw;
- bE += ssw;
- cE += ccw*jj2;
- dE += ssw*jj2;
- eE += ccw*jj4;
- fE += ssw*jj4;
+ aP += (x[j]-y[j]) * cw;
+ bP += (x[j]-y[j]) * sw;
+
+ }
- aP += (x[j]-y[j]) * cw;
- bP += (x[j]-y[j]) * sw;
+ aP = (aE>0.f ? aP/aE : 0.f);
+ bP = (bE>0.f ? bP/bE : 0.f);
+
+ for(j=0;j<len;j++){
+ float jj = j-len/2.+.5;
+
+ cP += (x[j]-y[j]-aP*cwork[j]) * cwork[j] * window[j] * jj;
+ dP += (x[j]-y[j]-bP*cwork[j]) * swork[j] * window[j] * jj;
+ }
+
+ cP = (cE>0.f ? cP/cE : 0.f);
+ dP = (dE>0.f ? dP/dE : 0.f);
+
+ for(j=0;j<len;j++){
+ float jj = j-len/2.+.5;
+ float jj2 = jj*jj;
+
+ eP += (x[j]-y[j]-(aP+cP*jj)*cwork[j]) * cwork[j]*window[j] * jj2;
+ fP += (x[j]-y[j]-(bP+dP*jj)*swork[j]) * swork[j]*window[j] * jj2;
+ }
+
+ eP = (eE>0.f ? eP/eE : 0.f);
+ fP = (fE>0.f ? fP/fE : 0.f);
+
+ {
+ float lA2 = aP*aP + bP*bP;
+ float lA = sqrt(aP*aP + bP*bP);
+ float lP;
+
+ if(lA!=0.){
+ if(bP>0){
+ lP = -acos(aP/lA);
+ }else{
+ lP = acos(aP/lA);
+ }
+ }else
+ lP=0.;
+
+ A[i] = lA;
+ P[i] = lP;
+ dA[i] = (cP*aP + bP*dP)/lA;
+ W[i] += (cP*bP - dP*aP)/lA2*iscale*.8;
+ dW[i] += (eP*bP - fP*aP)/lA2*iscale*.8;
+ }
+
+ if(todB(fabs(aP - aprev[i]))>-120. ||
+ todB(fabs(bP - bprev[i]))>-120. ||
+ todB(fabs(cP - cprev[i]))>-120. ||
+ todB(fabs(dP - dprev[i]))>-120. ||
+ todB(fabs(eP - eprev[i]))>-120. ||
+ todB(fabs(fP - fprev[i]))>-120.) flag=1;
- cwork[j] = c;
- swork[j] = s;
- }
-
- aP = (aE>0.f ? aP/aE : 0.f);
- bP = (bE>0.f ? bP/bE : 0.f);
-
- for(j=0;j<len;j++){
- float jj = j-len/2.+.5;
-
- cP += (x[j]-y[j]-aP*cwork[j]) * cwork[j] * window[j] * jj;
- dP += (x[j]-y[j]-bP*cwork[j]) * swork[j] * window[j] * jj;
- }
-
- cP = (cE>0.f ? cP/cE : 0.f);
- dP = (dE>0.f ? dP/dE : 0.f);
-
- for(j=0;j<len;j++){
- float jj = j-len/2.+.5;
- float jj2 = jj*jj;
-
- eP += (x[j]-y[j]-(aP+cP*jj)*cwork[j]) * cwork[j]*window[j] * jj2;
- fP += (x[j]-y[j]-(bP+dP*jj)*swork[j]) * swork[j]*window[j] * jj2;
-
- }
-
- eP = (eE>0.f ? eP/eE : 0.f);
- fP = (fE>0.f ? fP/fE : 0.f);
-
- {
- float lA2 = aP*aP + bP*bP;
- float lA = sqrt(aP*aP + bP*bP);
- float lP;
- float lW;
- float ldW;
+ fprintf(stderr,"%.0f %.0f %.0f %.0f %.0f %.0f\n",
+ todB(fabs(aP - aprev[i])),
+ todB(fabs(bP - bprev[i])),
+ todB(fabs(cP - cprev[i])),
+ todB(fabs(dP - dprev[i])),
+ todB(fabs(eP - eprev[i])),
+ todB(fabs(fP - fprev[i])));
- if(lA!=0.){
- if(bP>0){
- lP = -acos(aP/lA);
- }else{
- lP = acos(aP/lA);
- }
- }else
- lP=0.;
-
- lW = (cP*bP - dP*aP)/lA2*iscale;
- ldW = (eP*bP - fP*aP)/lA2*iscale;
-
- if( fabs(lW)>.001 || fabs(ldW)>.001) flag=1;
-
- A[i] = lA;
- P[i] = lP;
- dA[i] = (cP*aP + bP*dP)/lA;
- W[i] += lW;
- dW[i] += ldW;
+ aprev[i] = aP;
+ bprev[i] = bP;
+ cprev[i] = cP;
+ dprev[i] = dP;
+ eprev[i] = eP;
+ fprev[i] = fP;
+
+ for(j=0;j<len;j++){
+ float jj = (j-len*.5+.5)*ilen;
+ float a = A[i] + dA[i]*jj;
+ float w = W[i] + dW[i]*jj;
+ y[j] += a*cos(2.*M_PI*jj*w+P[i]);
+ }
}
- compute_sinusoid(A[i],W[i],P[i],dA[i],dW[i],y,len);
-
}
-
+ fprintf(stderr,"count remaining=%d\n",count);
return flag;
}
-void extract_modulated_sinusoids_nonlinear(float *x, float *w,
- float *A, float *W, float *P,
- float *dA, float *dW,
- float *y, int N, int len){
- float window[len];
- hanning(window,len);
- memset(y,0,sizeof(*y)*len);
- memset(A,0,sizeof(*A)*N);
- memset(P,0,sizeof(*A)*N);
- memset(dA,0,sizeof(*A)*N);
- memset(dW,0,sizeof(*A)*N);
- int count=5;
-
- while(extract_sinusoids_iteration(x, window, A, W, P, dA, dW, y, N, len) && --count);
-
-}
-
void extract_modulated_sinusoidsB(float *x, float *w,
float *Aout, float *Wout, float *Pout,
float *dAout, float *dWout, float *ddAout,
@@ -442,7 +454,7 @@
int i,j, iter;
hanning(window,len);
-
+
/* Build a table for the four basis functions at each frequency */
for (i=0;i<N;i++){
float tmpa=0;
@@ -451,14 +463,14 @@
float tmpd=0;
float tmpe=0;
float tmpf=0;
-
+
cos_table[i]=calloc(len,sizeof(**cos_table));
sin_table[i]=calloc(len,sizeof(**sin_table));
tcos_table[i]=calloc(len,sizeof(**tcos_table));
tsin_table[i]=calloc(len,sizeof(**tsin_table));
ttcos_table[i]=calloc(len,sizeof(**ttcos_table));
ttsin_table[i]=calloc(len,sizeof(**ttsin_table));
-
+
for (j=0;j<len;j++){
float jj = j-len/2.+.5;
@@ -492,14 +504,14 @@
tmpd = sqrt(tmpd);
tmpe = sqrt(tmpe);
tmpf = sqrt(tmpf);
-
+
cosE[i] = (tmpa>0.f ? 1.f/tmpa : 0.f);
sinE[i] = (tmpb>0.f ? 1.f/tmpb : 0.f);
tcosE[i] = (tmpc>0.f ? 1.f/tmpc : 0.f);
tsinE[i] = (tmpd>0.f ? 1.f/tmpd : 0.f);
ttcosE[i] = (tmpe>0.f ? 1.f/tmpe : 0.f);
ttsinE[i] = (tmpf>0.f ? 1.f/tmpf : 0.f);
-
+
for(j=0;j<len;j++){
cos_table[i][j] *= cosE[i];
sin_table[i][j] *= sinE[i];
@@ -518,14 +530,14 @@
for (iter=0;iter<5;iter++){
for (i=0;i<N;i++){
-
+
float tmpa=0, tmpb=0, tmpc=0;
float tmpd=0, tmpe=0, tmpf=0;
-
+
/* (Sort of) project the residual on the four basis functions */
for (j=0;j<len;j++){
- float jj = j-len/2.+.5;
-
+ float jj = j-len/2.+.5;
+
tmpa += (x[j]-y[j]) * cos_table[i][j] * window[j];
tmpb += (x[j]-y[j]) * sin_table[i][j] * window[j];
tmpc += (x[j]-y[j]) * tcos_table[i][j] * window[j];
@@ -538,15 +550,12 @@
/* Update the signal approximation for the next iteration */
for (j=0;j<len;j++){
- float jj = j-len/2.+.5;
-
y[j] += tmpa*cos_table[i][j];
y[j] += tmpb*sin_table[i][j];
y[j] += tmpc*tcos_table[i][j];
y[j] += tmpd*tsin_table[i][j];
y[j] += tmpe*ttcos_table[i][j];
y[j] += tmpf*ttsin_table[i][j];
-
}
ai[i] += tmpa;
@@ -555,10 +564,9 @@
di[i] += tmpd;
ei[i] += tmpe;
fi[i] += tmpf;
-
}
}
-
+
for (i=0;i<N;i++){
ai[i] *= cosE[i];
bi[i] *= sinE[i];
@@ -578,12 +586,12 @@
free(tsin_table[i]);
free(ttcos_table[i]);
free(ttsin_table[i]);
-
+
if(A!=0.){
if(bi[i]>0){
P = -acos(ai[i]/A);
}else{
- P = acos(ai[i]/A);
+ P = acos(ai[i]/A);
}
}else
P=0.;
@@ -597,4 +605,4 @@
}
}
-
+
Modified: trunk/ghost/monty/sinusoid/sinusoids.h
===================================================================
--- trunk/ghost/monty/sinusoid/sinusoids.h 2007-08-10 14:21:57 UTC (rev 13510)
+++ trunk/ghost/monty/sinusoid/sinusoids.h 2007-08-10 14:27:59 UTC (rev 13511)
@@ -18,6 +18,11 @@
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);
Modified: trunk/ghost/monty/sinusoid/test_sinusoids.c
===================================================================
--- trunk/ghost/monty/sinusoid/test_sinusoids.c 2007-08-10 14:21:57 UTC (rev 13510)
+++ trunk/ghost/monty/sinusoid/test_sinusoids.c 2007-08-10 14:27:59 UTC (rev 13511)
@@ -125,7 +125,7 @@
for (i=0;i<BLOCK_SIZE;i++)
float_in[i+BLOCK_SIZE] = short_in[i] * .000030517578125;
- {
+ if(frame==220){
/* generate a log spectrum */
float fft_buf[BLOCK_SIZE*2];
@@ -144,23 +144,22 @@
float window[BLOCK_SIZE*2];
hanning(fft_buf, float_in, BLOCK_SIZE*2);
- //dump_vec(float_in,BLOCK_SIZE*2,"data",frame);
+ dump_vec(float_in,BLOCK_SIZE*2,"data",frame);
drft_forward(&fft, fft_buf);
for(i=0;i<BLOCK_SIZE*2;i++)fft_buf[i] *= 1./BLOCK_SIZE;
mag_dB(log_fft,fft_buf,BLOCK_SIZE*2);
- //dump_vec(log_fft,BLOCK_SIZE+1,"logmag",frame);
+ dump_vec(log_fft,BLOCK_SIZE+1,"logmag",frame);
window_weight(log_fft,weight,BLOCK_SIZE+1, 0.f, 512,256, 30, 44100);
- //dump_vec(weight,BLOCK_SIZE+1,"weight",frame);
+ dump_vec(weight,BLOCK_SIZE+1,"weight",frame);
- //j=2;
- // w[0]=.0044*BLOCK_SIZE;
- // w[1]=.136*BLOCK_SIZE;
+ j=2;
+ w[0]=.0044*BLOCK_SIZE;
+ w[1]=.136*BLOCK_SIZE;
- int j,k;
- for(j=0;j<20;j++){
- /* largest weighted */
+ /*int j,k;
+ for(j=0;j<20;j++){
int best=-120;
int besti=-1;
for(i=0;i<BLOCK_SIZE+1;i++){
@@ -175,16 +174,11 @@
}
}
w[j] = besti;
- }
+ }*/
- //blackmann_harris(float_in, float_in, BLOCK_SIZE*2);
- /*hanningW(window,BLOCK_SIZE*2);
- extract_modulated_sinusoids_nonlinear(float_in, window, Aout, w, Pout, dAout, dWout, y, j, BLOCK_SIZE*2);
- extract_modulated_sinusoids_nonlinear(float_in, window, Aout, w, Pout, dAout, dWout, y, j, BLOCK_SIZE*2);
- extract_modulated_sinusoids_nonlinear(float_in, window, Aout, w, Pout, dAout, dWout, y, j, BLOCK_SIZE*2);
- extract_modulated_sinusoids_nonlinear(float_in, window, Aout, w, Pout, dAout, dWout, y, j, BLOCK_SIZE*2);
- extract_modulated_sinusoids_nonlinear(float_in, window, Aout, w, Pout, dAout, dWout, y, j, BLOCK_SIZE*2);*/
-
+
+ extract_sinusoids(float_in, Aout, w, Pout, dAout, dWout, y, j, BLOCK_SIZE*2, 200);
+ /*
extract_modulated_sinusoidsB(float_in, w, Aout, Wout, Pout, dAout, dWout, ddAout, y, j, BLOCK_SIZE*2);
for(i=0;i<j;i++)
w[i]=Wout[i];
@@ -195,52 +189,55 @@
for(i=0;i<j;i++)
w[i]=Wout[i];
extract_modulated_sinusoidsB(float_in, w, Aout, Wout, Pout, dAout, dWout, ddAout, y, j, BLOCK_SIZE*2);
-
-
for(i=0;i<j;i++)
- fprintf(stdout, "%d %f\n\n",frame,Wout[i]/BLOCK_SIZE*22050);
+ w[i]=Wout[i];
+ */
-
- /* for(i=0;i<BLOCK_SIZE*2;i++)
- fft_buf[i] = float_in[i]-y[i];
-
- dump_vec(fft_buf,BLOCK_SIZE*2,"res",0);
- hanning(fft_buf, fft_buf, BLOCK_SIZE*2);
- drft_forward(&fft, fft_buf);
- for(i=0;i<BLOCK_SIZE*2;i++)fft_buf[i] *= 1./BLOCK_SIZE;
- mag_dB(log_fft,fft_buf,BLOCK_SIZE*2);
-
- dump_vec(log_fft,BLOCK_SIZE+1,"exlogmag",0);
-
- for(i=0;i<j;i++){
- Aout[i]=todB(Aout[i]);
- dAout[i]=todB(dAout[i]);
- w[i] /= BLOCK_SIZE;
- }
-
- dump_vec2(w,Aout,j,"ex",0);
- dump_vec2(w,dAout,j,"dA",0);
- dump_vec2(w,dWout,j,"dW",0);
- dump_vec(y,BLOCK_SIZE*2,"extract",0);
- */
- //for (i=0;i<BLOCK_SIZE*2;i++)
- // float_out[i+outcount] += fft_buf[i];
- //outcount+=BLOCK_SIZE;
-
- //if(outcount+BLOCK_SIZE*2>262144){
-
- // hanning(float_out, float_out, 262144);
- //drft_init(&fft, 262144);
- //drft_forward(&fft, float_out);
- //for(i=0;i<262144;i++)float_out[i] *= 1./262144;
-
- //mag_dB(float_out,float_out,262144);
- //dump_vec(float_out,262144/2,"res",0);
- //exit(0);
-
+ for(i=0;i<j;i++)
+ fprintf(stdout, "%d %f\n\n",frame,w[i]/BLOCK_SIZE*22050);
+
+
+ for(i=0;i<BLOCK_SIZE*2;i++)
+ fft_buf[i] = float_in[i]-y[i];
+
+ dump_vec(fft_buf,BLOCK_SIZE*2,"res",0);
+
+ hanning(fft_buf, fft_buf, BLOCK_SIZE*2);
+ drft_forward(&fft, fft_buf);
+ for(i=0;i<BLOCK_SIZE*2;i++)fft_buf[i] *= 1./BLOCK_SIZE;
+ mag_dB(log_fft,fft_buf,BLOCK_SIZE*2);
+
+ dump_vec(log_fft,BLOCK_SIZE+1,"exlogmag",0);
+
+ for(i=0;i<j;i++){
+ Aout[i]=todB(Aout[i]);
+ dAout[i]=todB(dAout[i]);
+ w[i] /= BLOCK_SIZE;
+ }
+
+ dump_vec2(w,Aout,j,"ex",0);
+ dump_vec2(w,dAout,j,"dA",0);
+ dump_vec2(w,dWout,j,"dW",0);
+ dump_vec(y,BLOCK_SIZE*2,"extract",0);
+
+ //for (i=0;i<BLOCK_SIZE*2;i++)
+ // float_out[i+outcount] += fft_buf[i];
+ //outcount+=BLOCK_SIZE;
+
+ //if(outcount+BLOCK_SIZE*2>262144){
+
+ // hanning(float_out, float_out, 262144);
+ //drft_init(&fft, 262144);
+ //drft_forward(&fft, float_out);
+ //for(i=0;i<262144;i++)float_out[i] *= 1./262144;
+
+ //mag_dB(float_out,float_out,262144);
+ //dump_vec(float_out,262144/2,"res",0);
+ //exit(0);
+
+ //}
//}
- //}
}
More information about the commits
mailing list