[xiph-commits] r13248 - trunk/ghost/monty/sinusoid
xiphmont at svn.xiph.org
xiphmont at svn.xiph.org
Thu Jul 12 09:34:58 PDT 2007
Author: xiphmont
Date: 2007-07-12 09:34:57 -0700 (Thu, 12 Jul 2007)
New Revision: 13248
Modified:
trunk/ghost/monty/sinusoid/sinusoids.c
trunk/ghost/monty/sinusoid/sinusoids.h
trunk/ghost/monty/sinusoid/test_sinusoids.c
Log:
Nothing to see here
Modified: trunk/ghost/monty/sinusoid/sinusoids.c
===================================================================
--- trunk/ghost/monty/sinusoid/sinusoids.c 2007-07-12 14:53:21 UTC (rev 13247)
+++ trunk/ghost/monty/sinusoid/sinusoids.c 2007-07-12 16:34:57 UTC (rev 13248)
@@ -22,8 +22,8 @@
#include <stdlib.h>
#include <string.h>
-void level_mean(float *f, float *out, int n,
- int lowindow, int hiwindow, int min, int rate){
+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;
@@ -31,8 +31,11 @@
for(i=0;i<n;i++)
bark[i] = rint(toBark(i*binwidth) * 1024);
- float nacc = 0.f;
- float nadd = 0.f;
+ float gacc = 0.f;
+ float gadd = 0.f;
+ float aacc = 0.f;
+ float aadd = 0.f;
+
float dacc = 0.f;
float dadd = 0.f;
@@ -47,17 +50,20 @@
int c = hihead-i+1;
float d = (c<min?1./min:1./c);
- nadd += f[hihead]*d;
+ gadd += logf[hihead]*d;
+ aadd += fromdB(logf[hihead])*d;
dadd += d;
while(c<min){
- nacc += f[hihead]*d;
+ gacc += logf[hihead]*d;
+ aacc += fromdB(logf[hihead])*d;
dacc += d;
c++;
}
}
- nacc += nadd;
+ aacc += aadd;
+ gacc += gadd;
dacc += dadd;
if(lohead<n){
@@ -69,14 +75,16 @@
{
int c = lohead-i;
float d = 1./c;
- nadd -= f[i]*d;
+ gadd -= logf[i]*d;
+ aadd -= fromdB(logf[i])*d;
dadd -= d;
}
for( ;lotail<i && bark[lotail]<bark[i]-lowindow && lotail<=i-min; lotail++){
int c = i-lotail;
float d = 1./c;
- nadd += f[lotail]*d;
+ gadd += logf[lotail]*d;
+ aadd += fromdB(logf[lotail])*d;
dadd += d;
}
@@ -86,69 +94,21 @@
{
int c = i-hitail+1;
float d = (c<min?1./min:1./c);
- nadd -= f[i]*d;
+ gadd -= logf[i]*d;
+ aadd -= fromdB(logf[i])*d;
dadd -= d;
}
- out[i] = nacc / dacc;
+ {
+ float arith = todB(aacc / dacc);
+ float geom = gacc / dacc;
+ out[i] = logf[i] - arith + flatbias*(arith-geom);
+ }
+
}
}
-
-
-void extract_sinusoids(float *x, float *w, float *window, float *ai, float *bi, float *y, int N, int len)
-{
- float cos_table[N][len];
- float sin_table[N][len];
- float cosE[N], sinE[N];
- int i,j, iter;
- for (i=0;i<N;i++)
- {
- float tmp1=0, tmp2=0;
- for (j=0;j<len;j++)
- {
- cos_table[i][j] = cos(w[i]*j)*window[j];
- sin_table[i][j] = sin(w[i]*j)*window[j];
- tmp1 += cos_table[i][j]*cos_table[i][j];
- tmp2 += sin_table[i][j]*sin_table[i][j];
- }
- cosE[i] = tmp1;
- sinE[i] = tmp2;
- }
- for (j=0;j<len;j++)
- y[j] = 0;
- for (i=0;i<N;i++)
- ai[i] = bi[i] = 0;
-
- for (iter=0;iter<5;iter++)
- {
- for (i=0;i<N;i++)
- {
- float tmp1=0, tmp2=0;
- for (j=0;j<len;j++)
- {
- tmp1 += (x[j]-y[j])*cos_table[i][j];
- tmp2 += (x[j]-y[j])*sin_table[i][j];
- }
- tmp1 /= cosE[i];
- //Just in case it's a DC! Must fix that anyway
- tmp2 /= (.0001+sinE[i]);
- for (j=0;j<len;j++)
- {
- y[j] += tmp1*cos_table[i][j] + tmp2*sin_table[i][j];
- }
- ai[i] += tmp1;
- bi[i] += tmp2;
- if (iter==4)
- {
- printf ("%f %f ", w[i], (float)sqrt(ai[i]*ai[i] + bi[i]*bi[i]));
- }
- }
- }
- printf ("\n");
-}
-
/* Models the signal as modulated sinusoids
* x is the signal
* w are the frequencies
@@ -161,15 +121,23 @@
* N is the number of sinusoids
* len is the frame size
*/
-void extract_modulated_sinusoids(float *x, float *w, float *window, float *ai, float *bi, float *ci, float *di, float *y, int N, int len)
+void extract_modulated_sinusoids(float *x, float *w, float *ai, float *bi, float *ci, float *di, float *y, int N, int len)
{
- float cos_table[N][len];
- float sin_table[N][len];
- float tcos_table[N][len];
- float tsin_table[N][len];
- float cosE[N], sinE[N];
- float costE[N], sintE[N];
+ 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;
+
+ for(i=0;i<N;i++){
+ 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);
+ }
+
/* 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++)
{
@@ -177,9 +145,9 @@
float tmp3=0, tmp4=0;
for (j=0;j<len;j++)
{
- float jj = j-len/2+.5;
- cos_table[i][j] = cos(w[i]*jj)*window[j];
- sin_table[i][j] = sin(w[i]*jj)*window[j];
+ float jj = j-len/2.+.5;
+ cos_table[i][j] = cos(w[i]*jj);
+ sin_table[i][j] = sin(w[i]*jj);
tcos_table[i][j] = ((jj))*cos_table[i][j];
tsin_table[i][j] = ((jj))*sin_table[i][j];
/* The sinusoidal terms */
@@ -269,4 +237,70 @@
//printf ("0 0 0 0 0\n");
//printf ("\n");
+
+ for(i=0;i<N;i++){
+ free(cos_table[i]);
+ free(sin_table[i]);
+ free(tcos_table[i]);
+ free(tsin_table[i]);
+ }
+
+
}
+
+void compute_sinusoids(float *Aout, float *Wout, float *Pout,
+ float *delAout, float *delWout, float *y, int N, int len){
+ int i,j;
+ double ilen = 1./len;
+
+ for(j=0;j<len;j++){
+ 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]);
+ }
+ }
+}
+
+// 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;
+
+ /* w contains our frequency seeds; initialize Wout, Pout, Aout from this data */
+ /* delAout and delPout start at zero */
+ for(i=0;i<N;i++){
+ double re = 0.;
+ double im = 0.;
+
+ 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];
+ }
+
+ double ph = 0;
+ double m = hypot(re,im);
+ if(m!=0.){
+ if(im>0){
+ ph = acos(re/m)/M_PI;
+ }else{
+ ph = -acos(re/m)/M_PI;
+ }
+ }
+
+ Aout[i] = 2*m*ilen;
+ Pout[i] = ph;
+ Wout[i] = w[i];
+ delAout[i] = 0.f;
+ delWout[i] = 0.f;
+ }
+
+ 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-12 14:53:21 UTC (rev 13247)
+++ trunk/ghost/monty/sinusoid/sinusoids.h 2007-07-12 16:34:57 UTC (rev 13248)
@@ -15,4 +15,10 @@
********************************************************************/
-extern void level_mean(float *f,float *out,int n, int lowindow, int hiwindow, int min, int rate);
+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_sinusoids2(float *x, float *w, float *Aout, float *Wout, float *Pout,
+ float *delAout, float *delWout, float *y, int N, int len);
Modified: trunk/ghost/monty/sinusoid/test_sinusoids.c
===================================================================
--- trunk/ghost/monty/sinusoid/test_sinusoids.c 2007-07-12 14:53:21 UTC (rev 13247)
+++ trunk/ghost/monty/sinusoid/test_sinusoids.c 2007-07-12 16:34:57 UTC (rev 13248)
@@ -27,6 +27,24 @@
}
+static void dump_vec2(float *x, float *y, int n, char *base, int i){
+ char *filename;
+ FILE *f;
+
+ asprintf(&filename,"%s_%d.m",base,i);
+ f=fopen(filename,"w");
+
+ for(i=0;i<n;i++){
+ fprintf(f,"%f -120\n",
+ x[i]);
+ fprintf(f,"%f %f\n\n",
+ x[i],y[i]);
+ }
+
+ fclose(f);
+
+}
+
#define A0 .35875f
#define A1 .48829f
#define A2 .14128f
@@ -43,12 +61,12 @@
}
}
-void mag_dB(float *d, int n){
+void mag_dB(float *log,float *d, int n){
int i;
- d[0] = todB(d[0]*d[0])*.5;
+ log[0] = todB(d[0]*d[0])*.5;
for(i=2;i<n;i+=2)
- d[i>>1] = todB(d[i-1]*d[i-1] + d[i]*d[i])*.5;
- d[n/2] = todB(d[n-1]*d[n-1])*.5;
+ log[i>>1] = todB(d[i-1]*d[i-1] + d[i]*d[i])*.5;
+ log[n/2] = todB(d[n-1]*d[n-1])*.5;
}
int main(int argc, char **argv){
@@ -75,7 +93,7 @@
*/
while (1){
- int i;
+ int i,j;
memmove(float_in,float_in+BLOCK_SIZE,sizeof(*float_in)*BLOCK_SIZE);
fread(short_in, sizeof(short), BLOCK_SIZE, fin);
@@ -85,23 +103,89 @@
float_in[i+BLOCK_SIZE] = short_in[i] * .000030517578125;
{
+
/* generate a log spectrum */
float fft_buf[BLOCK_SIZE*2];
+ float log_fft[BLOCK_SIZE+1];
float weight[BLOCK_SIZE+1];
- blackmann_harris(fft_buf, float_in, BLOCK_SIZE*2);
- //dump_vec(float_in,BLOCK_SIZE*2,"data",frame);
- dump_vec(fft_buf,BLOCK_SIZE*2,"win_data",frame);
- drft_forward(&fft, fft_buf);
- for(i=0;i<BLOCK_SIZE*2;i++)fft_buf[i] *= 2./BLOCK_SIZE;
- //dump_vec(fft_buf,BLOCK_SIZE*2,"fft",frame);
+ memcpy(fft_buf,float_in,sizeof(float_in));
- mag_dB(fft_buf,BLOCK_SIZE*2);
- dump_vec(fft_buf,BLOCK_SIZE+1,"logmag",frame);
+ // polish the strongest peaks from weighting
+ if(frame==220){
+ float w[BLOCK_SIZE];
+ float Aout[BLOCK_SIZE]={0};
+ float Wout[BLOCK_SIZE]={0};
+ float Pout[BLOCK_SIZE]={0};
+ float delAout[BLOCK_SIZE]={0};
+ float delWout[BLOCK_SIZE]={0};
+ float y[BLOCK_SIZE*2];
- level_mean(fft_buf,weight,BLOCK_SIZE+1, 512,1024, 15, 44100);
- dump_vec(weight,BLOCK_SIZE+1,"weight",frame);
+ //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;
+ mag_dB(log_fft,fft_buf,BLOCK_SIZE*2);
+ //if(j==0)
+ dump_vec(log_fft,BLOCK_SIZE+1,"logmag",frame);
+
+ window_weight(log_fft,weight,BLOCK_SIZE+1, 2.f, 512,256, 30, 44100);
+
+ // 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];
+ }
+ 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);
+ 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++){
+ //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;
+ mag_dB(log_fft,fft_buf,BLOCK_SIZE*2);
+ dump_vec(log_fft,BLOCK_SIZE+1,"exlogmag",frame);
+
+
+ }
}
More information about the commits
mailing list