[xiph-commits] r13375 - trunk/ghost/monty/sinusoid
xiphmont at svn.xiph.org
xiphmont at svn.xiph.org
Fri Jul 27 18:39:06 PDT 2007
Author: xiphmont
Date: 2007-07-27 18:39:06 -0700 (Fri, 27 Jul 2007)
New Revision: 13375
Modified:
trunk/ghost/monty/sinusoid/scales.h
trunk/ghost/monty/sinusoid/sinusoids.c
trunk/ghost/monty/sinusoid/sinusoids.h
trunk/ghost/monty/sinusoid/test_sinusoids.c
Log:
More thought process.
Modified: trunk/ghost/monty/sinusoid/scales.h
===================================================================
--- trunk/ghost/monty/sinusoid/scales.h 2007-07-27 20:57:57 UTC (rev 13374)
+++ trunk/ghost/monty/sinusoid/scales.h 2007-07-28 01:39:06 UTC (rev 13375)
@@ -21,7 +21,7 @@
#define MAX(a,b) ((a)>(b) ? (a):(b))
static inline float todB(float x){
- return log(x*x+1e-80f)*4.34294480f;
+ return log(x*x+1e-24f)*4.34294480f;
}
static inline float fromdB(float x){
Modified: trunk/ghost/monty/sinusoid/sinusoids.c
===================================================================
--- trunk/ghost/monty/sinusoid/sinusoids.c 2007-07-27 20:57:57 UTC (rev 13374)
+++ trunk/ghost/monty/sinusoid/sinusoids.c 2007-07-28 01:39:06 UTC (rev 13375)
@@ -24,9 +24,9 @@
void window_weight(float *logf, float *out, int n, float flatbias,
int lowindow, int hiwindow, int min, int rate){
+
int bark[n],i;
float binwidth = rate*.5f/n;
- float ibinwidth = 1.f/binwidth;
for(i=0;i<n;i++)
bark[i] = rint(toBark(i*binwidth) * 1024);
@@ -102,7 +102,7 @@
{
float arith = todB(aacc / dacc);
float geom = gacc / dacc;
- out[i] = logf[i] - arith + flatbias*(arith-geom);
+ out[i] = arith - flatbias*(arith-geom);
}
}
@@ -112,7 +112,6 @@
/* 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
@@ -121,16 +120,21 @@
* N is the number of sinusoids
* len is the frame size
*/
-void extract_modulated_sinusoids(float *x, float *w, float *ai, float *bi, float *ci, float *di, float *y, int N, int len)
-{
+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){
float *cos_table[N];
float *sin_table[N];
float *tcos_table[N];
float *tsin_table[N];
float cosE[N], sinE[N];
float costE[N], sintE[N];
- int i,j, iter;
-
+ float ai[N];
+ float bi[N];
+ float ci[N];
+ float di[N];
+ int i,j, iter;
+
for(i=0;i<N;i++){
cos_table[i] = malloc(sizeof(**cos_table)*len);
sin_table[i] = malloc(sizeof(**sin_table)*len);
@@ -146,8 +150,8 @@
for (j=0;j<len;j++)
{
float jj = j-len/2.+.5;
- cos_table[i][j] = cos(w[i]*jj);
- sin_table[i][j] = sin(w[i]*jj);
+ cos_table[i][j] = cos((2.*M_PI*w[i]/len)*jj);
+ sin_table[i][j] = sin((2.*M_PI*w[i]/len)*jj);
tcos_table[i][j] = ((jj))*cos_table[i][j];
tsin_table[i][j] = ((jj))*sin_table[i][j];
/* The sinusoidal terms */
@@ -244,12 +248,217 @@
free(tcos_table[i]);
free(tsin_table[i]);
}
+
+ for(i=0;i<N;i++){
+ float A = hypot(ai[i],bi[i]);
+ double P;
+ if(A!=0.){
+ if(bi[i]>0){
+ P = -acos(ai[i]/A);
+ }else{
+ P = acos(ai[i]/A);
+ }
+ }else
+ P=0.;
+
+ Aout[i] = A;
+ Wout[i] = w[i] + ((ci[i]*bi[i] - di[i]*ai[i])/(A*A)/(2.*M_PI))*len;
+ Pout[i] = P;
+ delAout[i] = (ci[i]*ai[i] + bi[i]*di[i])/A;
+ delWout[i] = 0.;
+
+ }
+
+}
+/* Models the signal as modulated sinusoids
+ * x is the signal
+ * w are the frequencies
+ * 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_sinusoidsB(float *x, float *w,
+ float *Aout, float *Wout, float *Pout,
+ float *dAout, float *dWout, float *ddAout,
+ float *y, int N, int len){
+ float *cos_table[N];
+ float *sin_table[N];
+ float *tcos_table[N];
+ float *tsin_table[N];
+ float *ttcos_table[N];
+ float *ttsin_table[N];
+ float cosE[N], sinE[N];
+ float tcosE[N], tsinE[N];
+ float ttcosE[N], ttsinE[N];
+ float ai[N];
+ float bi[N];
+ float ci[N];
+ float di[N];
+ float ei[N];
+ float fi[N];
+ int i,j, iter;
+
+ /* Build a table for the four basis functions at each frequency */
+ for (i=0;i<N;i++){
+ float tmpa=0;
+ float tmpb=0;
+ float tmpc=0;
+ 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;
+
+ float c = cos((2.*M_PI*w[i]/len)*jj);
+ float s = sin((2.*M_PI*w[i]/len)*jj);
+
+ cos_table[i][j] = c;
+ sin_table[i][j] = s;
+ tcos_table[i][j] = jj*c;
+ tsin_table[i][j] = jj*s;
+ ttcos_table[i][j] = jj*jj*c;
+ ttsin_table[i][j] = jj*jj*s;
+
+ /* The sinusoidal terms */
+ tmpa += cos_table[i][j]*cos_table[i][j];
+ tmpb += sin_table[i][j]*sin_table[i][j];
+
+ /* The modulation terms */
+ tmpc += tcos_table[i][j]*tcos_table[i][j];
+ tmpd += tsin_table[i][j]*tsin_table[i][j];
+
+ /* The second order modulations */
+ tmpe += ttcos_table[i][j]*ttcos_table[i][j];
+ tmpf += ttsin_table[i][j]*ttsin_table[i][j];
+
+ }
+
+ tmpa = sqrt(tmpa);
+ tmpb = sqrt(tmpb);
+ tmpc = sqrt(tmpc);
+ 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];
+ tcos_table[i][j] *= tcosE[i];
+ tsin_table[i][j] *= tsinE[i];
+ ttcos_table[i][j] *= ttcosE[i];
+ ttsin_table[i][j] *= ttsinE[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] = ei[i] = fi[i] = 0;
+
+ 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;
+
+ tmpa += (x[j]-y[j]) * cos_table[i][j];
+ tmpb += (x[j]-y[j]) * sin_table[i][j];
+ tmpc += (x[j]-y[j]) * tcos_table[i][j];
+ tmpd += (x[j]-y[j]) * tsin_table[i][j];
+ tmpe += (x[j]-y[j]) * ttcos_table[i][j];
+ tmpf += (x[j]-y[j]) * ttsin_table[i][j];
+
+ }
+
+ /* 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;
+ bi[i] += tmpb;
+ ci[i] += tmpc;
+ di[i] += tmpd;
+ ei[i] += tmpe;
+ fi[i] += tmpf;
+
+ }
+ }
+
+ for (i=0;i<N;i++){
+ ai[i] *= cosE[i];
+ bi[i] *= sinE[i];
+ ci[i] *= tcosE[i];
+ di[i] *= tsinE[i];
+ ei[i] *= ttcosE[i];
+ fi[i] *= ttsinE[i];
+ }
+
+ for(i=0;i<N;i++){
+ float A = hypot(ai[i],bi[i]);
+ double P;
+
+ free(cos_table[i]);
+ free(sin_table[i]);
+ free(tcos_table[i]);
+ 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);
+ }
+ }else
+ P=0.;
+
+ Aout[i] = A;
+ Pout[i] = P;
+ dAout[i] = (ci[i]*ai[i] + bi[i]*di[i])/A;
+ Wout[i] = w[i] + ((ci[i]*bi[i] - di[i]*ai[i])/(A*A)/(2.*M_PI))*len;
+ ddAout[i] = (ei[i]*ai[i] + bi[i]*fi[i])/A;
+ dWout[i] = ((ei[i]*bi[i] - fi[i]*ai[i])/(A*A)/(2.*M_PI))*len;
+ }
+
}
-void compute_sinusoids(float *Aout, float *Wout, float *Pout,
- float *delAout, float *delWout, float *y, int N, int len){
+static void compute_sinusoids(float *A, float *W, float *P,
+ float *delA, float *delW, float *y, int N, int len){
int i,j;
double ilen = 1./len;
@@ -257,50 +466,307 @@
double jj = (j-len*.5+.5)*ilen;
y[j]=0.;
for(i=0;i<N;i++){
- double A = Aout[i] + delAout[i]*jj;
- double W = Wout[i] + delWout[i]*jj;
- y[j] += A*cos(M_2_PI*jj*W+Pout[i]);
+ double a = A[i] + delA[i]*jj;
+ double w = W[i] + delW[i]*jj;
+ y[j] += a*cos(2.*M_PI*jj*w+P[i]);
}
}
}
-// brute-force, simultaneous nonlinear system CG, meant to be a 'high anchor' of possible numeric performance
-void extract_modulated_sinusoids2(float *x, float *w, float *Aout, float *Wout, float *Pout,
- float *delAout, float *delWout, float *y, int N, int len){
- int i,j;
- double ilen = 1./len;
+/* 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
+*/
- /* w contains our frequency seeds; initialize Wout, Pout, Aout from this data */
- /* delAout and delPout start at zero */
+void extract_modulated_sinusoidsLSQR(float *x, float *w,
+ float *Aout, float *Wout, float *Pout, float *dAout, float *dWout, float *ddAout,
+ float *y, int N, int len){
+ float *cos_table[N];
+ float *sin_table[N];
+ float *tcos_table[N];
+ float *tsin_table[N];
+ float cosE[N], sinE[N];
+ float costE[N], sintE[N];
+ float ai[N];
+ float bi[N];
+ float ci[N];
+ float di[N];
+ int i,j, iter;
+
+ float e[len];
+ /* Symmetric and anti-symmetric components of the error */
+ int L2 = len/2;
+ float sym[L2], anti[L2];
+
for(i=0;i<N;i++){
- double re = 0.;
- double im = 0.;
+ cos_table[i] = malloc(sizeof(**cos_table)*len);
+ sin_table[i] = malloc(sizeof(**sin_table)*len);
+ tcos_table[i] = malloc(sizeof(**tcos_table)*len);
+ tsin_table[i] = malloc(sizeof(**tsin_table)*len);
+ }
- for(j=0;j<len;j++){
- double jj = M_2_PI*(j-len*.5+.5)*ilen*w[i];
- double s = sin(jj);
- double c = cos(jj);
- im += s*x[j];
- re += c*x[j];
+ /* 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;
+ float rotR, rotI;
+ rotR = cos(2*M_PI*w[i]/len);
+ rotI = sin(2*M_PI*w[i]/len);
+ /* Computing sin/cos table using complex rotations */
+ cos_table[i][0] = cos(M_PI*w[i]/len);
+ sin_table[i][0] = sin(M_PI*w[i]/len);
+ for (j=1;j<L2;j++){
+ float re, im;
+ re = cos_table[i][j-1]*rotR - sin_table[i][j-1]*rotI;
+ im = sin_table[i][j-1]*rotR + cos_table[i][j-1]*rotI;
+ cos_table[i][j] = re;
+ sin_table[i][j] = im;
}
+ /* Only need to compute the tables for half the length because of the symmetry.
+ Eventually, we'll have to replace the cos/sin with rotations */
+ for (j=0;j<L2;j++){
+ float jj = j+.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];
+ }
+ /* Double the energy because we only computed one half.
+ Eventually, we should be computing/tabulating these values directly
+ as a function of w[i]. */
+ cosE[i] = sqrt(2*tmp1);
+ sinE[i] = sqrt(2*tmp2);
+ costE[i] = sqrt(2*tmp3);
+ sintE[i] = sqrt(2*tmp4);
+ /* Normalise the basis (should multiply by the inverse instead) */
+ for (j=0;j<L2;j++){
+ cos_table[i][j] *= (1.f/cosE[i]);
+ sin_table[i][j] *= (1.f/sinE[i]);
+ tcos_table[i][j] *= (1.f/costE[i]);
+ tsin_table[i][j] *= (1.f/sintE[i]);
+ }
+ }
+ /* y is the initial approximation of the signal */
+ for (j=0;j<len;j++)
+ y[j] = 0;
+ for (j=0;j<len;j++)
+ e[j] = x[j];
+ /* Split the error into a symmetric component and an anti-symmetric component.
+ This speeds everything up by a factor of 2 */
+ for (j=0;j<L2;j++){
+ sym[j] = e[j+L2]+e[L2-j-1];
+ anti[j] = e[j+L2]-e[L2-j-1];
+ }
+
+ for (i=0;i<N;i++)
+ ai[i] = bi[i] = ci[i] = di[i] = 0;
+
+ /* y is the initial approximation of the signal */
+ for (j=0;j<len;j++)
+ y[j] = 0;
+ for (j=0;j<len;j++)
+ e[j] = x[j];
+ /* Split the error into a symmetric component and an anti-symmetric component.
+ This speeds everything up by a factor of 2 */
+ for (j=0;j<L2;j++){
+ sym[j] = e[j+L2]+e[L2-j-1];
+ anti[j] = e[j+L2]-e[L2-j-1];
+ }
+
+ for (i=0;i<N;i++)
+ ai[i] = bi[i] = ci[i] = di[i] = 0;{
+ float va[N];
+ float vb[N];
+ float vc[N];
+ float vd[N];
+ float wa[N];
+ float wb[N];
+ float wc[N];
+ float wd[N];
+ float alpha;
+ float beta;
+ float theta;
+ float phi;
+ float phibar;
+ float rho;
+ float rhobar;
+ float c;
+ float s;
+ /*beta*u=e*/
+ beta = phibar = 0;
+ for (j=0;j<L2;j++)
+ beta += sym[j]*sym[j] + anti[j]*anti[j];
+ phibar = beta = sqrtf(0.5f*beta);
+ if (beta > 0){
+ for (j=0;j<L2;j++){
+ sym[j] *= (1.f/beta);
+ anti[j] *= (1.f/beta);
+ }
+ /*alpha*v=A^T.u*/
+ alpha = 0;
+ for (i=0;i<N;i++) {
+ va[i] = 0;
+ for (j=0;j<L2;j++)
+ va[i] += sym[j]*cos_table[i][j];
+ vb[i] = 0;
+ for (j=0;j<L2;j++)
+ vb[i] += anti[j]*sin_table[i][j];
+ vc[i] = 0;
+ for (j=0;j<L2;j++)
+ vc[i] += anti[j]*tcos_table[i][j];
+ vd[i] = 0;
+ for (j=0;j<L2;j++)
+ vd[i] += sym[j]*tsin_table[i][j];
+ alpha += va[i]*va[i]+vb[i]*vb[i]+vc[i]*vc[i]+vd[i]*vd[i];
+ }
+ if (alpha > 0){
+ alpha = sqrtf(alpha);
+ for (i=0;i<N;i++){
+ wa[i] = va[i] *= (1.f/alpha);
+ wb[i] = vb[i] *= (1.f/alpha);
+ wc[i] = vc[i] *= (1.f/alpha);
+ wd[i] = vd[i] *= (1.f/alpha);
+ }
+ phibar = beta;
+ rhobar = alpha;
+ for (iter=0;iter<5;iter++){
+ float p;
+ float q;
+ /*beta*u=A.v-alpha*u*/
+ beta = 0;
+ for (j=0;j<L2;j++){
+ float da;
+ float ds;
+ ds = da = 0;
+ for (i=0;i<N;i++){
+ ds += va[i]*cos_table[i][j] + vd[i]*tsin_table[i][j];
+ da += vb[i]*sin_table[i][j] + vc[i]*tcos_table[i][j];
+ }
+ sym[j] = 2*ds-alpha*sym[j];
+ anti[j] = 2*da-alpha*anti[j];
+ beta += sym[j]*sym[j] + anti[j]*anti[j];
+ }
+ if (beta <= 0) break;
+ beta = sqrtf(0.5f*beta);
+ for (j=0;j<L2;j++){
+ sym[j] *= (1.f/beta);
+ anti[j] *= (1.f/beta);
+ }
+ /*alpha*v=A^T.u-beta*v*/
+ alpha = 0;
+ for (i=0;i<N;i++){
+ float v0;
+ float v1;
+ float v2;
+ float v3;
+ v0 = 0;
+ for (j=0;j<L2;j++)
+ v0 += sym[j]*cos_table[i][j];
+ v1 = 0;
+ for (j=0;j<L2;j++)
+ v1 += anti[j]*sin_table[i][j];
+ v2 = 0;
+ for (j=0;j<L2;j++)
+ v2 += anti[j]*tcos_table[i][j];
+ v3 = 0;
+ for (j=0;j<L2;j++)
+ v3 += sym[j]*tsin_table[i][j];
+ va[i] = v0 - beta*va[i];
+ vb[i] = v1 - beta*vb[i];
+ vc[i] = v2 - beta*vc[i];
+ vd[i] = v3 - beta*vd[i];
+ alpha += va[i]*va[i]+vb[i]*vb[i]+vc[i]*vc[i]+vd[i]*vd[i];
+ }
+ if (alpha <= 0) break;
+ alpha = sqrtf(alpha);
+ for (i=0;i<N;i++){
+ va[i] *= (1.f/alpha);
+ vb[i] *= (1.f/alpha);
+ vc[i] *= (1.f/alpha);
+ vd[i] *= (1.f/alpha);
+ }
+ rho = hypotf(rhobar,beta);
+ c = rhobar/rho;
+ s = beta/rho;
+ theta = s*alpha;
+ rhobar = -c*alpha;
+ phi = c*phibar;
+ phibar = s*phibar;
+ p = phi/rho;
+ q = theta/rho;
+ for (i=0;i<N;i++){
+ ai[i] += wa[i]*p;
+ bi[i] += wb[i]*p;
+ ci[i] += wc[i]*p;
+ di[i] += wd[i]*p;
+ wa[i] = va[i] - wa[i]*q;
+ wb[i] = vb[i] - wb[i]*q;
+ wc[i] = vc[i] - wc[i]*q;
+ wd[i] = vd[i] - wd[i]*q;
+ }
+ }
+ }
+ }
+ fprintf(stderr,"LSQR sq. err. = %lg\n",phibar*phibar);
+ phibar *= 0.5f;
+ for (j=0;j<L2;j++){
+ e[j+L2] = phibar*(sym[j]+anti[j]);
+ e[L2-j-1] = phibar*(sym[j]-anti[j]);
+ }
+ }
+ for (j=0;j<L2;j++){
+ float da;
+ float ds;
+ ds = da = 0;
+ for (i=0;i<N;i++){
+ ds += ai[i]*cos_table[i][j] + di[i]*tsin_table[i][j];
+ da += bi[i]*sin_table[i][j] + ci[i]*tcos_table[i][j];
+ }
+ y[j+L2] = ds+da;
+ y[L2-j-1] = ds-da;
+ }
+ for (i=0;i<N;i++){
+ ai[i] /= cosE[i];
+ bi[i] /= sinE[i];
+ ci[i] /= costE[i];
+ di[i] /= sintE[i];
+ }
- double ph = 0;
- double m = hypot(re,im);
- if(m!=0.){
- if(im>0){
- ph = acos(re/m)/M_PI;
+ for(i=0;i<N;i++){
+ float A = hypot(ai[i],bi[i]);
+ double P;
+
+ free(cos_table[i]);
+ free(sin_table[i]);
+
+ if(A!=0.){
+ if(bi[i]>0){
+ P = -acos(ai[i]/A);
}else{
- ph = -acos(re/m)/M_PI;
+ P = acos(ai[i]/A);
}
- }
-
- Aout[i] = 2*m*ilen;
- Pout[i] = ph;
- Wout[i] = w[i];
- delAout[i] = 0.f;
- delWout[i] = 0.f;
+ }else
+ P=0.;
+
+ Aout[i] = A;
+ Pout[i] = P;
+ dAout[i] = (ci[i]*ai[i] + bi[i]*di[i])/A;
+ Wout[i] = w[i] + ((ci[i]*bi[i] - di[i]*ai[i])/(A*A)/(2.*M_PI))*len;
+
}
- compute_sinusoids(Aout,Wout,Pout,delAout,delWout,y,N,len);
-
}
Modified: trunk/ghost/monty/sinusoid/sinusoids.h
===================================================================
--- trunk/ghost/monty/sinusoid/sinusoids.h 2007-07-27 20:57:57 UTC (rev 13374)
+++ trunk/ghost/monty/sinusoid/sinusoids.h 2007-07-28 01:39:06 UTC (rev 13375)
@@ -10,7 +10,7 @@
* *
********************************************************************
- function: research-grade sinusoidal extraction sode
+ function: research-grade sinusoidal extraction code
last mod: $Id$
********************************************************************/
@@ -18,7 +18,14 @@
extern void window_weight(float *logf, float *out, int n, float flatbias,
int lowindow, int hiwindow, int min, int rate);
-extern void extract_modulated_sinusoids(float *x, float *w, float *ai, float *bi, float *ci, float *di, float *y, int N, int len);
-
+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);
+extern void extract_modulated_sinusoidsLSQR(float *x, float *w, float *Aout, float *Wout, float *Pout,
+ float *delAout, float *delWout, float *ddAout, float *y, int N, int len);
+
extern void extract_modulated_sinusoids2(float *x, float *w, float *Aout, float *Wout, float *Pout,
float *delAout, float *delWout, float *y, int N, int len);
+extern void extract_modulated_sinusoidsB(float *x, float *w,
+ float *Aout, float *Wout, float *Pout,
+ float *dAout, float *dWout, float *ddAout,
+ float *y, int N, int len);
Modified: trunk/ghost/monty/sinusoid/test_sinusoids.c
===================================================================
--- trunk/ghost/monty/sinusoid/test_sinusoids.c 2007-07-27 20:57:57 UTC (rev 13374)
+++ trunk/ghost/monty/sinusoid/test_sinusoids.c 2007-07-28 01:39:06 UTC (rev 13375)
@@ -85,11 +85,11 @@
fout = fopen(argv[2], "w");
/*
- for(i=0;i<BLOCK_SIZE*2;i++)
+ 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));
+ 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){
@@ -109,90 +109,100 @@
float log_fft[BLOCK_SIZE+1];
float weight[BLOCK_SIZE+1];
- memcpy(fft_buf,float_in,sizeof(float_in));
-
// polish the strongest peaks from weighting
if(frame==220){
float w[BLOCK_SIZE];
float Aout[BLOCK_SIZE]={0};
+ float Pout[BLOCK_SIZE]={0};
+ float dAout[BLOCK_SIZE]={0};
float Wout[BLOCK_SIZE]={0};
- float Pout[BLOCK_SIZE]={0};
- float delAout[BLOCK_SIZE]={0};
- float delWout[BLOCK_SIZE]={0};
+ float ddAout[BLOCK_SIZE]={0};
+ float dWout[BLOCK_SIZE]={0};
float y[BLOCK_SIZE*2];
- //for(j=0;j<50;j++){
- j=BLOCK_SIZE;
- blackmann_harris(fft_buf, fft_buf, BLOCK_SIZE*2);
- //if(j==0)
- 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] *= 2./BLOCK_SIZE;
+ 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;
- mag_dB(log_fft,fft_buf,BLOCK_SIZE*2);
- //if(j==0)
- dump_vec(log_fft,BLOCK_SIZE+1,"logmag",frame);
+ 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);
- window_weight(log_fft,weight,BLOCK_SIZE+1, 2.f, 512,256, 30, 44100);
+ j=1;
+ w[0]=.0306*BLOCK_SIZE;
+
+ //for(j=0;j<BLOCK_SIZE;j++)
+ //w[j] = j+.5;
- // largest weighted
- /*int best=-120;
+ /* largest weighted
+ int best=-120;
int besti=-1;
- for(i=0;i<BLOCK_SIZE+1;i++)
- if(weight[i]>best){
- besti = i;
- best = weight[i];
+ for(i=0;i<BLOCK_SIZE+1;i++){
+ float v = log_fft[i] - weight[i];
+ if(v>best){
+ besti = i;
+ best = v;
+ }
}
- w[j] = M_PI*besti/BLOCK_SIZE;
+ //w[j] = M_PI*besti/BLOCK_SIZE;
*/
- for(i=0;i<j;i++)
- w[i] = i;//M_PI*(i)/j;
+ //blackmann_harris(fft_buf, float_in, BLOCK_SIZE*2);
+ extract_modulated_sinusoidsB(float_in, w, Aout, Wout, Pout, dAout, dWout, ddAout, y, j, BLOCK_SIZE*2);
+ w[0]=Wout[0];
+ extract_modulated_sinusoidsB(float_in, w, Aout, Wout, Pout, dAout, dWout, ddAout, y, j, BLOCK_SIZE*2);
+ w[0]=Wout[0];
+ extract_modulated_sinusoidsB(float_in, w, Aout, Wout, Pout, dAout, dWout, ddAout, y, j, BLOCK_SIZE*2);
+ w[0]=Wout[0];
+ extract_modulated_sinusoidsB(float_in, w, Aout, Wout, Pout, dAout, dWout, ddAout, y, j, BLOCK_SIZE*2);
+ w[0]=Wout[0];
+ extract_modulated_sinusoidsB(float_in, w, Aout, Wout, Pout, dAout, dWout, ddAout, y, j, BLOCK_SIZE*2);
+ w[0]=Wout[0];
+ extract_modulated_sinusoidsB(float_in, w, Aout, Wout, Pout, dAout, dWout, ddAout, y, j, BLOCK_SIZE*2);
- blackmann_harris(fft_buf, float_in, BLOCK_SIZE*2);
- dump_vec(fft_buf,BLOCK_SIZE*2,"win",frame);
- extract_modulated_sinusoids2(fft_buf, w, Aout, Wout, Pout, delAout, delWout, y, j, BLOCK_SIZE*2);
+ //for(i=0;i<j;i++)
+ //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);
- //for(i=0;i<j;i++){
- //float A = sqrt(a[i]*a[i] + b[i]*b[i]);
- //w[i] += (c[i]*b[i] - d[i]*a[i])/(A*A);
- //}
-
- memcpy(fft_buf,float_in,sizeof(float_in));
- for(i=0;i<BLOCK_SIZE*2;i++)
- fft_buf[i] -= y[i];
-
- //dump_vec(float_in,BLOCK_SIZE*2,"exdata",frame);
- dump_vec(y,BLOCK_SIZE*2,"extract",frame);
- //}
-
- for(i=0;i<j;i++){
- //float A = sqrt(a[i]*a[i] + b[i]*b[i]);
- Aout[i]=todB(Aout[i]);
- Wout[i] /= BLOCK_SIZE;
- }
-
- dump_vec2(Wout,Aout,j,"ex",frame);
-
- //dump_vec(float_in,BLOCK_SIZE*2,"exdata",frame);
- dump_vec(y,BLOCK_SIZE*2,"extract",frame);
-
blackmann_harris(fft_buf, fft_buf, BLOCK_SIZE*2);
drft_forward(&fft, fft_buf);
- for(i=0;i<BLOCK_SIZE*2;i++)fft_buf[i] *= 2./BLOCK_SIZE;
+ 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",frame);
+ 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]);
+ ddAout[i]=todB(ddAout[i]);
+ w[i] /= BLOCK_SIZE;
+ }
+
+ dump_vec2(w,Aout,j,"ex",0);
+ dump_vec2(w,dAout,j,"dA",0);
+ dump_vec2(w,ddAout,j,"ddA",0);
+ dump_vec2(w,dWout,j,"dW",0);
+ dump_vec(y,BLOCK_SIZE*2,"extract",0);
+
+
}
+
}
-
-
+
+
+
//fwrite(short_in, sizeof(short), BLOCK_SIZE, fout);
frame++;
- }
-
- return 0;
+ }
+
+ return 0;
}
More information about the commits
mailing list