[xiph-commits] r13465 - trunk/ghost/monty/sinusoid

xiphmont at svn.xiph.org xiphmont at svn.xiph.org
Wed Aug 8 07:48:12 PDT 2007


Author: xiphmont
Date: 2007-08-08 07:48:12 -0700 (Wed, 08 Aug 2007)
New Revision: 13465

Modified:
   trunk/ghost/monty/sinusoid/sinusoids.c
   trunk/ghost/monty/sinusoid/test_sinusoids.c
Log:
More thinking aloud, nothing to see here.



Modified: trunk/ghost/monty/sinusoid/sinusoids.c
===================================================================
--- trunk/ghost/monty/sinusoid/sinusoids.c	2007-08-08 14:21:19 UTC (rev 13464)
+++ trunk/ghost/monty/sinusoid/sinusoids.c	2007-08-08 14:48:12 UTC (rev 13465)
@@ -271,56 +271,42 @@
   
 }
 
-#define A0 .35875f
-#define A1 .48829f
-#define A2 .14128f
-#define A3 .01168f
-
-static void blackmann_harris(float *d, int n){
+static void hanning(float *d, int n){
   int i;
   float scale = 2*M_PI/n;
 
   for(i=0;i<n;i++){
     float i5 = i+.5;
-    d[i] = A0 - A1*cos(scale*i5) + A2*cos(scale*i5*2) - A3*cos(scale*i5*3);
+    d[i] = .5-.5*cos(scale*i5);
   }
 }
 
-static void hanning(float *d, int n){
-  int i;
-  float scale = 2*M_PI/n;
+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;
 
-  for(i=0;i<n;i++){
-    float i5 = i+.5;
-    d[i] = .5-.5*cos(scale*i5);
+  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);
   }
 }
 
-/* Models the signal as modulated sinusoids 
- * x  : input signal
- * w  : input  window (such as Hanning)
- * A  : output vector of amplitudes
- * W  : input/output vector of frequencies
- * P  : output phase vector
- * dA : output amplitude modulation 
- * dW : input/output frequency modulation; null to prevent conv.
- * N  : input number of sinusoids
- * len: input frame size
-*/
+int extract_sinusoids_iteration(float *x, float *window, 
+				float *A, float *W, float *P, 
+				float *dA, float *dW, 
+				float *y, int N, int len){
 
-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){
   float scale = 2.*M_PI/len ;
   float iscale = len/(2.*M_PI);
-  int i,j;
+  int i,j,flag=0;
   float cwork[len];
   float swork[len];
 
-  memset(y,0,sizeof(*y)*len);
-
-  for (i=0;i<N;i++){
+  for(i=0;i<N;i++){
     
     float aP=0, bP=0, cP=0;
     float dP=0, eP=0, fP=0;
@@ -328,6 +314,9 @@
     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(j=0;j<len;j++){
       float jj = j-len/2.+.5;
       float jj2 = jj*jj;
@@ -348,8 +337,6 @@
       
       aP += (x[j]-y[j]) * cw;
       bP += (x[j]-y[j]) * sw;
-      cP += (x[j]-y[j]) * cw * jj;
-      dP += (x[j]-y[j]) * sw * jj;
 
       cwork[j] = c;
       swork[j] = s;
@@ -357,45 +344,35 @@
 
     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;
-      float cw = cwork[j]*window[j];
-      float sw = swork[j]*window[j];
 
-      y[j] += aP*cwork[j];
-      y[j] += bP*swork[j];
-      y[j] += cP*jj*cwork[j];
-      y[j] += dP*jj*swork[j];
+      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;
 
-      if(dW){
-	eP += (x[j]-y[j]) * cw * jj2;
-	fP += (x[j]-y[j]) * sw * jj2;
-      }
     }
-
-    if(dW){
-      
-      eP = (eE>0.f ? eP/eE : 0.f);
-      fP = (fE>0.f ? fP/fE : 0.f);
-      
-      for(j=0;j<len;j++){
-	float jj = j-len/2.+.5;
-	float jj2 = jj*jj;
-	
-	y[j] += eP*jj2*cwork[j];
-	y[j] += fP*jj2*swork[j];
-      }
-    }
-
-
+    
+    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;
       
       if(lA!=0.){
 	if(bP>0){
@@ -405,16 +382,43 @@
 	}
       }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]  += (cP*bP - dP*aP)/lA2*iscale;
-      if(dW)dW[i] += (eP*bP - fP*aP)/lA2*iscale;
+      W[i]  += lW;
+      dW[i] += ldW;
     }
+
+    compute_sinusoid(A[i],W[i],P[i],dA[i],dW[i],y,len);
+
   }
