[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