[xiph-commits] r13370 - trunk/ghost/libghost
jm at svn.xiph.org
jm at svn.xiph.org
Thu Jul 26 20:33:55 PDT 2007
Author: jm
Date: 2007-07-26 20:33:55 -0700 (Thu, 26 Jul 2007)
New Revision: 13370
Modified:
trunk/ghost/libghost/sinusoids.c
Log:
proper indentation (no code change)
Modified: trunk/ghost/libghost/sinusoids.c
===================================================================
--- trunk/ghost/libghost/sinusoids.c 2007-07-27 03:33:49 UTC (rev 13369)
+++ trunk/ghost/libghost/sinusoids.c 2007-07-27 03:33:55 UTC (rev 13370)
@@ -264,160 +264,160 @@
/* Iterate on the non-linear part (frequency) a couple times */
while (1)
{
- /* 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(w[i]);
- rotI = sin(w[i]);
- /* Computing sin/cos table using complex rotations.
- Only need to compute the tables for half the length because of the symmetry.*/
- cos_table[i][0] = cos(.5*w[i]);
- sin_table[i][0] = sin(.5*w[i]);
- for (j=1;j<L2;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 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;
+ float tmp1=0, tmp2=0;
+ float tmp3=0, tmp4=0;
+ float rotR, rotI;
+ rotR = cos(w[i]);
+ rotI = sin(w[i]);
+ /* Computing sin/cos table using complex rotations.
+ Only need to compute the tables for half the length because of the symmetry.*/
+ cos_table[i][0] = cos(.5*w[i]);
+ sin_table[i][0] = sin(.5*w[i]);
+ 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;
+ }
+ /* Compute energy of basis functions. There may be a closed-form solution for that
+ or we could tabulate it. */
+ 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);
+ cosE_1[i] = 1.f/cosE[i];
+ sinE_1[i] = 1.f/sinE[i];
+ costE_1[i] = 1.f/costE[i];
+ sintE_1[i] = 1.f/sintE[i];
+ /* Normalise the basis (should multiply by the inverse instead) */
+ for (j=0;j<L2;j++)
+ {
+ cos_table[i][j] *= cosE_1[i];
+ sin_table[i][j] *= sinE_1[i];
+ tcos_table[i][j] *= costE_1[i];
+ tsin_table[i][j] *= sintE_1[i];
+ }
}
- /* Compute energy of basis functions. There may be a closed-form solution for that
- or we could tabulate it. */
+
+ /* 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++)
{
- 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];
+ sym[j] = x[j+L2]+x[L2-j-1];
+ anti[j] = x[j+L2]-x[L2-j-1];
}
- /* 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);
- cosE_1[i] = 1.f/cosE[i];
- sinE_1[i] = 1.f/sinE[i];
- costE_1[i] = 1.f/costE[i];
- sintE_1[i] = 1.f/sintE[i];
- /* Normalise the basis (should multiply by the inverse instead) */
- for (j=0;j<L2;j++)
+
+ /* Start the linear optimisation where we left last time */
+ for (i=0;i<N;i++)
{
- cos_table[i][j] *= cosE_1[i];
- sin_table[i][j] *= sinE_1[i];
- tcos_table[i][j] *= costE_1[i];
- tsin_table[i][j] *= sintE_1[i];
+ ai[i] = A[i]*cos(phi[i])*cosE[i];
+ bi[i] = -A[i]*sin(phi[i])*sinE[i];
+ ci[i] = dA[i]*cos(phi[i])*costE[i];
+ di[i] = -dA[i]*sin(phi[i])*sintE[i];
+ for (j=0;j<L2;j++)
+ sym[j] -= (2*ai[i])*cos_table[i][j] + (2*di[i])*tsin_table[i][j];
+ for (j=0;j<L2;j++)
+ anti[j] -= (2*bi[i])*sin_table[i][j] + (2*ci[i])*tcos_table[i][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] = x[j+L2]+x[L2-j-1];
- anti[j] = x[j+L2]-x[L2-j-1];
- }
-
- /* Start the linear optimisation where we left last time */
- for (i=0;i<N;i++)
- {
- ai[i] = A[i]*cos(phi[i])*cosE[i];
- bi[i] = -A[i]*sin(phi[i])*sinE[i];
- ci[i] = dA[i]*cos(phi[i])*costE[i];
- di[i] = -dA[i]*sin(phi[i])*sintE[i];
- for (j=0;j<L2;j++)
- sym[j] -= (2*ai[i])*cos_table[i][j] + (2*di[i])*tsin_table[i][j];
- for (j=0;j<L2;j++)
- anti[j] -= (2*bi[i])*sin_table[i][j] + (2*ci[i])*tcos_table[i][j];
- }
- /* Stop condition. We stop here, just after updating the error using the new frequency */
- if (nonlin-- <= 0)
- break;
+ /* Stop condition. We stop here, just after updating the error using the new frequency */
+ if (nonlin-- <= 0)
+ break;
- /* This is an iterative linear solution -- much quicker than inverting a matrix.
- It's pretty much the Gauss-Seidel algorithm */
- for (iter=0;iter<2;iter++)
- {
- for (i=0;i<N;i++)
+ /* This is an iterative linear solution -- much quicker than inverting a matrix.
+ It's pretty much the Gauss-Seidel algorithm */
+ for (iter=0;iter<2;iter++)
{
- float tmp1=0, tmp2=0;
- float tmp3=0, tmp4=0;
- /* For each of the four basis functions, project the residual (symmetric or
- anti-symmetric) onto the basis function, then update the residual. */
- for (j=0;j<L2;j++)
- tmp1 += sym[j]*cos_table[i][j];
- for (j=0;j<L2;j++)
- sym[j] -= (2*tmp1)*cos_table[i][j];
+ for (i=0;i<N;i++)
+ {
+ float tmp1=0, tmp2=0;
+ float tmp3=0, tmp4=0;
+ /* For each of the four basis functions, project the residual (symmetric or
+ anti-symmetric) onto the basis function, then update the residual. */
+ for (j=0;j<L2;j++)
+ tmp1 += sym[j]*cos_table[i][j];
+ for (j=0;j<L2;j++)
+ sym[j] -= (2*tmp1)*cos_table[i][j];
- for (j=0;j<L2;j++)
- tmp2 += anti[j]*sin_table[i][j];
- for (j=0;j<L2;j++)
- anti[j] -= (2*tmp2)*sin_table[i][j];
+ for (j=0;j<L2;j++)
+ tmp2 += anti[j]*sin_table[i][j];
+ for (j=0;j<L2;j++)
+ anti[j] -= (2*tmp2)*sin_table[i][j];
- for (j=0;j<L2;j++)
- tmp3 += anti[j]*tcos_table[i][j];
- for (j=0;j<L2;j++)
- anti[j] -= (2*tmp3)*tcos_table[i][j];
+ for (j=0;j<L2;j++)
+ tmp3 += anti[j]*tcos_table[i][j];
+ for (j=0;j<L2;j++)
+ anti[j] -= (2*tmp3)*tcos_table[i][j];
- for (j=0;j<L2;j++)
- tmp4 += sym[j]*tsin_table[i][j];
- for (j=0;j<L2;j++)
- sym[j] -= (2*tmp4)*tsin_table[i][j];
+ for (j=0;j<L2;j++)
+ tmp4 += sym[j]*tsin_table[i][j];
+ for (j=0;j<L2;j++)
+ sym[j] -= (2*tmp4)*tsin_table[i][j];
- ai[i] += tmp1;
- bi[i] += tmp2;
- ci[i] += tmp3;
- di[i] += tmp4;
+ ai[i] += tmp1;
+ bi[i] += tmp2;
+ ci[i] += tmp3;
+ di[i] += tmp4;
+ }
}
- }
- for (i=0;i<N;i++)
- {
- ai[i] *= cosE_1[i];
- bi[i] *= sinE_1[i];
- ci[i] *= costE_1[i];
- di[i] *= sintE_1[i];
- }
+ for (i=0;i<N;i++)
+ {
+ ai[i] *= cosE_1[i];
+ bi[i] *= sinE_1[i];
+ ci[i] *= costE_1[i];
+ di[i] *= sintE_1[i];
+ }
- /* Update real sinusoid parameters */
- for (i=0;i<N;i++)
- {
- //float A, phi, dA, dw;
- A[i] = sqrt(ai[i]*ai[i] + bi[i]*bi[i]);
- if (A[i]<.01)
- phi[i] = 0;
- else
- phi[i] = atan2(-bi[i], ai[i]);
- dA[i] = (ci[i]*ai[i] + bi[i]*di[i])/(.01f+A[i]);
- dw[i] = (ci[i]*bi[i] - di[i]*ai[i])/(.01f+A[i]*A[i]);
- //printf ("%f %f %f %f %f %f %f %f %f\n", w[i], ai[i], bi[i], ci[i], di[i], A, phi, dA, dw);
- /* Update the frequency with only a fraction of the dw term to prevent instabilities */
- w[i] = w[i] + .5*dw[i];
- /* Make sure we didn't deviate too much from the original frequency */
- if (w[i] > w0[i]+.025)
- w[i] = w0[i]+.025;
- if (w[i] < w0[i]-.025)
- w[i] = w0[i]-.025;
- //Prevent stupid frequencies (i.e. outsite of ]0,pi[
- if (w[i] < .002)
- w[i] = .002;
- if (w[i] > M_PI-.01)
- w[i] = M_PI-.01;
- //if (nonlin == 4)
- // 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]);
+ /* Update real sinusoid parameters */
+ for (i=0;i<N;i++)
+ {
+ //float A, phi, dA, dw;
+ A[i] = sqrt(ai[i]*ai[i] + bi[i]*bi[i]);
+ if (A[i]<.01)
+ phi[i] = 0;
+ else
+ phi[i] = atan2(-bi[i], ai[i]);
+ dA[i] = (ci[i]*ai[i] + bi[i]*di[i])/(.01f+A[i]);
+ dw[i] = (ci[i]*bi[i] - di[i]*ai[i])/(.01f+A[i]*A[i]);
+ //printf ("%f %f %f %f %f %f %f %f %f\n", w[i], ai[i], bi[i], ci[i], di[i], A, phi, dA, dw);
+ /* Update the frequency with only a fraction of the dw term to prevent instabilities */
+ w[i] = w[i] + .5*dw[i];
+ /* Make sure we didn't deviate too much from the original frequency */
+ if (w[i] > w0[i]+.025)
+ w[i] = w0[i]+.025;
+ if (w[i] < w0[i]-.025)
+ w[i] = w0[i]-.025;
+ //Prevent stupid frequencies (i.e. outsite of ]0,pi[
+ if (w[i] < .002)
+ w[i] = .002;
+ if (w[i] > M_PI-.01)
+ w[i] = M_PI-.01;
+ //if (nonlin == 4)
+ // 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]);
+ }
}
- }
for (j=0;j<L2;j++)
{
More information about the commits
mailing list