+
+  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, 
@@ -594,18 +598,3 @@
   
 }
 
-static void compute_sinusoids(float *A, float *W, float *P, 
-			      float *delA, float *delW, float *y, int N, int len){
-  int i,j;
-  float ilen = 1./len;
-
-  for(j=0;j<len;j++){
-    float jj = (j-len*.5+.5)*ilen;
-    y[j]=0.;
-    for(i=0;i<N;i++){
-      float a = A[i] + delA[i]*jj;
-      float w = W[i] + delW[i]*jj;
-      y[j] += a*cos(2.*M_PI*jj*w+P[i]);
-    }
-  }
-}

Modified: trunk/ghost/monty/sinusoid/test_sinusoids.c
===================================================================
--- trunk/ghost/monty/sinusoid/test_sinusoids.c	2007-08-08 14:21:19 UTC (rev 13464)
+++ trunk/ghost/monty/sinusoid/test_sinusoids.c	2007-08-08 14:48:12 UTC (rev 13465)
@@ -133,69 +133,75 @@
       float weight[BLOCK_SIZE+1];
 
       // polish the strongest peaks from weighting
-      if(frame==220){
-	float w[BLOCK_SIZE];
-	float Wout[BLOCK_SIZE];
-	float Aout[BLOCK_SIZE]={0};
-	float Pout[BLOCK_SIZE]={0};
-	float dAout[BLOCK_SIZE]={0};
-	float ddAout[BLOCK_SIZE]={0};
-	float dWout[BLOCK_SIZE]={0};
-	float y[BLOCK_SIZE*2];
-	float window[BLOCK_SIZE*2];
+      float w[BLOCK_SIZE];
+      float Wout[BLOCK_SIZE];
+      float Aout[BLOCK_SIZE]={0};
+      float Pout[BLOCK_SIZE]={0};
+      float dAout[BLOCK_SIZE]={0};
+      float ddAout[BLOCK_SIZE]={0};
+      float dWout[BLOCK_SIZE]={0};
+      float y[BLOCK_SIZE*2];
+      float window[BLOCK_SIZE*2];
+      
+      hanning(fft_buf, float_in, BLOCK_SIZE*2);
+      //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);
+      
+      window_weight(log_fft,weight,BLOCK_SIZE+1, 0.f, 512,256, 30, 44100);
+      //dump_vec(weight,BLOCK_SIZE+1,"weight",frame);
 
-	blackmann_harris(fft_buf, float_in, BLOCK_SIZE*2);
-	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;
+      //j=2;
+      //	w[0]=.0044*BLOCK_SIZE;
+      //	w[1]=.136*BLOCK_SIZE;
 
-	mag_dB(log_fft,fft_buf,BLOCK_SIZE*2);
-	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);
-
-	j=2;
-	w[0]=.0306*BLOCK_SIZE;
-	w[1]=.136*BLOCK_SIZE;
-
-	//for(j=0;j<BLOCK_SIZE;j++)
-	//w[j] = j+.5;
-	  
-	/* largest weighted
-	   int best=-120;
-	  int besti=-1;
-	  for(i=0;i<BLOCK_SIZE+1;i++){
+      int j,k;
+      for(j=0;j<20;j++){
+	/* largest weighted */
+	int best=-120;
+	int besti=-1;
+	for(i=0;i<BLOCK_SIZE+1;i++){
 	  float v = log_fft[i] - weight[i];
-	    if(v>best){
-	    besti = i;
+	  if(v>best){
+	    for(k=0;k<j;k++)
+	      if(i == w[j])break;
+	    if(k==j){
+	      besti = i;
 	      best = v;
-	      }
 	    }
-	  //w[j] = M_PI*besti/BLOCK_SIZE;
-	  */
+	  }
+	}
+	w[j] = besti;
+      }
 
-	//blackmann_harris(float_in, float_in, BLOCK_SIZE*2);
-	/*hanningW(window,BLOCK_SIZE*2);
+      //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_modulated_sinusoidsB(float_in, w, Aout, Wout, Pout, dAout, dWout, ddAout, y, j, BLOCK_SIZE*2);
-	w[0]=Wout[0];
+	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);
-	w[0]=Wout[0];
+	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);
-	w[0]=Wout[0];
+	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++)
-	//w[i]=Wout[i];
+	for(i=0;i<j;i++)
+	  fprintf(stdout, "%d %f\n\n",frame,Wout[i]/BLOCK_SIZE*22050);
+
 	
-	for(i=0;i<BLOCK_SIZE*2;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);
@@ -217,7 +223,7 @@
 	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;
@@ -234,7 +240,7 @@
 	//exit(0);
 
 	//}
-      }
+	//}
       
     }
   



More information about the commits mailing list