[xiph-commits] r17920 - trunk/ghost/monty/chirp
xiphmont at svn.xiph.org
xiphmont at svn.xiph.org
Mon Apr 11 23:45:19 PDT 2011
Author: xiphmont
Date: 2011-04-11 23:45:19 -0700 (Mon, 11 Apr 2011)
New Revision: 17920
Added:
trunk/ghost/monty/chirp/chirpgraph.c
Modified:
trunk/ghost/monty/chirp/chirp.c
trunk/ghost/monty/chirp/chirp.h
trunk/ghost/monty/chirp/harness.c
Log:
Get a few things in SVN to avoid loss. None should be inspected, all is work in progress.
Modified: trunk/ghost/monty/chirp/chirp.c
===================================================================
--- trunk/ghost/monty/chirp/chirp.c 2011-03-28 18:22:35 UTC (rev 17919)
+++ trunk/ghost/monty/chirp/chirp.c 2011-04-12 06:45:19 UTC (rev 17920)
@@ -10,9 +10,17 @@
* *
********************************************************************
- function: research-grade chirp extraction code
+ function: research-grade chirp estimation code
last mod: $Id$
+ Provides a set of chirp estimation algorithm variants for comparison
+ and characterization. This is canned code meant for testing and
+ exploration of sinusoidal estimation. It's designed to be held still
+ and used for reference and comparison against later improvement.
+ Please, don't optimize this code. Test optimizations elsewhere in
+ code to be compared to this code. Roll new technique improvements
+ into this reference set only after they have been completed.
+
********************************************************************/
#define _GNU_SOURCE
@@ -29,7 +37,6 @@
return (A->A<B->A) - (B->A<A->A);
}
-
static int cmp_ascending_W(const void *p1, const void *p2){
chirp *A = (chirp *)p1;
chirp *B = (chirp *)p2;
@@ -37,216 +44,610 @@
return (B->W<A->W) - (A->W<B->W);
}
-/* performs a second-order nonlinear chirp fit. Passed in chirps are
- used as initial estimation inputs to the iteration. Applies the
- square of passed in window to the input/basis functions.
+/* The stimators project a residue vecotr onto basis functions; the
+ results of the accumulation map directly to the parameters we're
+ interested in (A, P, W, dA, dW, ddA) and vice versa. The mapping
+ equations are broken out below. */
- x: input signal (unwindowed)
- y: output reconstruction (unwindowed)
- window: window to apply to input/basis during fitting process
- len: length of x and window vectors
- c: chirp estimation inputs, chirp fit outputs (sorted in frequency order)
- n: number of chirps
- iter_limit: maximum allowable number of iterations
- fit_limit: desired fit limit
+static float toA(float Ai, float Bi){
+ return hypotf(Ai,Bi);
+}
+static float toP(float Ai, float Bi){
+ return -atan2f(Bi, Ai);
+}
+static float toW(float Ai, float Bi, float Ci, float Di){
+ return (Ai*Ai || Bi*Bi) ? (Ci*Bi - Di*Ai)/(Ai*Ai + Bi*Bi) : 0;
+}
+static float todA(float Ai, float Bi, float Ci, float Di){
+ return (Ai*Ai || Bi*Bi) ? (Ci*Ai + Di*Bi)/hypotf(Ai,Bi) : 0;
+}
+static float todW(float Ai, float Bi, float Ei, float Fi){
+ return (Ai*Ai || Bi*Bi) ? (Ei*Bi - Fi*Ai)/(Ai*Ai + Bi*Bi) : 0;
+}
+static float toddA(float Ai, float Bi, float Ei, float Fi){
+ return (Ai*Ai || Bi*Bi) ? (Ei*Ai + Fi*Bi)/hypotf(Ai,Bi) : 0;
+}
- Chirp frequency must be within ~ 1 FFT bin of fit frequency (see
- paper). Amplitude, amplitude' and phase are fit wouthout a need
- for initial estimate. Frequency' fit will converge from any point,
- though it will fit to the closest lobe (wheter mainlobe or
- sidelobe). All parameters fit within 3 iterations with the
- exception of frequency' which shows approximately linear
- convergence outside of ~ .5 frequency bins.
+static float toAi(float A, float P){
+ return A * cosf(-P);
+}
+static float toBi(float A, float P){
+ return A * sinf(-P);
+}
- Fitting terminates when no fit parameter changes by more than
- fit_limit in an iteration or when the fit loop reaches the
- iteration limit. The fit parameters as checked are scaled over the
- length of the block.
+static float WtoCi(float Ai, float Bi, float W){
+ return W*Bi;
+}
+static float WtoDi(float Ai, float Bi, float W){
+ return -W*Ai;
+}
+static float dWtoEi(float Ai, float Bi, float dW){
+ return dW*Bi;
+}
+static float dWtoFi(float Ai, float Bi, float dW){
+ return -dW*Ai;
+}
-*/
+static float dAtoCi(float Ai, float Bi, float dA){
+ return (Ai*Ai || Bi*Bi) ? dA*Ai/hypotf(Ai,Bi) : 0.;
+}
+static float dAtoDi(float Ai, float Bi, float dA){
+ return (Ai*Ai || Bi*Bi) ? dA*Bi/hypotf(Ai,Bi) : 0.;
+}
+static float ddAtoEi(float Ai, float Bi, float ddA){
+ return (Ai*Ai || Bi*Bi) ? ddA*Ai/hypotf(Ai,Bi) : 0.;
+}
+static float ddAtoFi(float Ai, float Bi, float ddA){
+ return (Ai*Ai || Bi*Bi) ? ddA*Bi/hypotf(Ai,Bi) : 0.;
+}
-int estimate_chirps(const float *x, float *y, float *r, const float *window, int len,
- chirp *c, int n, int iter_limit, float fit_limit){
+/* Nonlinear estimation iterator; functionally, it should differ from
+ the linear version only in that it restarts the basis functions for
+ each chirp after each new fit. */
+static int nonlinear_iterate(const float *x,
+ float *y,
+ const float *window,
+ int len,
+ chirp *cc,
+ int n,
- int i,j,flag=1;
- float E0=0;
- float E1=0;
- float E2=0;
+ float fit_limit,
+ int iter_limit,
- /* Build a starting reconstruction based on the passied in initial
- chirp estimates */
- memset(y,0,sizeof(*y)*len);
- for(i=0;i<n;i++){
- /* initialize phase state for chirp oscillator */
- float dc = cosf(c[i].dW*2.);
- float ds = sinf(c[i].dW*2.);
- float wc = cosf(c[i].W - c[i].dW*len);
- float ws = sinf(c[i].W - c[i].dW*len);
- float cc = cosf((.25*c[i].dW*(len-1) - .5*c[i].W)*(len-1) + c[i].P);
- float cs = sinf((.25*c[i].dW*(len-1) - .5*c[i].W)*(len-1) + c[i].P);
- float tc;
+ int fitW,
+ int fitdA,
+ int fitdW,
+ int fitddA,
+ int symm_norm,
+ float E0,
+ float E1,
+ float E2,
+ int fit_compound,
+ int bound_edges){
+ int i,j;
+ int flag=1;
+ float r[len];
- for(j=0;j<len;j++){
- float jj = j-len*.5+.5;
- float a = c[i].A + c[i].dA*jj;
- y[j] += a*cc;
-
- /* chirp oscillator iteration */
- tc = wc;
- wc = wc*dc - ws*ds;
- ws = tc*ds + ws*dc;
- tc = cc;
- cc = cc*wc - cs*ws;
- cs = tc*ws + cs*wc;
- }
- }
-
- if(iter_limit==0)
- for(j=0;j<len;j++)
- r[j]=(x[j]-y[j])*window[j];
-
- /* this can be moved out/elsewhere/handled differently, as the
- window would normally not change */
- for(i=0;i<len;i++){
- float jj = i-len*.5+.5;
- float w2 = window[i]*window[i];
- float w2j2 = w2*jj*jj;
- float w2j4 = w2j2*jj*jj;
- E0 += w2;
- E1 += w2j2;
- E2 += w2j4;
- }
- E0=2/E0;
- E1=2/E1;
- E2=2/E2;
-
/* outer fit iteration */
- while(flag && iter_limit--){
+ while(flag && iter_limit>0){
+ iter_limit--;
flag=0;
/* precompute the portion of the projection/fit estimate shared by
- the zero, first and second order fits. Subtractsthe current
- best fits from the input signal and widows the result. */
+ the zero, first and second order fits. Subtracts the current
+ best fits from the input signal and windows the result. */
for(j=0;j<len;j++)
r[j]=(x[j]-y[j])*window[j];
memset(y,0,sizeof(*y)*len);
/* Sort chirps by descending amplitude */
- qsort(c, n, sizeof(*c), cmp_descending_A);
+ qsort(cc, n, sizeof(*cc), cmp_descending_A);
- /* reestimate the next chirp fit */
for(i=0;i<n;i++){
- float lA2,lA,lW,ldA,lP,ldW;
+ chirp *c = cc+i;
float aP=0, bP=0;
float cP=0, dP=0;
float eP=0, fP=0;
float cP2=0, dP2=0;
float eP2=0, fP2=0;
float eP3=0, fP3=0;
- float aC = cos(c[i].P);
- float aS = sin(c[i].P);
+ float aE=0, bE=0;
+ float cE=0, dE=0;
+ float eE=0, fE=0;
+ float aC = cos(c->P);
+ float aS = sin(c->P);
- /* initialize phase state for chirp oscillator */
- float dc = cosf(c[i].dW*2.);
- float ds = sinf(c[i].dW*2.);
- float wc = cosf(c[i].W - c[i].dW*len);
- float ws = sinf(c[i].W - c[i].dW*len);
- float cc = cosf((.25*c[i].dW*(len-1) - .5*c[i].W)*(len-1));
- float cs = sinf((.25*c[i].dW*(len-1) - .5*c[i].W)*(len-1));
- float tc;
for(j=0;j<len;j++){
- float jj = j-len*.5+.5;
- float jj2 = jj*jj;
- float co = cc*window[j];
- float si = cs*window[j];
- float c2 = co*co*jj;
- float s2 = si*si*jj;
+ float jj = j-len*.5+.5;
+ float jj2 = jj*jj;
+ float co,si,c2,s2;
+ float yy=r[j];
+ sincosf((c->W + c->dW*jj)*jj,&si,&co);
+ si*=window[j];
+ co*=window[j];
+ c2 = co*co*jj;
+ s2 = si*si*jj;
+
/* add the current estimate back to the residue vector */
- float yy = r[j] += (aC*co-aS*si) * (c[i].A + c[i].dA*jj);
+ r[j] += (aC*co-aS*si) * (c->A + (c->dA + c->ddA*jj)*jj);
- /* chirp oscillator iteration */
- tc = wc;
- wc = wc*dc - ws*ds;
- ws = tc*ds + ws*dc;
- tc = cc;
- cc = cc*wc - cs*ws;
- cs = tc*ws + cs*wc;
-
- /* partial zero order projection */
- aP += co*yy;
- bP += si*yy;
- /* partial first order projection */
+ /* zero order projection */
+ aP += co*yy;
+ bP += si*yy;
+ /* first order projection */
cP += co*yy*jj;
dP += si*yy*jj;
- cP2 += c2;
- dP2 += s2;
- /* partial second order projection */
+ /* second order projection */
eP += co*yy*jj2;
fP += si*yy*jj2;
- eP2 += c2*jj;
- fP2 += s2*jj;
- eP3 += c2*jj2;
- fP3 += s2*jj2;
+
+ if(fit_compound){
+ /* subtract zero order estimate from first */
+ cP2 += c2;
+ dP2 += s2;
+ /* subtract zero order estimate from second */
+ eP2 += c2*jj;
+ fP2 += s2*jj;
+ /* subtract first order estimate from second */
+ eP3 += c2*jj2;
+ fP3 += s2*jj2;
+ }
+
+ if(!symm_norm){
+ /* accumulate per-basis energy for normalization */
+ aE += co*co;
+ bE += si*si;
+ cE += c2*jj;
+ dE += s2*jj;
+ eE += c2*jj*jj2;
+ fE += s2*jj*jj2;
+ }
}
- /* finish projections, scale in terms of total basis energy */
- aP *= E0;
- bP *= E0;
- cP = (cP-aP*cP2)*E1;
- dP = (dP-bP*dP2)*E1;
- eP = (eP-aP*eP2-cP*eP3)*E2;
- fP = (fP-bP*fP2-dP*fP3)*E2;
- /* compute new chirp fit estimate from basis projections */
- lA2 = aP*aP + bP*bP;
- lA = sqrt(lA2);
- lP = -atan2(bP, aP);
- lW = (cP*bP - dP*aP)/lA2;
- ldA = (cP*aP + dP*bP)/lA;
- ldW = 1.75*(eP*bP - fP*aP)/lA2;
+ /* normalize, complete projections */
+ if(symm_norm){
+ aP *= E0;
+ bP *= E0;
+ cP = (cP-aP*cP2)*E1;
+ dP = (dP-bP*dP2)*E1;
+ eP = (eP-aP*eP2-cP*eP3)*E2;
+ fP = (fP-bP*fP2-dP*fP3)*E2;
+ }else{
+ aP = (aE?aP/aE:0.);
+ bP = (bE?bP/bE:0.);
+ cP = (cE?(cP-aP*cP2)/cE:0.);
+ dP = (dE?(dP-bP*dP2)/dE:0.);
+ eP = (eE?(eP-aP*eP2-cP*eP3)/eE:0.);
+ fP = (fE?(fP-bP*fP2-dP*fP3)/fE:0.);
+ }
- /* have we converged to within a requested limit */
- if(c[i].A==0 || fabs((lA-c[i].A)/c[i].A)>fit_limit) flag=1;
- if(fabs(lP - c[i].P)>fit_limit) flag=1;
- if(fabs(lW*len)>fit_limit) flag=1;
- if(c[i].A==0 || fabs((ldA-c[i].dA)/c[i].A)>fit_limit) flag=1;
- if(fabs(ldW*len*len)>fit_limit*2.*M_PI) flag=1;
+ if(!fitW && !fitdA)
+ cP=dP=0;
+ if(!fitdW && !fitddA)
+ eP=fP=0;
- /* save new fit estimate */
- c[i].A = lA;
- c[i].P = lP;
- c[i].W += lW;
- c[i].dA = ldA;
- c[i].dW += ldW;
+ /* we're fitting to the remaining error; add the fit to date
+ back in to relate our newest incremental results to the
+ global fit so far. Note that this does not include W or dW,
+ as they're already 'subtracted' when the bases are recentered
+ each iteration */
+ {
+ float A = toAi(c->A, c->P);
+ float B = toBi(c->A, c->P);
+ float C = dAtoCi(A,B,c->dA) + cP;
+ float D = dAtoDi(A,B,c->dA) + dP;
+ float E = ddAtoEi(A,B,c->ddA) + eP;
+ float F = ddAtoFi(A,B,c->ddA) + fP;
+ A += aP;
+ B += bP;
- /* reinitialize phase state for chirp oscillator */
- dc = cosf(c[i].dW*2.);
- ds = sinf(c[i].dW*2.);
- wc = cosf(c[i].W - c[i].dW*len);
- ws = sinf(c[i].W - c[i].dW*len);
- cc = cosf((.25*c[i].dW*(len-1) - .5*c[i].W)*(len-1) + c[i].P);
- cs = sinf((.25*c[i].dW*(len-1) - .5*c[i].W)*(len-1) + c[i].P);
+ /* base convergence on relative projection movement this
+ iteration */
+ float move = (aP*aP + bP*bP)/(A*A + B*B + fit_limit*fit_limit) +
+ (cP*cP + dP*dP)/(A*A + B*B + fit_limit*fit_limit) +
+ (eP*eP + fP*fP)/(A*A + B*B + fit_limit*fit_limit);
+
+ if(fit_limit>0 && move>fit_limit*fit_limit)flag=1;
+ if(fit_limit<0 && move>1e-14)flag=1;
+
+ aP = A;
+ bP = B;
+ cP = C;
+ dP = D;
+ eP = E;
+ fP = F;
+
+ /* guard overflow; if we're this far out, assume we're never
+ coming back. drop out now. */
+ if((aP*aP + bP*bP)>1e20 ||
+ (cP*cP + dP*dP)>1e20 ||
+ (eP*eP + fP*fP)>1e20){
+ iter_limit=0;
+ i=n;
+ }
+ }
+
+ /* save new estimate */
+ c->A = toA(aP,bP);
+ c->P = toP(aP,bP);
+ c->W += (fitW ? toW(aP,bP,cP,dP) : 0);
+ c->dA = (fitdA ? todA(aP,bP,cP,dP) : 0);
+ c->dW += 1.75*(fitdW ? todW(aP,bP,eP,fP) : 0);
+ c->ddA = (fitddA ? toddA(aP,bP,eP,fP) : 0);
+
+ if(bound_edges){
+ /* do not allow negative or trans-Nyquist center frequencies */
+ if(c->W<0){
+ c->W = 0; /* clamp frequency to 0 (DC) */
+ c->dW = 0; /* if W is 0, the chirp rate must also be 0 to
+ avoid negative frequencies */
+ }
+ if(c->W>M_PI){
+ c->W = M_PI; /* clamp frequency to Nyquist */
+ c->dW = 0; /* if W is at Nyquist, the chirp rate must
+ also be 0 to avoid trans-Nyquist
+ frequencies */
+ }
+ /* Just like with the center frequency, don't allow the
+ chirp rate (dW) parameter to cross over to a negative
+ instantaneous frequency */
+ if(c->W+c->dW*len < 0)
+ c->dW = -c->W/len;
+ if(c->W-c->dW*len < 0)
+ c->dW = c->W/len;
+ /* ...or exceed Nyquist */
+ if(c->W + c->dW*len > M_PI)
+ c->dW = (M_PI - c->W)/len;
+ if(c->W - c->dW*len > M_PI)
+ c->dW = (M_PI + c->W)/len;
+ }
+
/* update the reconstruction/residue vectors with new fit */
for(j=0;j<len;j++){
float jj = j-len*.5+.5;
- float a = c[i].A + c[i].dA*jj;
- float v = a*cc;
-
- /* chirp oscillator iteration */
- tc = wc;
- wc = wc*dc - ws*ds;
- ws = tc*ds + ws*dc;
- tc = cc;
- cc = cc*wc - cs*ws;
- cs = tc*ws + cs*wc;
-
+ float a = c->A + (c->dA + c->ddA*jj)*jj;
+ float v = a*cosf((c->W + c->dW*jj)*jj + c->P);
r[j] -= v*window[j];
y[j] += v;
}
}
}
+ return iter_limit;
+}
+
+/* linear estimation iterator; sets fixed basis functions for each
+ chirp according to the passed in initial estimates. This code looks
+ quite different from the nonlinear estimator above as I wanted to
+ keep the structure very similar to the way it was as originally
+ coded by JM. New techniques are bolted on rather than fully
+ integrated. */
+
+static int linear_iterate(const float *x,
+ float *y,
+ const float *window,
+ int len,
+ chirp *ch,
+ int n,
+
+ float fit_limit,
+ int iter_limit,
+
+ int fitW,
+ int fitdA,
+ int fitdW,
+ int fitddA,
+ int symm_norm,
+ float E0,
+ float E1,
+ float E2,
+ int fit_compound){
+
+ 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 tcosC2[n], tsinC2[n];
+ float ttcosC2[n], ttsinC2[n];
+ float ttcosC3[n], ttsinC3[n];
+
+ float ai[n], bi[n];
+ float ci[n], di[n];
+ float ei[n], fi[n];
+ int i,j;
+ int flag=1;
+
+ for (i=0;i<n;i++){
+ float tmpa=0;
+ float tmpb=0;
+ float tmpc=0;
+ float tmpd=0;
+ float tmpe=0;
+ float tmpf=0;
+
+ float tmpc2=0;
+ float tmpd2=0;
+ float tmpe3=0;
+ float tmpf3=0;
+
+ /* Build a table for the basis functions at each frequency */
+ 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,s;
+ sincosf((ch[i].W + ch[i].dW*jj)*jj,&s,&c);
+
+ /* save basis funcs */
+ 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]*window[j]*window[j];
+ tmpb += sin_table[i][j]*sin_table[i][j]*window[j]*window[j];
+
+ /* The modulation terms */
+ tmpc += tcos_table[i][j]*tcos_table[i][j]*window[j]*window[j];
+ tmpd += tsin_table[i][j]*tsin_table[i][j]*window[j]*window[j];
+
+ /* The second order modulations */
+ tmpe += ttcos_table[i][j]*ttcos_table[i][j]*window[j]*window[j];
+ tmpf += ttsin_table[i][j]*ttsin_table[i][j]*window[j]*window[j];
+
+ /* compound fit terms */
+ tmpc2 += cos_table[i][j]*tcos_table[i][j]*window[j]*window[j];
+ tmpd2 += sin_table[i][j]*tsin_table[i][j]*window[j]*window[j];
+ tmpe3 += tcos_table[i][j]*ttcos_table[i][j]*window[j]*window[j];
+ tmpf3 += tsin_table[i][j]*ttsin_table[i][j]*window[j]*window[j];
+
+ }
+
+ /* Set basis normalizations */
+ if(symm_norm){
+ cosE[i] = sinE[i] = E0;
+ tcosE[i] = tsinE[i] = E1;
+ ttcosE[i] = ttsinE[i] = E2;
+ }else{
+ cosE[i] = (tmpa>0.f ? 1./tmpa : 0);
+ sinE[i] = (tmpb>0.f ? 1./tmpb : 0);
+ tcosE[i] = (tmpc>0.f ? 1./tmpc : 0);
+ tsinE[i] = (tmpd>0.f ? 1./tmpd : 0);
+ ttcosE[i] = (tmpe>0.f ? 1./tmpe : 0);
+ ttsinE[i] = (tmpf>0.f ? 1./tmpf : 0);
+ }
+
+ /* set compound fit terms */
+ tcosC2[i] = tmpc2;
+ tsinC2[i] = tmpd2;
+ ttcosC2[i] = tmpc;
+ ttsinC2[i] = tmpd;
+ ttcosC3[i] = tmpe3;
+ ttsinC3[i] = tmpf3;
+
+ /* seed the basis accumulators from our passed in estimated parameters */
+ /* Don't add in W or dW; these are already included by the basis */
+ {
+ float A = toAi(ch[i].A, ch[i].P);
+ float B = toBi(ch[i].A, ch[i].P);
+ ai[i] = A;
+ bi[i] = B;
+ ci[i] = dAtoCi(A,B,ch[i].dA);
+ di[i] = dAtoDi(A,B,ch[i].dA);
+ ei[i] = ddAtoEi(A,B,ch[i].ddA);
+ fi[i] = ddAtoFi(A,B,ch[i].ddA);
+ }
+ }
+
+ while(flag && iter_limit){
+ iter_limit--;
+ flag=0;
+ for (i=0;i<n;i++){
+
+ float tmpa=0, tmpb=0;
+ float tmpc=0, tmpd=0;
+ float tmpe=0, tmpf=0;
+
+ /* (Sort of) project the residual on the four basis functions */
+ for (j=0;j<len;j++){
+ float r = (x[j]-y[j])*window[j]*window[j];
+ tmpa += r*cos_table[i][j];
+ tmpb += r*sin_table[i][j];
+ tmpc += r*tcos_table[i][j];
+ tmpd += r*tsin_table[i][j];
+ tmpe += r*ttcos_table[i][j];
+ tmpf += r*ttsin_table[i][j];
+ }
+
+ tmpa*=cosE[i];
+ tmpb*=sinE[i];
+
+ if(fit_compound){
+ tmpc -= tmpa*tcosC2[i];
+ tmpd -= tmpb*tsinC2[i];
+ }
+
+ tmpc*=tcosE[i];
+ tmpd*=tsinE[i];
+
+ if(fit_compound){
+ tmpe -= tmpa*ttcosC2[i] + tmpc*ttcosC3[i];
+ tmpf -= tmpb*ttsinC2[i] + tmpd*ttsinC3[i];
+ }
+
+ tmpe*=ttcosE[i];
+ tmpf*=ttsinE[i];
+
+ if(!fitW && !fitdA)
+ tmpc=tmpd=0;
+ if(!fitdW && !fitddA)
+ tmpe=tmpf=0;
+
+ /* Update the signal approximation for the next iteration */
+ for (j=0;j<len;j++){
+ y[j] +=
+ tmpa*cos_table[i][j]+
+ tmpb*sin_table[i][j]+
+ tmpc*tcos_table[i][j]+
+ tmpd*tsin_table[i][j]+
+ tmpe*ttcos_table[i][j]+
+ tmpf*ttsin_table[i][j];
+ }
+
+ ai[i] += tmpa;
+ bi[i] += tmpb;
+ ci[i] += tmpc;
+ di[i] += tmpd;
+ ei[i] += tmpe;
+ fi[i] += tmpf;
+
+ /* base convergence on basis projection movement this
+ iteration */
+ if( fit_limit<0 ||
+ (tmpa*tmpa + tmpb*tmpb)/(ai[i]*ai[i]+bi[i]*bi[i]+fit_limit*fit_limit) +
+ (tmpc*tmpc + tmpd*tmpd)/(ci[i]*ci[i]+di[i]*di[i]+fit_limit*fit_limit) +
+ (tmpe*tmpe + tmpf*tmpf)/(ei[i]*ei[i]+fi[i]*fi[i]+fit_limit*fit_limit) > fit_limit*fit_limit) flag=1;
+
+ }
+ }
+
+ for(i=0;i<n;i++){
+ ch[i].A = toA(ai[i],bi[i]);
+ ch[i].P = toP(ai[i],bi[i]);
+ ch[i].W += (fitW ? toW(ai[i],bi[i],ci[i],di[i]) : 0);
+ ch[i].dA = (fitdA ? todA(ai[i],bi[i],ci[i],di[i]) : 0);
+ ch[i].dW += (fitdW ? todW(ai[i],bi[i],ei[i],fi[i]) : 0);
+ ch[i].ddA = (fitddA ? toddA(ai[i],bi[i],ei[i],fi[i]) : 0);
+
+ free(cos_table[i]);
+ free(sin_table[i]);
+ free(tcos_table[i]);
+ free(tsin_table[i]);
+ free(ttcos_table[i]);
+ free(ttsin_table[i]);
+ }
+
+ return iter_limit;
+}
+
+/* Performs an iterative chirp estimation using the passed
+ in chirp paramters as an initial estimate.
+
+ x: input signal (unwindowed)
+ y: output reconstruction (unwindowed)
+ window: window to apply to input and bases during fitting process
+ len: length of x,y,r and window vectors
+ c: chirp initial estimation inputs, chirp fit outputs (sorted in frequency order)
+ n: number of chirps
+ iter_limit: maximum allowable number of iterations
+ fit_limit: desired fit quality limit
+
+ fitW : flag indicating that fitting should include the W parameter.
+ fitdA : flag indicating that fitting should include the dA parameter.
+ fitdW : flag indicating that fitting should include the dW parameter.
+ fitddA : flag indicating that fitting should include the ddA parameter.
+ symm_norm
+ int fit_compound
+ int bound_zero
+
+ Input estimates affect convergence region and speed and fit
+ quality; precise effects depend on the fit estimation algorithm
+ selected. Note that non-zero initial values for dA, dW or ddA are
+ ignored when fitdA, fitdW, fitddA are unset. An initial value for
+ W is always used even when W is left unfitted.
+
+ Fitting terminates when no fit parameter changes by more than
+ fit_limit in an iteration or when the fit loop reaches the
+ iteration limit. The fit parameters as checked are scaled over the
+ length of the block.
+
+*/
+
+int estimate_chirps(const float *x,
+ float *y,
+ const float *window,
+ int len,
+ chirp *c,
+ int n,
+ int iter_limit,
+ float fit_limit,
+
+ int linear,
+ int fitW,
+ int fitdA,
+ int fitdW,
+ int fitddA,
+ int symm_norm,
+ int fit_compound,
+ int bound_zero
+ ){
+
+ int i,j;
+ float E0=0, E1=0, E2=0;
+
+ if(symm_norm){
+ for(i=0;i<len;i++){
+ float jj = i-len*.5+.5;
+ float w2 = window[i]*window[i];
+ float w2j2 = w2*jj*jj;
+ float w2j4 = w2j2*jj*jj;
+ E0 += w2;
+ E1 += w2j2;
+ E2 += w2j4;
+ }
+ E0=2/E0;
+ E1=2/E1;
+ E2=2/E2;
+ }
+
+ /* zero out initial estimates for params we'll not fit */
+ for(i=0;i<n;i++){
+ if(!fitdA) c[i].dA=0.;
+ if(!fitdW) c[i].dW=0.;
+ if(!fitddA) c[i].ddA=0.;
+ }
+
+ /* Build a starting reconstruction based on the passed-in initial
+ chirp estimates */
+ memset(y,0,sizeof(*y)*len);
+ for(i=0;i<n;i++){
+ for(j=0;j<len;j++){
+ float jj = j-len*.5+.5;
+ float a = c[i].A + (c[i].dA + c[i].ddA*jj)*jj;
+ y[j] += a*cosf((c[i].W + c[i].dW*jj)*jj + c[i].P);
+ }
+ }
+
+ if(linear){
+ if(bound_zero) return -1;
+ iter_limit = linear_iterate(x,y,window,len,c,n,
+ fit_limit,iter_limit,
+ fitW,fitdA,fitdW,fitddA,
+ symm_norm,E0,E1,E2,
+ fit_compound);
+ }else{
+ iter_limit = nonlinear_iterate(x,y,window,len,c,n,
+ fit_limit,iter_limit,
+ fitW,fitdA,fitdW,fitddA,
+ symm_norm,E0,E1,E2,
+ fit_compound,bound_zero);
+ }
+
/* Sort by ascending frequency */
qsort(c, n, sizeof(*c), cmp_ascending_W);
@@ -265,50 +666,19 @@
}
}
-/* OMG! An example! */
-
-#if 0
-
-void hanning(float *x, int n){
- float scale = 2*M_PI/n;
- int i;
- for(i=0;i<n;i++){
- float i5 = i+.5;
- x[i] = .5-.5*cos(scale*i5);
+/* dead, but leave it around; does a slow, brute force DCFT */
+void slow_dcft_row(float *x, float *y, int k, int n){
+ int i,j;
+ for(i=0;i<n/2+1;i++){ /* frequency loop */
+ float real=0.,imag=0.;
+ for(j=0;j<n;j++){ /* multiply loop */
+ float s,c;
+ float jj = (j-n/2+.5)/n;
+ sincosf(2.*M_PI*(i+k*jj)*jj,&s,&c);
+ real += c*x[j];
+ imag += s*x[j];
+ }
+ y[i] = hypot(real,imag);
}
}
-int main(){
- int BLOCKSIZE=1024;
- int i;
- float w[1024];
- float x[1024];
- float y[1024];
- float r[1024];
-
- chirp c[]={
- {.3, 100.*2.*M_PI/BLOCKSIZE, .2, 0, 0},
- {1., 95.*2.*M_PI/BLOCKSIZE, .2, 0, 2.5*2.*M_PI/BLOCKSIZE/BLOCKSIZE}
- };
-
- hanning(w,BLOCKSIZE);
-
- for(i=0;i<BLOCKSIZE;i++){
- float jj = i-BLOCKSIZE*.5+.5;
- x[i]=(.3+.1*jj/BLOCKSIZE)*cos(jj*2.*M_PI/BLOCKSIZE * 100.2 +1.2);
- x[i]+=cos(jj*2.*M_PI/BLOCKSIZE * (96. + 3./BLOCKSIZE*jj) +1.2);
- }
-
- int ret=estimate_chirps(x,y,r,w,BLOCKSIZE,c,2,55,.01);
-
- for(i=0;i<sizeof(c)/sizeof(*c);i++)
- fprintf(stderr,"W:%f\nP:%f\nA:%f\ndW:%f\ndA:%f\nreturned=%d\n\n",
- c[i].W*BLOCKSIZE/2./M_PI,
- c[i].P,
- c[i].A,
- c[i].dW*BLOCKSIZE*BLOCKSIZE/2./M_PI,
- c[i].dA*BLOCKSIZE,ret);
-
- return 0;
-}
-#endif
Modified: trunk/ghost/monty/chirp/chirp.h
===================================================================
--- trunk/ghost/monty/chirp/chirp.h 2011-03-28 18:22:35 UTC (rev 17919)
+++ trunk/ghost/monty/chirp/chirp.h 2011-04-12 06:45:19 UTC (rev 17920)
@@ -21,12 +21,27 @@
float P; /* phase (radians) */
float dA; /* amplitude modulation (linear change per sample) */
float dW; /* frequency modulation (radians per sample^2) */
+ float ddA; /* amplitude modulation (linear change per sample^2) */
int label;/* used for tracking by outside code */
} chirp;
-extern int estimate_chirps(const float *x, float *y, float *r,
- const float *window, int len,
- chirp *c, int n, int iter_limit,
- float fit_limit);
+extern int estimate_chirps(const float *x,
+ float *y,
+ const float *window,
+ int len,
+ chirp *c,
+ int n,
+ int iter_limit,
+ float fit_limit,
+
+ int linear,
+ int fitW,
+ int fitdA,
+ int fitdW,
+ int fitddA,
+ int symm_norm,
+ int fit_compound,
+ int bound_zero);
+
extern void advance_chirps(chirp *c, int n, int len);
Added: trunk/ghost/monty/chirp/chirpgraph.c
===================================================================
--- trunk/ghost/monty/chirp/chirpgraph.c (rev 0)
+++ trunk/ghost/monty/chirp/chirpgraph.c 2011-04-12 06:45:19 UTC (rev 17920)
@@ -0,0 +1,830 @@
+/********************************************************************
+ * *
+ * THIS FILE IS PART OF THE OggGhost SOFTWARE CODEC SOURCE CODE. *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS *
+ * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
+ * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. *
+ * *
+ * THE OggGhost SOURCE CODE IS (C) COPYRIGHT 2007-2011 *
+ * by the Xiph.Org Foundation http://www.xiph.org/ *
+ * *
+ ********************************************************************
+
+ function: graphing code for chirp tests
+ last mod: $Id$
+
+ ********************************************************************/
+
+#define _GNU_SOURCE
+#include <math.h>
+#include "chirp.h"
+#include "scales.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <cairo/cairo.h>
+#include <pthread.h>
+
+/* configuration globals. */
+int BLOCKSIZE=2048;
+int cores=8;
+int xstepd=100;
+int ystepd=100;
+int xstepg=25;
+int ystepg=25;
+int x_n=801;
+int y_n=300;
+int fontsize=12;
+
+int randomize_est_A=0;
+int randomize_est_P=0;
+int randomize_est_dA=0;
+int randomize_est_dW=0;
+int randomize_est_ddA=0;
+
+int randomize_chirp_P=0;
+int randomize_chirp_dA=0;
+int randomize_chirp_dW=0;
+int randomize_chirp_ddA=0;
+
+float initial_est_A = 0.;
+float initial_est_P = 0.;
+float initial_est_dA = 0.;
+float initial_est_dW = 0.;
+float initial_est_ddA = 0.;
+
+float initial_chirp_A = 1.;
+float initial_chirp_P = M_PI/4;
+float initial_chirp_dA = 0.;
+float initial_chirp_dW = 0.;
+float initial_chirp_ddA = 0.;
+
+int fit_linear=0;
+int fit_W=1;
+int fit_dA=1;
+int fit_dW=1;
+int fit_ddA=1;
+int fit_symm_norm=1;
+int fit_gauss_seidel=1;
+int fit_bound_zero=0;
+float fit_tolerance=.000001;
+int fit_hanning=1;
+int iterations=-1;
+
+void hanning(float *x, int n){
+ float scale = 2*M_PI/n;
+ int i;
+ for(i=0;i<n;i++){
+ float i5 = i+.5;
+ x[i] = .5-.5*cos(scale*i5);
+ }
+}
+
+/********************** colors for graphs *********************************/
+
+void set_iter_color(cairo_t *cC, int ret, float a){
+ if (ret>100){
+ cairo_set_source_rgba(cC,1,1,1,a); /* white */
+ }else if (ret>50){ /* 51(black)-100(white) */
+ cairo_set_source_rgba(cC,(ret-50)*.02,(ret-50)*.02,(ret-50)*.02,a);
+ }else if (ret>30){ /* 31(red) - 50(black) */
+ cairo_set_source_rgba(cC,1-(ret-30)*.05,0,0,a);
+ }else if (ret>10){ /* 11(yellow) - 30(red) */
+ cairo_set_source_rgba(cC,1,1-(ret-10)*.05,0,a);
+ }else if (ret>4){ /* 5 (green) - 10(yellow) */
+ cairo_set_source_rgba(cC,(ret-5)*.15+.25,1,0,a);
+ }else if (ret==4){ /* 4 brighter cyan */
+ cairo_set_source_rgba(cC,0,.8,.4,a);
+ }else if (ret==3){ /* 3 cyan */
+ cairo_set_source_rgba(cC,0,.6,.6,a);
+ }else if (ret==2){ /* 2 dark cyan */
+ cairo_set_source_rgba(cC,0,.4,.8,a);
+ }else{ /* 1 blue */
+ cairo_set_source_rgba(cC,.1,.1,1,a);
+ }
+}
+
+void set_error_color(cairo_t *c, float err,float a){
+ if(isnan(err) || fabs(err)>1.){
+ cairo_set_source_rgba(c,1,1,1,a); /* white */
+ }else{
+ err=fabs(err);
+ if(err>.1){
+ cairo_set_source_rgba(c,(err-.1)/.9,(err-.1)/.9,(err-.1)/.9,a); /* white->black */
+ }else if(err>.01){
+ cairo_set_source_rgba(c,1.-((err-.01)/.09),0,0,a); /* black->red */
+ }else if(err>.001){
+ cairo_set_source_rgba(c,1,1.-((err-.001)/.009),0,a); /* red->yellow */
+ }else if(err>.0001){
+ cairo_set_source_rgba(c,((err-.0001)/.0009),1,0,a); /* yellow->green */
+ }else if(err>.00001){
+ cairo_set_source_rgba(c,0,1,1.-((err-.00001)/.00009),a); /* green->cyan */
+ }else if(err>.000001){
+ cairo_set_source_rgba(c,.2-((err-.000001)/.000009)*.2,
+ ((err-.000001)/.000009)*.8+.2,1,a); /* cyan->blue */
+ }else{
+ cairo_set_source_rgba(c,.1,.1,1,a); /* blue */
+ }
+ }
+}
+
+/********* draw everything in the graph except the graph data itself *******/
+
+#define DT_iterations 0
+#define DT_abserror 1
+#define DT_percent 2
+
+/* computed configuration globals */
+int leftpad=0;
+int rightpad=0;
+int toppad=0;
+int bottompad=0;
+float legendy=0;
+float legendh=0;
+int pic_w=0;
+int pic_h=0;
+float titleheight=0.;
+
+/* determines padding, etc. Call several times to determine maximums. */
+void setup_graph(char *title1,
+ char *title2,
+ char *title3,
+ char *xaxis_label,
+ char *yaxis_label,
+ char *legend_label,
+ int datatype){
+
+ /* determine ~ padding needed */
+ cairo_surface_t *cs=cairo_image_surface_create(CAIRO_FORMAT_RGB24,100,100);
+ cairo_t *ct=cairo_create(cs);
+ int minx=-leftpad,maxx=x_n+rightpad;
+ int miny=-toppad,maxy=2*y_n+1+bottompad;
+ int textpad=fontsize;
+ cairo_text_extents_t extents;
+ int x,y;
+
+ /* y labels */
+ cairo_set_font_size(ct, fontsize);
+ for(y=0;y<=y_n;y+=ystepd){
+ char buf[80];
+ snprintf(buf,80,"%.0f",(float)y/ystepd);
+ cairo_text_extents(ct, buf, &extents);
+ if(-textpad-extents.height*.5<miny)miny=-textpad-extents.height*.5;
+ if(-textpad-extents.width-extents.height*1.5<minx)minx=-textpad-extents.width-extents.height*1.5;
+
+ if(y>0 && y<y_n){
+ snprintf(buf,80,"%.0f",(float)-y/ystepd);
+ cairo_text_extents(ct, buf, &extents);
+ if(-textpad-extents.height*.5<miny)miny=-textpad-extents.height*.5;
+ if(-textpad-extents.width<minx)minx=-textpad-extents.width;
+ }
+ }
+
+ /* x labels */
+ for(x=0;x<=x_n;x+=xstepd){
+ char buf[80];
+ snprintf(buf,80,"%.0f",(float)x/xstepd);
+ cairo_text_extents(ct, buf, &extents);
+
+ if(y_n*2+1+extents.height*3>maxy)maxy=y_n*2+1+extents.height*3;
+ if(x-(extents.width/2)-textpad<minx)minx=x-(extents.width/2)-textpad;
+ if(x+(extents.width/2)+textpad>maxx)maxx=x+(extents.width/2)+textpad;
+ }
+
+ /* center horizontally */
+ if(maxx-x_n < -minx)maxx = x_n-minx;
+ if(maxx-x_n > -minx)minx = x_n-maxx;
+
+ /* top legend */
+ {
+ float sofar=0;
+ if(title1){
+ cairo_save(ct);
+ cairo_select_font_face (ct, "",
+ CAIRO_FONT_SLANT_NORMAL,
+ CAIRO_FONT_WEIGHT_BOLD);
+ cairo_set_font_size(ct,fontsize*1.25);
+ cairo_text_extents(ct, title1, &extents);
+
+ if(extents.width+textpad > maxx-minx){
+ int moar = extents.width+textpad-maxx+minx;
+ minx-=moar/2;
+ maxx+=moar/2;
+ }
+ sofar+=extents.height;
+ cairo_restore(ct);
+ }
+ if(title2){
+ cairo_text_extents(ct, title2, &extents);
+
+ if(extents.width+textpad > maxx-minx){
+ int moar = extents.width+textpad-maxx+minx;
+ minx-=moar/2;
+ maxx+=moar/2;
+ }
+ sofar+=extents.height;
+ }
+ if(title3){
+ cairo_text_extents(ct, title3, &extents);
+
+ if(extents.width+textpad > maxx-minx){
+ int moar = extents.width+textpad-maxx+minx;
+ minx-=moar/2;
+ maxx+=moar/2;
+ }
+ sofar+=extents.height;
+ }
+ if(sofar>titleheight)titleheight=sofar;
+ if(-miny<titleheight+textpad*2)miny=-titleheight-textpad*2;
+
+ }
+
+ /* color legend */
+ {
+ float width=0.;
+ float w1,w2,w3;
+ cairo_text_extents(ct, legend_label, &extents);
+ width=extents.width;
+
+ cairo_text_extents(ct, "100", &extents);
+ w1=extents.width*2*11;
+ cairo_text_extents(ct, ".000001", &extents);
+ w2=extents.width*1.5*7;
+ cairo_text_extents(ct, ".0001%", &extents);
+ w3=extents.width*1.5*7;
+
+ width+=MAX(w1,MAX(w2,w3));
+
+ legendy = y_n*2+1+extents.height*4;
+ legendh = extents.height*1.5;
+
+ if(width + textpad > x_n-minx){
+ int moar = width+textpad-x_n;
+ minx-=moar;
+ maxx+=moar;
+ }
+ if(legendy+legendh*4>maxy)
+ maxy=legendy+legendh*2.5;
+ }
+
+ toppad = -miny;
+ leftpad = -minx;
+ rightpad = maxx-x_n;
+ bottompad = maxy-(y_n*2+1);
+ if(toppad<leftpad*.75)toppad = leftpad*.75;
+ if(leftpad<toppad)leftpad=toppad;
+ if(rightpad<toppad)rightpad=toppad;
+
+ pic_w = maxx-minx;
+ pic_h = maxy-miny;
+}
+
+/* draws the page surrounding the graph data itself */
+void draw_page(cairo_t *c,
+ char *title1,
+ char *title2,
+ char *title3,
+ char *xaxis_label,
+ char *yaxis_label,
+ char *legend_label,
+ int datatype){
+
+ int i;
+ cairo_text_extents_t extents;
+ cairo_save(c);
+
+ /* clear page to white */
+ cairo_set_source_rgb(c,1,1,1);
+ cairo_rectangle(c,0,0,pic_w,pic_h);
+ cairo_fill(c);
+
+ /* set graph area to transparent */
+ cairo_set_source_rgba(c,0,0,0,0);
+ cairo_set_operator(c,CAIRO_OPERATOR_SOURCE);
+ cairo_rectangle(c,leftpad,toppad,x_n,y_n*2+1);
+ cairo_fill(c);
+ cairo_restore(c);
+
+ /* Y axis numeric labels */
+ cairo_set_source_rgb(c,0,0,0);
+ for(i=0;i<=y_n;i+=ystepd){
+ int y = toppad+y_n;
+ char buf[80];
+
+ snprintf(buf,80,"%.0f",(float)i/ystepd);
+ cairo_text_extents(c, buf, &extents);
+ cairo_move_to(c,leftpad - fontsize*.5 - extents.width,y-i+extents.height*.5);
+ cairo_show_text(c,buf);
+
+ if(i>0 && i<y_n){
+
+ snprintf(buf,80,"%.0f",(float)-i/ystepd);
+ cairo_text_extents(c, buf, &extents);
+ cairo_move_to(c,leftpad - fontsize*.5 - extents.width,y+i+extents.height*.5);
+ cairo_show_text(c,buf);
+
+ }
+ }
+
+ /* X axis labels */
+ for(i=0;i<=x_n;i+=xstepd){
+ char buf[80];
+
+ if(i==0){
+ snprintf(buf,80,"DC");
+ }else{
+ snprintf(buf,80,"%.0f",(float)i/xstepd);
+ }
+ cairo_text_extents(c, buf, &extents);
+ cairo_move_to(c,leftpad + i - extents.width/2,y_n*2+1+toppad+extents.height+fontsize*.5);
+ cairo_show_text(c,buf);
+ }
+
+ /* Y axis label */
+ {
+ cairo_matrix_t a;
+ cairo_matrix_t b = {0.,-1., 1.,0., 0.,0.}; // account for border!
+ cairo_matrix_t d;
+ cairo_text_extents(c, yaxis_label, &extents);
+ cairo_move_to(c,leftpad-extents.height*2,y_n+toppad+extents.width*.5);
+
+ cairo_save(c);
+ cairo_get_matrix(c,&a);
+ cairo_matrix_multiply(&d,&a,&b);
+ cairo_set_matrix(c,&d);
+ cairo_show_text(c,yaxis_label);
+ cairo_restore(c);
+ }
+
+ /* X axis caption */
+ {
+ cairo_text_extents(c, xaxis_label, &extents);
+ cairo_move_to(c,pic_w/2-extents.width/2,y_n*2+1+toppad+extents.height*2+fontsize*.5);
+ cairo_show_text(c,xaxis_label);
+ }
+
+ /* top title(s) */
+ {
+ float y = (toppad-titleheight);
+ if(title1){
+ cairo_save(c);
+ cairo_select_font_face (c, "",
+ CAIRO_FONT_SLANT_NORMAL,
+ CAIRO_FONT_WEIGHT_BOLD);
+ cairo_set_font_size(c,fontsize*1.25);
+ cairo_text_extents(c, title1, &extents);
+ cairo_move_to(c,pic_w/2-extents.width/2,y);
+ cairo_show_text(c,title1);
+ y+=extents.height;
+ cairo_restore(c);
+ }
+ if(title2){
+ cairo_text_extents(c, title2, &extents);
+ cairo_move_to(c,pic_w/2-extents.width/2,y);
+ cairo_show_text(c,title2);
+ y+=extents.height;
+ }
+ if(title3){
+ cairo_text_extents(c, title3, &extents);
+ cairo_move_to(c,pic_w/2-extents.width/2,y);
+ cairo_show_text(c,title3);
+ y+=extents.height;
+ }
+ }
+
+ /* color legend */
+ {
+ float cw;
+ cairo_text_extents(c, "100", &extents);
+ cw=extents.width*2*11;
+ cairo_text_extents(c, ".000001", &extents);
+ cw=MAX(cw,extents.width*1.5*7);
+ cairo_text_extents(c, ".0001%", &extents);
+ cw=MAX(cw,extents.width*1.5*7);
+
+ if(datatype==DT_iterations){
+ char buf[80];
+ int i;
+ float w=cw/11,px=leftpad+x_n;
+
+ for(i=1;i<=100;i++){
+ switch(i){
+ case 1:
+ case 2:
+ case 3:
+ case 4:
+ case 5:
+ case 10:
+ case 20:
+ case 30:
+ case 40:
+ case 50:
+ case 100:
+ px-=w;
+ cairo_rectangle(c,px,legendy+toppad,w,legendh*1.25);
+ cairo_set_source_rgb(c,0,0,0);
+ cairo_stroke_preserve(c);
+ set_iter_color(c, i, 1.);
+ cairo_fill(c);
+
+ snprintf(buf,80,"%d",i);
+ cairo_text_extents(c, buf, &extents);
+ cairo_move_to(c,px+w/2-extents.width*.5,legendy+toppad+legendh*.625+extents.height/2);
+ cairo_set_source_rgba(c,1,1,1,.7);
+ cairo_text_path (c, buf);
+ cairo_stroke_preserve(c);
+ cairo_set_source_rgb(c,0,0,0);
+ cairo_fill(c);
+
+ break;
+ }
+ }
+ cairo_text_extents(c, legend_label, &extents);
+ cairo_move_to(c,px-extents.width-legendh*.75,toppad+legendy+legendh*.625+extents.height*.5);
+ cairo_show_text(c,legend_label);
+ }else{
+ int per_p = (datatype==DT_percent);
+ int i;
+ float w=cw/7,px=leftpad+x_n;
+
+ for(i=0;i<7;i++){
+ char *buf;
+
+ px-=w;
+ cairo_rectangle(c,px,legendy+toppad,w,legendh*1.25);
+ cairo_set_source_rgb(c,0,0,0);
+ cairo_stroke_preserve(c);
+
+ switch(i){
+ case 0:
+ buf=(per_p?".0001%":".000001");
+ set_error_color(c, 0., 1.);
+ break;
+ case 1:
+ buf=(per_p?".001%":".00001");
+ set_error_color(c, .00001, 1.);
+ break;
+ case 2:
+ buf=(per_p?".01%":".0001");
+ set_error_color(c, .0001, 1.);
+ break;
+ case 3:
+ buf=(per_p?".1%":".001");
+ set_error_color(c, .001, 1.);
+ break;
+ case 4:
+ buf=(per_p?"1%":".01");
+ set_error_color(c, .01, 1.);
+ break;
+ case 5:
+ buf=(per_p?"10%":".1");
+ set_error_color(c, .1, 1.);
+ break;
+ case 6:
+ buf=(per_p?"100%":"1");
+ set_error_color(c, 1., 1.);
+ break;
+ }
+ cairo_fill(c);
+
+ cairo_text_extents(c, buf, &extents);
+ cairo_move_to(c,px+w/2-extents.width*.5,legendy+toppad+legendh*.625+extents.height/2);
+ cairo_set_source_rgba(c,1,1,1,.7);
+ cairo_text_path (c, buf);
+ cairo_stroke_preserve(c);
+ cairo_set_source_rgb(c,0,0,0);
+ cairo_fill(c);
+ }
+ cairo_text_extents(c, legend_label, &extents);
+ cairo_move_to(c,px-extents.width-legendh*.75,toppad+legendy+legendh*.625+extents.height*.5);
+ cairo_show_text(c,legend_label);
+ }
+ }
+}
+
+void to_png(cairo_surface_t *c,char *base, char *name){
+ if(c){
+ char buf[320];
+ snprintf(buf,320,"%s-%s.png",name,base);
+ cairo_surface_write_to_png(c,buf);
+ }
+}
+
+/*********************** Plot w estimate error against w **********************/
+
+typedef struct {
+ float *in;
+ float *window;
+
+ chirp *c;
+ int *ret;
+ int max_y;
+} vearg;
+
+
+pthread_mutex_t ymutex = PTHREAD_MUTEX_INITIALIZER;
+int next_y=0;
+
+#define _GNU_SOURCE
+#include <fenv.h>
+
+void *w_e_column(void *in){
+ float rec[BLOCKSIZE];
+ vearg *arg = (vearg *)in;
+ int y;
+ chirp save;
+
+ while(1){
+ pthread_mutex_lock(&ymutex);
+ y=next_y;
+ if(y>=arg->max_y){
+ pthread_mutex_unlock(&ymutex);
+ return NULL;
+ }
+ next_y++;
+ pthread_mutex_unlock(&ymutex);
+
+ int except=feenableexcept(FE_ALL_EXCEPT);
+ fedisableexcept(FE_INEXACT);
+ fedisableexcept(FE_UNDERFLOW);
+ save=arg->c[y];
+
+ /* iterate only until convergence */
+ arg->ret[y]=
+ estimate_chirps(arg->in,rec,arg->window,BLOCKSIZE,
+ arg->c+y,1,iterations,fit_tolerance,
+ fit_linear,
+ fit_W,
+ fit_dA,
+ fit_dW,
+ fit_ddA,
+ fit_symm_norm,
+ fit_gauss_seidel,
+ fit_bound_zero);
+
+ /* continue iterating to get error numbers for a fixed large
+ number of iterations. The linear estimator must be restarted
+ from the beginning, else 'continuing' causes it to recenter the
+ basis-- which renders it nonlinear
+ if(fit_linear) arg->c[y]=save;
+ int ret=estimate_chirps(arg->in,rec,arg->window,BLOCKSIZE,
+ arg->c+y,1,fit_linear?iterations:arg->ret[y],-1.,
+ fit_linear,
+ fit_W,
+ fit_dA,
+ fit_dW,
+ fit_ddA,
+ fit_symm_norm,
+ fit_gauss_seidel,
+ fit_bound_zero);*/
+ arg->ret[y] = iterations-arg->ret[y];
+ feclearexcept(FE_ALL_EXCEPT);
+ feenableexcept(except);
+
+ }
+ return NULL;
+}
+
+void w_e(){
+ float window[BLOCKSIZE];
+ float in[BLOCKSIZE];
+ int i,x,y;
+
+ cairo_surface_t *converge=NULL;
+ cairo_surface_t *Aerror=NULL;
+ cairo_surface_t *Perror=NULL;
+ cairo_surface_t *Werror=NULL;
+ cairo_surface_t *dAerror=NULL;
+ cairo_surface_t *dWerror=NULL;
+ cairo_surface_t *ddAerror=NULL;
+
+ cairo_t *cC=NULL;
+ cairo_t *cA=NULL;
+ cairo_t *cP=NULL;
+ cairo_t *cW=NULL;
+ cairo_t *cdA=NULL;
+ cairo_t *cdW=NULL;
+ cairo_t *cddA=NULL;
+
+ pthread_t threads[cores];
+ vearg arg[cores];
+
+ char *filebase="W-vs-Westimate";
+ char *yaxis_label = "initial distance from W (cycles/block)";
+ char *xaxis_label = "W (cycles/block)";
+ char *title2=NULL;
+ char *title3=NULL;
+
+ pic_w = x_n;
+ pic_h = y_n*2+1;
+ if(iterations<0)iterations=100;
+
+ /* determine ~ padding needed */
+
+ setup_graph("Convergence",title2,title3,xaxis_label,yaxis_label,
+ "Iterations:",DT_iterations);
+ setup_graph("A (Amplitude) Error",title2,title3,xaxis_label,yaxis_label,
+ "Percentage Error:",DT_percent);
+ setup_graph("P (Phase) Error",title2,title3,xaxis_label,yaxis_label,
+ "Error (rads/block)",DT_abserror);
+ if(fit_W)
+ setup_graph("W (Frequency) Error",title2,title3,xaxis_label,yaxis_label,
+ "Error (cycles/block)",DT_abserror);
+ if(fit_dA)
+ setup_graph("dA (Amplitude Modulation) Error",title2,title3,xaxis_label,yaxis_label,
+ "Percentage Error:",DT_percent);
+ if(fit_dW)
+ setup_graph("dW (Chirp Rate) Error",title2,title3,xaxis_label,yaxis_label,
+ "Error (cycles/block)",DT_abserror);
+ if(fit_ddA)
+ setup_graph("ddA (Amplitude Modulation Squared) Error",title2,title3,xaxis_label,yaxis_label,
+ "Percentage Error:",DT_percent);
+
+ /* Make cairo drawables */
+ converge=cairo_image_surface_create(CAIRO_FORMAT_ARGB32,pic_w, pic_h);
+ if(!converge || cairo_surface_status(converge)!=CAIRO_STATUS_SUCCESS){
+ fprintf(stderr,"Could not set up Cairo surface.\n\n");
+ exit(1);
+ }
+ cC=cairo_create(converge);
+ draw_page(cC,"Convergence",title2,title3,xaxis_label,yaxis_label,
+ "Iterations:",DT_iterations);
+
+ Aerror=cairo_image_surface_create(CAIRO_FORMAT_ARGB32,pic_w, pic_h);
+ if(!Aerror || cairo_surface_status(Aerror)!=CAIRO_STATUS_SUCCESS){
+ fprintf(stderr,"Could not set up Cairo surface.\n\n");
+ exit(1);
+ }
+ cA=cairo_create(Aerror);
+ draw_page(cA,"A (Amplitude) Error",title2,title3,xaxis_label,yaxis_label,
+ "Percentage Error:",DT_percent);
+
+ Perror=cairo_image_surface_create(CAIRO_FORMAT_ARGB32,pic_w, pic_h);
+ if(!Perror || cairo_surface_status(Perror)!=CAIRO_STATUS_SUCCESS){
+ fprintf(stderr,"Could not set up Cairo surface.\n\n");
+ exit(1);
+ }
+ cP=cairo_create(Perror);
+ draw_page(cP,"P (Phase) Error",title2,title3,xaxis_label,yaxis_label,
+ "Error (rads/block)",DT_abserror);
+
+ if(fit_W){
+ Werror=cairo_image_surface_create(CAIRO_FORMAT_ARGB32,pic_w, pic_h);
+ if(!Werror || cairo_surface_status(Werror)!=CAIRO_STATUS_SUCCESS){
+ fprintf(stderr,"Could not set up Cairo surface.\n\n");
+ exit(1);
+ }
+ cW=cairo_create(Werror);
+ draw_page(cW,"W (Frequency) Error",title2,title3,xaxis_label,yaxis_label,
+ "Error (cycles/block)",DT_abserror);
+ }
+
+ if(fit_dA){
+ dAerror=cairo_image_surface_create(CAIRO_FORMAT_ARGB32,pic_w, pic_h);
+ if(!dAerror || cairo_surface_status(dAerror)!=CAIRO_STATUS_SUCCESS){
+ fprintf(stderr,"Could not set up Cairo surface.\n\n");
+ exit(1);
+ }
+ cdA=cairo_create(dAerror);
+ draw_page(cdA,"dA (Amplitude Modulation) Error",
+ title2,title3,xaxis_label,yaxis_label,
+ "Percentage Error:",DT_percent);
+ }
+
+ if(fit_dW){
+ dWerror=cairo_image_surface_create(CAIRO_FORMAT_ARGB32,pic_w, pic_h);
+ if(!dWerror || cairo_surface_status(dWerror)!=CAIRO_STATUS_SUCCESS){
+ fprintf(stderr,"Could not set up Cairo surface.\n\n");
+ exit(1);
+ }
+ cdW=cairo_create(dWerror);
+ draw_page(cdW,"dW (Chirp Rate) Error",title2,title3,xaxis_label,yaxis_label,
+ "Error (cycles/block)",DT_abserror);
+ }
+
+ if(fit_ddA){
+ ddAerror=cairo_image_surface_create(CAIRO_FORMAT_ARGB32,pic_w, pic_h);
+ if(!ddAerror || cairo_surface_status(ddAerror)!=CAIRO_STATUS_SUCCESS){
+ fprintf(stderr,"Could not set up Cairo surface.\n\n");
+ exit(1);
+ }
+ cddA=cairo_create(ddAerror);
+ draw_page(cddA,"ddA (Amplitude Modulation Squared) Error",
+ title2,title3,xaxis_label,yaxis_label,
+ "Percentage Error:",DT_percent);
+ }
+
+ for(i=0;i<BLOCKSIZE;i++)
+ window[i]=1.;
+ if(fit_hanning)
+ hanning(window,BLOCKSIZE);
+
+ /* graph computation */
+ for(x=0;x<x_n;x++){
+ float w = (float)x/xstepd;
+ chirp ch[y_n*2+1];
+ int ret[y_n*2+1];
+
+ float rand_P_ch = 2*M_PI*(float)drand48()-M_PI;
+
+ for(i=0;i<BLOCKSIZE;i++){
+ float jj = i-BLOCKSIZE*.5+.5;
+ in[i]=initial_chirp_A * cos(w*jj*2.*M_PI/BLOCKSIZE+initial_chirp_P);
+ }
+
+ fprintf(stderr,"\rW estimate distance vs. W graphs: %d%%...",100*x/x_n);
+
+ for(y=y_n;y>=-y_n;y--){
+ float rand_A_est = (float)drand48();
+ float rand_P_est = 2*M_PI*(float)drand48()-M_PI;
+ int yi=y_n-y;
+ float we=(float)y/ystepd+w;
+ ch[yi]=(chirp){0,
+ (we)*2.*M_PI/BLOCKSIZE,
+ 0,
+ initial_est_dA,
+ initial_est_dW,
+ initial_est_ddA,
+ yi};
+ }
+ next_y=0;
+ for(y=0;y<cores;y++){
+ arg[y].in = in;
+ arg[y].window=window;
+ arg[y].c=ch;
+ arg[y].ret=ret;
+ arg[y].max_y=y_n*2+1;
+ pthread_create(threads+y,NULL,w_e_column,arg+y);
+ }
+ for(y=0;y<cores;y++)
+ pthread_join(threads[y],NULL);
+
+ for(y=-y_n;y<=y_n;y++){
+ int yi=y+y_n;
+ float a = (x%xstepg==0 || y%ystepg==0 ?
+ (x%xstepd==0 || y%ystepd==0 ? .3 : .8) : 1.);
+
+ set_iter_color(cC,ret[yi],a);
+ cairo_rectangle(cC,x+leftpad,yi+toppad,1,1);
+ cairo_fill(cC);
+
+ set_error_color(cA,fabs(initial_chirp_A-ch[yi].A),a);
+ cairo_rectangle(cA,x+leftpad,yi+toppad,1,1);
+ cairo_fill(cA);
+
+ set_error_color(cP,fabs(initial_chirp_P-ch[yi].P),a);
+ cairo_rectangle(cP,x+leftpad,yi+toppad,1,1);
+ cairo_fill(cP);
+
+ if(cW){
+ set_error_color(cW,fabs(w-(ch[yi].W/2./M_PI*BLOCKSIZE)),a);
+ cairo_rectangle(cW,x+leftpad,yi+toppad,1,1);
+ cairo_fill(cW);
+ }
+
+ if(cdA){
+ set_error_color(cdA,ch[yi].dA*BLOCKSIZE,a);
+ cairo_rectangle(cdA,x+leftpad,yi+toppad,1,1);
+ cairo_fill(cdA);
+ }
+
+ if(cdW){
+ set_error_color(cdW,ch[yi].dW/M_PI*BLOCKSIZE*BLOCKSIZE,a);
+ cairo_rectangle(cdW,x+leftpad,yi+toppad,1,1);
+ cairo_fill(cdW);
+ }
+
+ if(cddA){
+ set_error_color(cddA,ch[yi].ddA/BLOCKSIZE*BLOCKSIZE,a);
+ cairo_rectangle(cddA,x+leftpad,yi+toppad,1,1);
+ cairo_fill(cddA);
+ }
+ }
+
+ if((x&15)==0){
+ to_png(converge,filebase,"converge");
+ to_png(Aerror,filebase,"Aerror");
+ to_png(Perror,filebase,"Perror");
+ to_png(Werror,filebase,"Werror");
+ to_png(dAerror,filebase,"dAerror");
+ to_png(dWerror,filebase,"dWerror");
+ to_png(ddAerror,filebase,"ddAerror");
+ }
+
+ }
+
+ to_png(converge,filebase,"converge");
+ to_png(Aerror,filebase,"Aerror");
+ to_png(Perror,filebase,"Perror");
+ to_png(Werror,filebase,"Werror");
+ to_png(dAerror,filebase,"dAerror");
+ to_png(dWerror,filebase,"dWerror");
+ to_png(ddAerror,filebase,"ddAerror");
+
+}
+
+int main(){
+ w_e();
+ return 0;
+}
+
Modified: trunk/ghost/monty/chirp/harness.c
===================================================================
--- trunk/ghost/monty/chirp/harness.c 2011-03-28 18:22:35 UTC (rev 17919)
+++ trunk/ghost/monty/chirp/harness.c 2011-04-12 06:45:19 UTC (rev 17920)
@@ -5,11 +5,13 @@
#include <string.h>
#include <cairo/cairo.h>
#include <squishyio/squishyio.h>
-#include "sinusoids.h"
#include "smallft.h"
#include "scales.h"
+#include "chirp.h"
+#include "lpc.h"
-static int BLOCKSIZE=1024;
+static int BLOCKSIZE=2048;
+static int MAX_CPC=200; /* maximum of 200 chirp candidates per channel */
static float MINdB=-100;
static float MAXdB=0;
@@ -17,40 +19,23 @@
#define AREAS 2
static float AREA_SIZE[AREAS]={1.,1.};
static int AREA_HZOOM[AREAS]={1,1};
-static int AREA_VZOOM[AREAS]={2,2};
+static int AREA_VZOOM[AREAS]={16,1};
static int AREA_FRAMEOUT_TRIGGER=1;
+static int WAVEFORM_AREA=-1;
static int SPECTROGRAM_AREA=0;
+static int CHIRPOGRAM_AREA=-1;
-static int SPECTRUM_AREA=-1;
+static int SPECTRUM_AREA=1;
static int SPECTRUM_GRIDFONTSIZE=8;
-static int ENVELOPE_AREA=-1;
+static int ENVELOPE_AREA=-11;
+static int TRACK_AREA=1;
+
static int pic_w,pic_h;
static drft_lookup fft;
-
-#define A0 .35875f
-#define A1 .48829f
-#define A2 .14128f
-#define A3 .01168f
-
-void spread_bark(float *in, float *out, int n, int rate,
- float loslope, float width, float hislope);
-
-
-void blackmann_harris(float *out, float *in, int n){
- int i;
- float scale = 2*M_PI/n;
-
- for(i=0;i<n;i++){
- float i5 = i+.5;
- float w = A0 - A1*cos(scale*i5) + A2*cos(scale*i5*2) - A3*cos(scale*i5*3);
- out[i] = in[i]*w;
- }
-}
-
static void hanning(float *out, float *in, int n){
int i;
float scale = 2*M_PI/n;
@@ -167,7 +152,9 @@
*r=1.;
break;
case 1:
- *g=1.;
+ *r=.7;
+ *g=.7;
+ *b=1;
break;
case 2:
*b=1.;
@@ -183,6 +170,7 @@
case 0:
*r=.7;
*b=.7;
+ *g=1;
break;
case 6:
*r=.45;
@@ -196,12 +184,14 @@
}
}
-void spectrogram(pcm_t *pcm, cairo_t *c,off_t pos){
+void spectrogram(pcm_t *pcm,
+ float **current_dB,
+ float **current_spread,
+ cairo_t *c,
+ off_t pos){
int i,j;
float r,g,b;
float work[BLOCKSIZE];
- float work2[BLOCKSIZE];
- float pc[BLOCKSIZE];
int y = area_y(SPECTROGRAM_AREA)+area_h(SPECTROGRAM_AREA)-1;
int hzoom = AREA_HZOOM[SPECTROGRAM_AREA];
cairo_save(c);
@@ -211,25 +201,8 @@
cairo_set_operator(c,CAIRO_OPERATOR_ADD);
for(i=0;i<pcm->ch;i++){
- if(pos+BLOCKSIZE<=pcm->samples)
- memcpy(work, &pcm->data[i][pos],BLOCKSIZE*sizeof(*work));
- else{
- memset(work, 0, sizeof(work));
- memcpy(work, &pcm->data[i][pos],(pcm->samples-pos)*sizeof(*work));
- }
-
- /* window and drft */
- hanning(work, work,BLOCKSIZE);
- drft_forward(&fft, work);
- for(j=0;j<BLOCKSIZE;j++)
- work[j]*=4./BLOCKSIZE;
-
- mag_dB(work,work,BLOCKSIZE);
- spread_bark(work, work2, BLOCKSIZE/2+1, pcm->rate, -20, 1., -20);
for(j=0;j<BLOCKSIZE/2;j++)
- work[j] = work[j]-work2[j]-50.;
- for(j=0;j<BLOCKSIZE/2;j++)
- pc[j] = pc[j]-work2[j]-50.;
+ work[j] = current_dB[i][j]-current_spread[i][j]-50.;
pick_a_color(i,&r,&g,&b);
for(j=0;j<BLOCKSIZE/2;j++){
@@ -245,17 +218,64 @@
cairo_restore(c);
}
+extern void slow_dcft_row(float *x, float *y, int k, int n);
+void chirpogram(pcm_t *pcm,
+ float *window,
+ float **current_block,
+ float **current_dB,
+ float **current_spread,
+ cairo_t *c,
+ off_t pos){
+ int i,j,k;
+ float r,g,b;
+ float work[BLOCKSIZE];
+ float work2[BLOCKSIZE];
+ int h = area_h(CHIRPOGRAM_AREA);
+ int y = area_y(CHIRPOGRAM_AREA);
+ int hzoom = AREA_HZOOM[CHIRPOGRAM_AREA];
+ cairo_save(c);
+ clip_area(c,CHIRPOGRAM_AREA);
+ clear_area(c,CHIRPOGRAM_AREA);
+
+ cairo_set_operator(c,CAIRO_OPERATOR_ADD);
+
+ for(i=0;i<pcm->ch;i++){
+ memcpy(work2,current_block[i],sizeof(work2));
+ for(j=0;j<BLOCKSIZE;j++)work2[j]*=window[j];
+
+ for(k=0;k<h;k++){
+ slow_dcft_row(work2,work, k, BLOCKSIZE);
+
+ for(j=0;j<BLOCKSIZE/2;j++)
+ work[j] = todB(work[j]*2/BLOCKSIZE);//-current_spread[i][j]-50.;
+
+ pick_a_color(i+1,&r,&g,&b);
+ for(j=0;j<BLOCKSIZE/2;j++){
+ float alpha = work[j]*(1./-MINdB)+1.f;
+ if(alpha<0.)alpha=0.;
+ if(alpha>1.)alpha=1.;
+ cairo_set_source_rgba(c,r,g,b,alpha);
+ cairo_rectangle(c,j*hzoom-hzoom/2,y+k,hzoom,1);
+ cairo_fill(c);
+ }
+ }
+ }
+
+ cairo_restore(c);
+}
+
+
void grid_lines(cairo_t *c,int area){
int y1 = area_y(area);
int h = area_h(area);
int y2 = y1+h;
int j;
-
+ int spacing = (AREA_HZOOM[area]<4 ? 10: AREA_HZOOM[area]);
cairo_save(c);
clip_area(c,area);
cairo_set_line_width(c,1.);
cairo_set_source_rgb(c,.2,.2,.2);
- for(j=0;j<pic_w;j+=AREA_HZOOM[area]){
+ for(j=0;j<pic_w;j+=spacing){
cairo_move_to(c,j+.5,y1);
cairo_line_to(c,j+.5,y2);
}
@@ -319,145 +339,67 @@
}
-void spectrum(pcm_t *pcm,cairo_t *c,off_t pos){
- float work[BLOCKSIZE];
- int i,j;
+void spectrum(pcm_t *pcm,float **current_dB,cairo_t *c,off_t pos){
+ int i;
cairo_save(c);
clip_area(c,SPECTRUM_AREA);
clear_area(c,SPECTRUM_AREA);
grid_lines(c,SPECTRUM_AREA);
cairo_set_line_width(c,.7);
- for(i=0;i<pcm->ch;i++){
- if(pos+BLOCKSIZE<=pcm->samples)
- memcpy(work, &pcm->data[i][pos],BLOCKSIZE*sizeof(*work));
- else{
- memset(work, 0, sizeof(work));
- memcpy(work, &pcm->data[i][pos],(pcm->samples-pos)*sizeof(*work));
- }
+ for(i=0;i<pcm->ch;i++)
+ frequency_vector_dB(current_dB[i],i,BLOCKSIZE/2+1,SPECTRUM_AREA,c,0);
- /* window and drft */
- hanning(work, work,BLOCKSIZE);
- drft_forward(&fft, work);
- for(j=0;j<BLOCKSIZE;j++)
- work[j]*=4./BLOCKSIZE;
-
- mag_dB(work,work,BLOCKSIZE);
- frequency_vector_dB(work,i,BLOCKSIZE/2+1,SPECTRUM_AREA,c,1);
- }
cairo_restore(c);
}
-void spread_bark_naive(float *in, float *out, int n, int rate,
- float loslope, float width, float hislope){
+void waveform_vector(float *vec,int ch,int n,
+ int area, cairo_t *c,int fillp){
+ float r,g,b;
+ int y1 = area_y(area);
+ int h = area_h(area);
+ float hzoom = AREA_HZOOM[area]*.5;
+ int j;
- static float *spread_integral=NULL;
- static float *bark=NULL;
- static int spread_integral_n=0;
-
- int i,j;
- memset(out,0,sizeof(*out)*n);
-
- if(spread_integral_n!=n){
- float lobarkwidth=toBark(.5/n*rate)-toBark(0);
- float hibarkwidth=toBark(.5*rate)-toBark(.5*(n-1)/n*rate);
-
- if(spread_integral)free(spread_integral);
- if(bark)free(bark);
- spread_integral = calloc(n,sizeof(*spread_integral));
- bark = calloc(n,sizeof(*bark));
- spread_integral_n = n;
-
- for(i=0;i<n;i++){
- float total=1.;
- float att=0.;
- float centerbark = toBark(.5*i/n*rate);
- j=i;
- bark[i]=centerbark;
-
- /* loop for computation that assumes we've not fallen off bottom
- edge of the spectrum */
- while(--j>=0){
- float b = toBark(.5*j/n*rate);
- if(b>centerbark-width*.5){
- /* still inside the F3 point */
- att = (centerbark-b)/(width*.5)*-3.;
- }else{
- /* outside the F3 point */
- att = ((centerbark-width*.5)-b)*loslope-3.;
- }
- total += fromdB(att);
- }
-
- /* avoid area contraction at edges of spectrum */
- while(att>-120.){
- att += (lobarkwidth*loslope);
- total += fromdB(att);
- }
-
- j=i;
- /* loop for computation that assumes we've not fallen off top
- edge of the spectrum */
- while(++j<n){
- float b = toBark(.5*j/n*rate);
- if(b<centerbark+width*.5){
- /* still inside the F3 point */
- att = (b-centerbark)/(width*.5)*-3.;
- }else{
- /* outside the F3 point */
- att = (b-(width*.5+centerbark))*hislope-3.;
- }
- total += fromdB(att);
- }
-
- /* avoid area contraction at edges of spectrum */
- while(att>-120.){
- att += (hibarkwidth*hislope);
- total += fromdB(att);
- }
-
- /* numerical integration complete to well beyond limits of reason */
- spread_integral[i]=total;
+ for(j=0;j<n;j++){
+ float y = (-vec[j]+1)*h*.5;
+ if(j==0){
+ cairo_move_to(c,j*hzoom+.5,y1+y+.5);
+ }else{
+ cairo_line_to(c,j*hzoom+.5,y1+y+.5);
}
}
+ pick_a_color(ch,&r,&g,&b);
+ if(fillp){
+ cairo_line_to(c,pic_w-.5,y1+h-.5);
+ cairo_line_to(c,.5,y1+h-.5);
+ cairo_close_path(c);
+ cairo_set_source_rgba(c,r,g,b,.2);
+ cairo_fill(c);
+ cairo_new_path(c);
+ }else{
+ cairo_set_source_rgb(c,r,g,b);
+ cairo_stroke(c);
+ }
- for(i=0;i<n;i++){
- float centerbark = bark[i];
- float e = fromdB(in[i])*fromdB(in[i])/spread_integral[i];
- j=i;
+}
- out[j] += e;
- while(--j>=0){
- float b = bark[j];
- float att;
- if(b>centerbark-width*.5){
- /* still inside the F3 point */
- att = (centerbark-b)/(width*.5)*-3.;
- }else{
- /* outside the F3 point */
- att = ((centerbark-width*.5)-b)*loslope-3.;
- }
- out[j] += fromdB(att)*e;
- }
+void waveform(pcm_t *pcm,float *window,float **current,cairo_t *c,off_t pos){
+ int i,j;
+ float work[BLOCKSIZE];
+ cairo_save(c);
+ clip_area(c,WAVEFORM_AREA);
+ clear_area(c,WAVEFORM_AREA);
+ cairo_set_line_width(c,1.);
- j=i;
- while(++j<n){
- float b = bark[j];
- float att;
- if(b<centerbark+width*.5){
- /* still inside the F3 point */
- att = (b-centerbark)/(width*.5)*-3.;
- }else{
- /* outside the F3 point */
- att = (b-(width*.5+centerbark))*hislope-3.;
- }
- out[j] += fromdB(att)*e;
- }
+ for(i=0;i<pcm->ch;i++){
+ for(j=0;j<BLOCKSIZE;j++)
+ work[j]= -current[i][j]*window[j];
+ waveform_vector(work,i,BLOCKSIZE,WAVEFORM_AREA,c,0);
+ waveform_vector(window,i,BLOCKSIZE,WAVEFORM_AREA,c,0);
}
- for(i=0;i<n;i++)
- out[i] = todB(sqrt(out[i]));
-
+ cairo_restore(c);
}
void spread_bark(float *in, float *out, int n, int rate,
@@ -630,48 +572,284 @@
dBlo *= fromdB(hislope*barkw[i]);
}
-
-
-
for(i=0;i<n;i++)
out[i]=todB(sqrt(out[i]));
}
-void envelope(pcm_t *pcm,cairo_t *c,off_t pos){
- float work[BLOCKSIZE];
- float work2[BLOCKSIZE];
- int i,j;
+void envelope(pcm_t *pcm,float **current_spread,cairo_t *c,off_t pos){
+ int i;
cairo_save(c);
clip_area(c,SPECTRUM_AREA);
cairo_set_line_width(c,.7);
+ for(i=0;i<pcm->ch;i++)
+ frequency_vector_dB(current_spread[i],i+pcm->ch,BLOCKSIZE/2+1,ENVELOPE_AREA,c,0);
+
+ cairo_restore(c);
+}
+
+static int seed_chirps(float *fft,float *dB,float *spread,
+ chirp *chirps,int *n){
+ int i;
+ int ret=0;
+
+#if 0
+ if(dB[0]>spread[0]+5 &&
+ dB[0]>dB[1] &&
+ dB[0]>-100){
+ if(*n<MAX_CPC){
+ chirps[*n].A=fft[0];
+ chirps[*n].W=.000001;
+ chirps[*n].P=0;
+ chirps[*n].dA=0;
+ chirps[*n].dW=0;
+ chirps[*n].label=-1;
+ ret=1;
+ (*n)++;
+ }
+ }
+
+ if(dB[BLOCKSIZE/2]>spread[BLOCKSIZE/2]+5 &&
+ dB[BLOCKSIZE/2]>dB[BLOCKSIZE/2-1]&&
+ dB[BLOCKSIZE/2]>-100){
+ if(*n<MAX_CPC){
+ chirps[*n].A=fft[BLOCKSIZE-1];
+ chirps[*n].W=M_PI;
+ chirps[*n].P=M_PI*.5;
+ chirps[*n].dA=0;
+ chirps[*n].dW=0;
+ chirps[*n].label=-1;
+ ret=1;
+ (*n)++;
+ }
+ }
+#endif
+ for(i=3;i<BLOCKSIZE/2;i++){
+ if(dB[i]>spread[i]+5 &&
+ dB[i]>dB[i-1] &&
+ dB[i]>dB[i+1] &&
+ dB[i]>-100){
+ if(*n<MAX_CPC){
+ chirps[*n].A=fromdB(dB[i]);
+ chirps[*n].W=(float)i/BLOCKSIZE*2.*M_PI;
+ chirps[*n].P=-atan2f(fft[i*2],fft[i*2-1]);
+ chirps[*n].dA=0;
+ chirps[*n].dW=0;
+ chirps[*n].label=-1;
+ ret=1;
+ (*n)++;
+ }else
+ break;
+ }
+ }
+
+ return ret;
+}
+
+void dump_waveform(char *base,int no,float *d,int n){
+ char buf[80];
+ FILE *f;
+ snprintf(buf,80,"%s_%d.m",base,no);
+ f=fopen(buf,"w");
+ if(f){
+ int i;
+ for(i=0;i<n;i++)
+ fprintf(f,"%f\n",d[i]);
+ fclose(f);
+ }
+}
+
+
+static int fitno=0;
+void track(pcm_t *pcm,
+ float *window,
+ float **current_block,
+ float **current_fft,
+ float **current_dB,
+ float **current_spread,
+ chirp **chirps,
+ int *chirps_used,
+ cairo_t *c,
+ off_t off,
+ float **y,
+ float **r){
+ int i,j,k;
+
+ static int count=0;
+ fprintf(stderr,"pos: %ld\n",(long)off);
+ float rec[BLOCKSIZE];
+ for(i=0;i<BLOCKSIZE;i++)rec[i]=1.;
+
+
for(i=0;i<pcm->ch;i++){
- if(pos+BLOCKSIZE<=pcm->samples)
- memcpy(work, &pcm->data[i][pos],BLOCKSIZE*sizeof(*work));
- else{
- memset(work, 0, sizeof(work));
- memcpy(work, &pcm->data[i][pos],(pcm->samples-pos)*sizeof(*work));
+ float lfft[BLOCKSIZE];
+ float ldB[BLOCKSIZE/2+1];
+ float lspread[BLOCKSIZE/2+1];
+ chirp prev[MAX_CPC];
+ int prev_chirps;
+ int flag=1;
+ int ret;
+ memcpy(prev,chirps[i],chirps_used[i]*sizeof(**chirps));
+ prev_chirps = chirps_used[i];
+
+ /* advance preexisting peaks */
+ advance_chirps(chirps[i], chirps_used[i], BLOCKSIZE/AREA_VZOOM[TRACK_AREA]);
+ for(j=0;j<chirps_used[i];j++)chirps[i][j].label=j;
+
+ fprintf(stderr,"FITNO: %d\n",fitno);
+ fprintf(stderr,"chirps: %d\n",chirps_used[i]);
+ for(j=0;j<chirps_used[i];j++)
+ fprintf(stderr,"in>>%f:%f:%f::%f:%f ",chirps[i][j].A,chirps[i][j].W*BLOCKSIZE/2./M_PI,chirps[i][j].P,chirps[i][j].dA*BLOCKSIZE,chirps[i][j].dW*BLOCKSIZE*BLOCKSIZE/2/M_PI);
+
+ ret=estimate_chirps(current_block[i],y[i],r[i],window,BLOCKSIZE,
+ chirps[i],chirps_used[i],50,.01);
+ for(j=0;j<chirps_used[i];j++)
+ fprintf(stderr,"out(%d)>>%f:%f:%f::%f:%f ",ret,chirps[i][j].A,chirps[i][j].W*BLOCKSIZE/2./M_PI,chirps[i][j].P,chirps[i][j].dA*BLOCKSIZE,chirps[i][j].dW*BLOCKSIZE*BLOCKSIZE/2./M_PI);
+ dump_waveform("x",fitno,current_block[i],BLOCKSIZE);
+ dump_waveform("y",fitno,y[i],BLOCKSIZE);
+ dump_waveform("r",fitno,r[i],BLOCKSIZE);
+ fitno++;
+ fprintf(stderr,"\n");
+
+ /* perform spread on residue */
+ memcpy(lfft, r[i], BLOCKSIZE*sizeof(*lfft));
+ drft_forward(&fft, lfft);
+ for(j=0;j<BLOCKSIZE;j++) lfft[j]*=4./BLOCKSIZE;
+ mag_dB(ldB,lfft,BLOCKSIZE);
+ spread_bark(ldB, lspread, BLOCKSIZE/2+1,
+ pcm->rate,-20., 1., -20.);
+
+ /* 5?: scan residue for peaks exceeding theshold; add to end of
+ chirps array if no advanced chirp is within a bin */
+ flag=seed_chirps(lfft,ldB,current_spread[i],
+ chirps[i],chirps_used+i);
+
+ fprintf(stderr,"FITNO: %d\n",fitno);
+ fprintf(stderr,"chirps: %d\n",chirps_used[i]);
+ for(j=0;j<chirps_used[i];j++)
+ fprintf(stderr,"in>>%f:%f:%f::%f:%f ",chirps[i][j].A,chirps[i][j].W*BLOCKSIZE/2./M_PI,chirps[i][j].P,chirps[i][j].dA*BLOCKSIZE,chirps[i][j].dW*BLOCKSIZE*BLOCKSIZE/2/M_PI);
+
+ ret=estimate_chirps(current_block[i],y[i],r[i],(count>10?rec:window),BLOCKSIZE,
+ chirps[i],chirps_used[i],50,.01);
+
+ for(j=0;j<chirps_used[i];j++)
+ fprintf(stderr,"out(%d)>>%f:%f:%f::%f:%f ",ret,chirps[i][j].A,chirps[i][j].W*BLOCKSIZE/2./M_PI,chirps[i][j].P,chirps[i][j].dA*BLOCKSIZE,chirps[i][j].dW*BLOCKSIZE*BLOCKSIZE/2./M_PI);
+ dump_waveform("x",fitno,current_block[i],BLOCKSIZE);
+ dump_waveform("y",fitno,y[i],BLOCKSIZE);
+ dump_waveform("r",fitno,r[i],BLOCKSIZE);
+ fitno++;
+ fprintf(stderr,"\n");
+
+#if 0
+ for(k=0;k<2;k++){
+ /* fit/track */
+ estimate_chirps(current_block[i],y[i],r[i],window,BLOCKSIZE,
+ chirps[i],chirps_used[i],50,.01);
+
+ /* perform spread on residue */
+ memcpy(lfft, r[i], BLOCKSIZE*sizeof(*lfft));
+ drft_forward(&fft, lfft);
+ for(j=0;j<BLOCKSIZE;j++) lfft[j]*=4./BLOCKSIZE;
+ mag_dB(ldB,lfft,BLOCKSIZE);
+ //spread_bark(ldB, lspread, BLOCKSIZE/2+1,
+ // pcm->rate,-40., 1., -40.);
+
+ /* 5?: scan residue for peaks exceeding theshold; add to end of
+ chirps array if no advanced chirp is within a bin */
+ flag=seed_chirps(lfft,ldB,current_spread[i],
+ chirps[i],chirps_used+i);
+
+ fprintf(stderr,"chirps: %d\n",chirps_used[i]);
+ /* 6?: if residue peaks, loop */
}
+#endif
+ /* 7: compare all chirps against original spread; cull those
+ below threshold, also remove them from reconstruction */
+ //cull_chirps(current_fft[i],current_dB[i],current_spread[i],
+ // chirps[i],chirps_used+i);
- /* window and drft */
- hanning(work, work,BLOCKSIZE);
- drft_forward(&fft, work);
- for(j=0;j<BLOCKSIZE;j++)
- work[j]*=4./BLOCKSIZE;
+ /* 8: draw telemetry */
+ /* The spectrum gets lollypops */
+ {
+ float r,g,b;
+ int ytop = area_y(SPECTRUM_AREA);
+ int h = area_h(SPECTRUM_AREA);
+ int hzoom = AREA_HZOOM[SPECTRUM_AREA];
+ int xw = (hzoom<4?4:hzoom);
+ cairo_save(c);
+ clip_area(c,SPECTRUM_AREA);
+ cairo_set_line_width(c,.7);
- mag_dB(work,work,BLOCKSIZE);
- spread_bark(work, work2, BLOCKSIZE/2+1, pcm->rate,
- -20., 1., -20.);
- frequency_vector_dB(work2,i+pcm->ch,BLOCKSIZE/2+1,ENVELOPE_AREA,c,0);
+ for(j=0;j<chirps_used[i];j++){
+ float x = chirps[i][j].W*BLOCKSIZE*hzoom/2./M_PI;
+ float y1 = todB(chirps[i][j].A - chirps[i][j].dA*.5)*(1./MINdB);
+ float y2 = todB(chirps[i][j].A + chirps[i][j].dA*.5)*(1./MINdB);
+ cairo_move_to(c,x+.5,ytop+y1*h+.5);
+ cairo_line_to(c,x+.5,ytop+h-.5);
+
+ cairo_move_to(c,x+.5-xw,ytop+y1*h+.5);
+ cairo_line_to(c,x+.5+xw,ytop+y1*h+.5);
+ cairo_move_to(c,x+.5-xw,ytop+y2*h+.5);
+ cairo_line_to(c,x+.5+xw,ytop+y2*h+.5);
+
+ }
+ cairo_set_source_rgb(c,0.,1.,1.);
+ //pick_a_color(i,&r,&g,&b);
+ cairo_stroke(c);
+ cairo_restore(c);
+ }
+
+ /* the spectrogram gets tracer lines */
+ {
+ float r,g,b;
+ int ytop = area_y(SPECTROGRAM_AREA);
+ int h = area_h(SPECTROGRAM_AREA);
+ int hzoom = AREA_HZOOM[SPECTROGRAM_AREA];
+ int vzoom = AREA_VZOOM[SPECTROGRAM_AREA]/AREA_VZOOM[SPECTROGRAM_AREA];
+ cairo_save(c);
+ clip_area(c,SPECTROGRAM_AREA);
+ cairo_set_line_width(c,.7);
+
+ for(j=0;j<chirps_used[i];j++){
+ if(chirps[i][j].label>=0){
+ float x1 = prev[chirps[i][j].label].W*BLOCKSIZE*hzoom/2./M_PI;
+ float x2 = chirps[i][j].W*BLOCKSIZE*hzoom/2./M_PI;
+ cairo_move_to(c,x1+.5,ytop+h-vzoom);
+ cairo_line_to(c,x2+.5,ytop+h);
+ }else{
+ float x1 = (chirps[i][j].W-chirps[i][j].dW*.5)*BLOCKSIZE*hzoom/2./M_PI;
+ float x2 = chirps[i][j].W*BLOCKSIZE*hzoom/2./M_PI;
+ cairo_move_to(c,x1+.5,ytop+h-vzoom);
+ cairo_line_to(c,x2+.5,ytop+h);
+ }
+ }
+ cairo_set_source_rgb(c,1.,1.,1.);
+ //pick_a_color(i,&r,&g,&b);
+ cairo_stroke(c);
+ cairo_restore(c);
+ }
}
- cairo_restore(c);
+ count++;
}
+
int main(int argc, char **argv){
+ float **current_block;
+ float **current_fft;
+ float **current_dB;
+ float **current_spread;
+ float **current_reconstruct;
+ float **current_residue;
+ float window[BLOCKSIZE];
+
+ int i,j;
+
pcm_t *pcm;
cairo_t *c;
+ chirp **chirps=NULL;
+ int *chirps_used;
int areas_active[AREAS]={1,1};
off_t areas_pos[AREAS]={0,0};
@@ -679,12 +857,42 @@
pic_w = BLOCKSIZE/2;
pic_h = pic_w*9/16;
-
drft_init(&fft, BLOCKSIZE);
/* squishyload a file. */
pcm = squishyio_load_file(argv[1]);
+ current_block = calloc(pcm->ch, sizeof(*current_block));
+ current_fft = calloc(pcm->ch, sizeof(*current_fft));
+ current_dB = calloc(pcm->ch, sizeof(*current_dB));
+ current_spread = calloc(pcm->ch, sizeof(*current_spread));
+ current_residue = calloc(pcm->ch, sizeof(*current_residue));
+ current_reconstruct = calloc(pcm->ch, sizeof(*current_reconstruct));
+ for(i=0;i<pcm->ch;i++){
+#if 0
+ pcm->data[i]=realloc(pcm->data[i],
+ (pcm->samples+BLOCKSIZE*2)*sizeof(**pcm->data));
+ memmove(pcm->data[i]+BLOCKSIZE,pcm->data[i],
+ sizeof(**pcm->data)*pcm->samples);
+
+ preextrapolate(pcm->data[i]+BLOCKSIZE,BLOCKSIZE,
+ pcm->data[i],BLOCKSIZE);
+ postextrapolate(pcm->data[i]+pcm->samples,BLOCKSIZE,
+ pcm->data[i]+pcm->samples+BLOCKSIZE,BLOCKSIZE);
+#endif
+ current_block[i] = calloc(BLOCKSIZE, sizeof(**current_block));
+ current_fft[i] = calloc(BLOCKSIZE, sizeof(**current_fft));
+ current_dB[i] = calloc(BLOCKSIZE/2+1, sizeof(**current_dB));
+ current_spread[i] = calloc(BLOCKSIZE/2+1, sizeof(**current_spread));
+ current_residue[i] = calloc(BLOCKSIZE, sizeof(**current_residue));
+ current_reconstruct[i] = calloc(BLOCKSIZE, sizeof(**current_reconstruct));
+
+ }
+#if 0
+ pcm->samples+=BLOCKSIZE*2-1;
+#endif
+ pcm->samples=(pcm->samples/BLOCKSIZE)*BLOCKSIZE;
+
/* Make a cairo drawable */
cairo_surface_t *pane=cairo_image_surface_create(CAIRO_FORMAT_RGB24,pic_w, pic_h);
if(!pane || cairo_surface_status(pane)!=CAIRO_STATUS_SUCCESS){
@@ -703,8 +911,15 @@
pic_w,pic_h,pcm->rate,BLOCKSIZE/AREA_VZOOM[AREA_FRAMEOUT_TRIGGER]);
printf("AUDIO R%d C%d\n",pcm->rate,pcm->ch);
+ chirps=calloc(pcm->ch,sizeof(*chirps));
+ chirps_used=calloc(pcm->ch,sizeof(*chirps_used));
+ for(i=0;i<pcm->ch;i++)
+ chirps[i]=calloc(MAX_CPC,sizeof(**chirps));
+
+ for(j=0;j<BLOCKSIZE;j++)window[j]=1.;
+ hanning(window,window,BLOCKSIZE);
+
while(1){
- int i;
/* which area is next? */
off_t min=0;
@@ -716,17 +931,41 @@
min_area=i;
}
- if(min>=pcm->samples) break;
+ if(min+BLOCKSIZE>pcm->samples) break;
+ /* 'pos' is the center of the current window of BLOCKSIZE */
+ for(i=0;i<pcm->ch;i++){
+ memcpy(current_block[i], &pcm->data[i][min],
+ BLOCKSIZE*sizeof(**current_block));
+
+ /* window and drft */
+ for(j=0;j<BLOCKSIZE;j++)current_fft[i][j]=window[j]*current_block[i][j];
+ drft_forward(&fft, current_fft[i]);
+ for(j=0;j<BLOCKSIZE;j++) current_fft[i][j]*=4./BLOCKSIZE;
+ mag_dB(current_dB[i],current_fft[i],BLOCKSIZE);
+ spread_bark(current_dB[i], current_spread[i], BLOCKSIZE/2+1,
+ pcm->rate,-20., 1., -20.);
+ }
+
if(min_area == SPECTROGRAM_AREA)
- spectrogram(pcm,c,areas_pos[min_area]);
+ spectrogram(pcm,current_dB,current_spread,c,min);
+ if(min_area == CHIRPOGRAM_AREA)
+ chirpogram(pcm,window,current_block,current_dB,current_spread,c,min);
+
if(min_area == SPECTRUM_AREA)
- spectrum(pcm,c,areas_pos[min_area]);
+ spectrum(pcm,current_dB,c,min);
+ if(min_area == WAVEFORM_AREA)
+ waveform(pcm,window,current_block,c,min);
+
if(min_area == ENVELOPE_AREA)
- envelope(pcm,c,areas_pos[min_area]);
+ envelope(pcm,current_spread,c,min);
+ if(min_area == TRACK_AREA)
+ track(pcm,window,current_block,current_fft,current_dB,current_spread,
+ chirps,chirps_used,c,min,
+ current_reconstruct,current_residue);
/* lastly */
if(min_area == AREA_FRAMEOUT_TRIGGER){
More information about the commits
mailing list