[xiph-commits] r13420 - trunk/ghost/libghost
jm at svn.xiph.org
jm at svn.xiph.org
Tue Jul 31 23:24:12 PDT 2007
Author: jm
Date: 2007-07-31 23:24:12 -0700 (Tue, 31 Jul 2007)
New Revision: 13420
Modified:
trunk/ghost/libghost/sinusoids.c
trunk/ghost/libghost/sinusoids.h
Log:
new constrained MP algo
Modified: trunk/ghost/libghost/sinusoids.c
===================================================================
--- trunk/ghost/libghost/sinusoids.c 2007-08-01 06:24:05 UTC (rev 13419)
+++ trunk/ghost/libghost/sinusoids.c 2007-08-01 06:24:12 UTC (rev 13420)
@@ -98,10 +98,10 @@
//printf ("\n");
}
-#define MP_RES 4096
+#define MP_RES 8192
#define MP_LEN 256
/* This is a trivial matching pursuit implementation. It's awfully slow, so beware */
-void extract_sinusoids_mp(float *x, float *w, float *y, int N, int len)
+void extract_sinusoids_mp(float *x, float *w, float *window, float *y, int N, int len)
{
static float cos_table[MP_RES][MP_LEN];
static float sin_table[MP_RES][MP_LEN];
@@ -119,8 +119,8 @@
for (j=0;j<MP_LEN;j++)
{
float jj=j-(MP_LEN/2)+.5;
- cos_table[i][j] = cos((M_PI*(float)i/(float)MP_RES)*jj);
- sin_table[i][j] = sin((M_PI*(float)i/(float)MP_RES)*jj);
+ cos_table[i][j] = window[j]*cos((M_PI*(float)i/(float)MP_RES)*jj);
+ sin_table[i][j] = window[j]*sin((M_PI*(float)i/(float)MP_RES)*jj);
cosE[i] += cos_table[i][j]*cos_table[i][j];
sinE[i] += sin_table[i][j]*sin_table[i][j];
}
@@ -171,6 +171,84 @@
}
+/* This is a trivial matching pursuit implementation. It's awfully slow, so beware */
+void extract_sinusoids_mp_constrained(float *x, float *w, float *window, float *y, int N, int len)
+{
+ static float cos_table[MP_RES][MP_LEN];
+ static float sin_table[MP_RES][MP_LEN];
+ static float sinE[MP_RES];
+ static float cosE[MP_RES];
+ static init = 0;
+ float res[len];
+ int i,j;
+ if (!init)
+ {
+ for (i=0;i<MP_RES;i++)
+ {
+ sinE[i] = 0;
+ cosE[i] = 0;
+ for (j=0;j<MP_LEN;j++)
+ {
+ float jj=j-(MP_LEN/2)+.5;
+ cos_table[i][j] = window[j]*cos((M_PI*(float)i/(float)MP_RES)*jj);
+ sin_table[i][j] = window[j]*sin((M_PI*(float)i/(float)MP_RES)*jj);
+ cosE[i] += cos_table[i][j]*cos_table[i][j];
+ sinE[i] += sin_table[i][j]*sin_table[i][j];
+ }
+ }
+ init = 1;
+ }
+
+ for (i=0;i<len;i++)
+ res[i] = x[i];
+ //Loop on all sinusoids we want to extract
+ for (i=0;i<N;i++)
+ {
+ int k;
+ float best_fit=-1;
+ int best_id=1;
+ float best_gainC=0;
+ float best_gainS=0;
+ float fit = 0;
+ int start, end;
+
+ start = floor(MP_RES*(w[i]-.025f)/M_PI);
+ end = ceil(MP_RES*(w[i]+.025f)/M_PI);
+ if (start<1)
+ start = 1;
+ if (end > MP_RES-1)
+ end = MP_RES-1;
+ //Find best sinusoid
+ for (j=start;j<end;j++)
+ {
+ float sumC = 0, sumS = 0;
+ for (k=0;k<len;k++)
+ {
+ sumC += res[k]*cos_table[j][k];
+ sumS += res[k]*sin_table[j][k];
+ }
+ fit = sumC/cosE[j]*sumC + sumS/sinE[j]*sumS;
+ if (fit > best_fit)
+ {
+ best_id = j;
+ best_fit = fit;
+ best_gainC = sumC/cosE[j];
+ best_gainS = sumS/sinE[j];
+ }
+ }
+
+ w[i] = M_PI*best_id/MP_RES;
+
+ //Remove the sinusoid we found
+ for (k=0;k<len;k++)
+ res[k] = res[k] - best_gainC*cos_table[best_id][k] - best_gainS*sin_table[best_id][k];
+ //printf ("%d %f %f\n", best_id, best_gainC, best_gainS);
+ }
+ for (i=0;i<len;i++)
+ y[i] = x[i] - res[i];
+
+}
+
void extract_sinusoids(float *x, float *w, float *window, float *ai, float *bi, float *y, int N, int len)
{
float cos_table[N][len];
@@ -260,7 +338,7 @@
for (i=0;i<N;i++)
ai[i] = bi[i] = ci[i] = di[i] = 0;
- nonlin = 10;
+ nonlin = 3;
/* Iterate on the non-linear part (frequency) a couple times */
while (1)
{
@@ -345,7 +423,9 @@
for (j=0;j<L2;j++)
anti[j] -= (2*bi[i])*sin_table[i][j] + (2*ci[i])*tcos_table[i][j];
}
-
+ //for (i=0;i<N;i++)
+ // printf ("%g ", w[i]);
+
/* Stop condition. We stop here, just after updating the error using the new frequency */
if (nonlin-- <= 0)
break;
@@ -422,7 +502,7 @@
// printf("%d %f %f %f %f %f %f %f %f %f %f\n", count, w[i], ai[i], bi[i], ci[i], di[i], A[i], phi[i], dA[i], dw[i], w0[i]);
}
}
-
+ //printf ("\n");
for (j=0;j<L2;j++)
{
e[j+L2] = .5*(sym[j]+anti[j]);
Modified: trunk/ghost/libghost/sinusoids.h
===================================================================
--- trunk/ghost/libghost/sinusoids.h 2007-08-01 06:24:05 UTC (rev 13419)
+++ trunk/ghost/libghost/sinusoids.h 2007-08-01 06:24:12 UTC (rev 13420)
@@ -24,7 +24,7 @@
#define _SINUSOIDS_H
void find_sinusoids(float *psd, float *w, int *N, int length);
-
+void extract_sinusoids_mp(float *x, float *w, float *window, float *y, int N, int len);
void extract_sinusoids(float *x, float *w, float *window, float *ai, float *bi, float *y, int N, int len);
void extract_modulated_sinusoids(float *x, float *w, float *window, float *ai, float *bi, float *ci, float *di, float *y, int N, int len);
More information about the commits
mailing list