celt/bands.c

@@ -41,6 +41,21 @@
 #include "mathops.h"
 #include "rate.h"
 
+int hysteresis_decision(opus_val16 val, const opus_val16 *thresholds, const opus_val16 *hysteresis, int N, int prev)
+{
+   int i;
+   for (i=0;i<N;i++)
+   {
+      if (val < thresholds[i])
+         break;
+   }
+   if (i>prev && val < thresholds[prev]+hysteresis[prev])
+      i=prev;
+   if (i<prev && val > thresholds[prev-1]-hysteresis[prev-1])
+      i=prev;
+   return i;
+}
+
 opus_uint32 celt_lcg_rand(opus_uint32 seed)
 {
    return 1664525 * seed + 1013904223;

celt/bands.h

@@ -92,4 +92,6 @@ void anti_collapse(const CELTMode *m, celt_norm *X_, unsigned char *collapse_mas
 
 opus_uint32 celt_lcg_rand(opus_uint32 seed);
 
+int hysteresis_decision(opus_val16 val, const opus_val16 *thresholds, const opus_val16 *hysteresis, int N, int prev);
+
 #endif /* BANDS_H */

celt/celt.h

@@ -50,7 +50,20 @@ extern "C" {
 #define CELTDecoder OpusCustomDecoder
 #define CELTMode OpusCustomMode
 
-#define _celt_check_mode_ptr_ptr(ptr) ((ptr) + ((ptr) - (const CELTMode**)(ptr)))
+typedef struct {
+   int valid;
+   opus_val16 tonality;
+   opus_val16 tonality_slope;
+   opus_val16 noisiness;
+   opus_val16 activity;
+   int boost_band[2];
+   opus_val16 boost_amount[2];
+   opus_val16 music_prob;
+}AnalysisInfo;
+
+#define __celt_check_mode_ptr_ptr(ptr) ((ptr) + ((ptr) - (const CELTMode**)(ptr)))
+
+#define __celt_check_analysis_ptr(ptr) ((ptr) + ((ptr) - (const AnalysisInfo*)(ptr)))
 
 /* Encoder/decoder Requests */
 
@@ -81,11 +94,18 @@ extern "C" {
 
 #define CELT_GET_MODE_REQUEST    10015
 /** Get the CELTMode used by an encoder or decoder */
-#define CELT_GET_MODE(x) CELT_GET_MODE_REQUEST, _celt_check_mode_ptr_ptr(x)
+#define CELT_GET_MODE(x) CELT_GET_MODE_REQUEST, __celt_check_mode_ptr_ptr(x)
 
 #define CELT_SET_SIGNALLING_REQUEST    10016
 #define CELT_SET_SIGNALLING(x) CELT_SET_SIGNALLING_REQUEST, __opus_check_int(x)
 
+#define CELT_SET_TONALITY_REQUEST    10018
+#define CELT_SET_TONALITY(x) CELT_SET_TONALITY_REQUEST, __opus_check_int(x)
+#define CELT_SET_TONALITY_SLOPE_REQUEST    10020
+#define CELT_SET_TONALITY_SLOPE(x) CELT_SET_TONALITY_SLOPE_REQUEST, __opus_check_int(x)
+
+#define CELT_SET_ANALYSIS_REQUEST    10022
+#define CELT_SET_ANALYSIS(x) CELT_SET_ANALYSIS_REQUEST, __celt_check_analysis_ptr(x)
 
 
 /* Encoder stuff */

celt/mdct.c

@@ -109,12 +109,14 @@ void clt_mdct_forward(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_scalar
    int N, N2, N4;
    kiss_twiddle_scalar sine;
    VARDECL(kiss_fft_scalar, f);
+   VARDECL(kiss_fft_scalar, f2);
    SAVE_STACK;
    N = l->n;
    N >>= shift;
    N2 = N>>1;
    N4 = N>>2;
    ALLOC(f, N2, kiss_fft_scalar);
+   ALLOC(f2, N2, kiss_fft_scalar);
    /* sin(x) ~= x here */
 #ifdef FIXED_POINT
    sine = TRIG_UPSCALE*(QCONST16(0.7853981f, 15)+N2)/N;
@@ -180,12 +182,12 @@ void clt_mdct_forward(const mdct_lookup *l, kiss_fft_scalar *in, kiss_fft_scalar
    }
 
    /* N/4 complex FFT, down-scales by 4/N */
-   opus_fft(l->kfft[shift], (kiss_fft_cpx *)f, (kiss_fft_cpx *)in);
+   opus_fft(l->kfft[shift], (kiss_fft_cpx *)f, (kiss_fft_cpx *)f2);
 
    /* Post-rotate */
    {
       /* Temp pointers to make it really clear to the compiler what we're doing */
-      const kiss_fft_scalar * restrict fp = in;
+      const kiss_fft_scalar * restrict fp = f2;
       kiss_fft_scalar * restrict yp1 = out;
       kiss_fft_scalar * restrict yp2 = out+stride*(N2-1);
       const kiss_twiddle_scalar *t = &l->trig[0];

celt/pitch.c

@@ -321,6 +321,7 @@ opus_val16 remove_doubling(opus_val16 *x, int maxperiod, int minperiod,
       int T1, T1b;
       opus_val16 g1;
       opus_val16 cont=0;
+      opus_val16 thresh;
       T1 = (2*T0+k)/(2*k);
       if (T1 < minperiod)
          break;
@@ -362,7 +363,14 @@ opus_val16 remove_doubling(opus_val16 *x, int maxperiod, int minperiod,
          cont = HALF32(prev_gain);
       else
          cont = 0;
-      if (g1 > QCONST16(.3f,15) + MULT16_16_Q15(QCONST16(.4f,15),g0)-cont)
+      thresh = MAX16(QCONST16(.3f,15), MULT16_16_Q15(QCONST16(.7,15),g0)-cont);
+      /* Bias against very high pitch (very short period) to avoid false-positives
+         due to short-term correlation */
+      if (T1<3*minperiod)
+         thresh = MAX16(QCONST16(.4f,15), MULT16_16_Q15(QCONST16(.85,15),g0)-cont);
+      else if (T1<2*minperiod)
+         thresh = MAX16(QCONST16(.5f,15), MULT16_16_Q15(QCONST16(.9,15),g0)-cont);
+      if (g1 > thresh)
       {
          best_xy = xy;
          best_yy = yy;

src/analysis.c

+/* Copyright (c) 2011 Xiph.Org Foundation
+   Written by Jean-Marc Valin */
+/*
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions
+   are met:
+
+   - Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+
+   - Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
+   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "kiss_fft.h"
+#include "celt.h"
+#include "modes.h"
+#include "arch.h"
+#include "quant_bands.h"
+#include <stdio.h>
+#ifndef FIXED_POINT
+#include "mlp.c"
+#include "mlp_data.c"
+#endif
+
+#ifndef M_PI
+#define M_PI 3.141592653
+#endif
+
+float dct_table[128] = {
+        0.250000, 0.250000, 0.250000, 0.250000, 0.250000, 0.250000, 0.250000, 0.250000,
+        0.250000, 0.250000, 0.250000, 0.250000, 0.250000, 0.250000, 0.250000, 0.250000,
+        0.351851, 0.338330, 0.311806, 0.273300, 0.224292, 0.166664, 0.102631, 0.034654,
+        -0.034654, -0.102631, -0.166664, -0.224292, -0.273300, -0.311806, -0.338330, -0.351851,
+        0.346760, 0.293969, 0.196424, 0.068975, -0.068975, -0.196424, -0.293969, -0.346760,
+        -0.346760, -0.293969, -0.196424, -0.068975, 0.068975, 0.196424, 0.293969, 0.346760,
+        0.338330, 0.224292, 0.034654, -0.166664, -0.311806, -0.351851, -0.273300, -0.102631,
+        0.102631, 0.273300, 0.351851, 0.311806, 0.166664, -0.034654, -0.224292, -0.338330,
+        0.326641, 0.135299, -0.135299, -0.326641, -0.326641, -0.135299, 0.135299, 0.326641,
+        0.326641, 0.135299, -0.135299, -0.326641, -0.326641, -0.135299, 0.135299, 0.326641,
+        0.311806, 0.034654, -0.273300, -0.338330, -0.102631, 0.224292, 0.351851, 0.166664,
+        -0.166664, -0.351851, -0.224292, 0.102631, 0.338330, 0.273300, -0.034654, -0.311806,
+        0.293969, -0.068975, -0.346760, -0.196424, 0.196424, 0.346760, 0.068975, -0.293969,
+        -0.293969, 0.068975, 0.346760, 0.196424, -0.196424, -0.346760, -0.068975, 0.293969,
+        0.273300, -0.166664, -0.338330, 0.034654, 0.351851, 0.102631, -0.311806, -0.224292,
+        0.224292, 0.311806, -0.102631, -0.351851, -0.034654, 0.338330, 0.166664, -0.273300,
+};
+
+#define NB_FRAMES 8
+
+#define NB_TBANDS 18
+static const int tbands[NB_TBANDS+1] = {
+       2,  4,  6,  8, 10, 12, 14, 16, 20, 24, 28, 32, 40, 48, 56, 68, 80, 96, 120
+};
+
+static const float tweight[NB_TBANDS+1] = {
+      .3, .4, .5, .6, .7, .8, .9, 1., 1., 1., 1., 1., 1., 1., .8, .7, .6, .5
+};
+
+#define NB_TONAL_SKIP_BANDS 9
+
+typedef struct {
+   float angle[240];
+   float d_angle[240];
+   float d2_angle[240];
+   float prev_band_tonality[NB_TBANDS];
+   float prev_tonality;
+   float E[NB_FRAMES][NB_TBANDS];
+   float lowE[NB_TBANDS], highE[NB_TBANDS];
+   float meanE[NB_TBANDS], meanRE[NB_TBANDS];
+   float mem[32];
+   float cmean[8];
+   float std[9];
+   float music_prob;
+   float Etracker;
+   float lowECount;
+   int E_count;
+   int last_music;
+   int last_transition;
+   int count;
+   int opus_bandwidth;
+} TonalityAnalysisState;
+
+void tonality_analysis(TonalityAnalysisState *tonal, AnalysisInfo *info, CELTEncoder *celt_enc, const opus_val16 *x, int C)
+{
+    int i, b;
+    const CELTMode *mode;
+    const kiss_fft_state *kfft;
+    kiss_fft_cpx in[480], out[480];
+    int N = 480, N2=240;
+    float * restrict A = tonal->angle;
+    float * restrict dA = tonal->d_angle;
+    float * restrict d2A = tonal->d2_angle;
+    float tonality[240];
+    float noisiness[240];
+    float band_tonality[NB_TBANDS];
+    float logE[NB_TBANDS];
+    float BFCC[8];
+    float features[100];
+    float frame_tonality;
+    float max_frame_tonality;
+    float tw_sum=0;
+    float frame_noisiness;
+    const float pi4 = M_PI*M_PI*M_PI*M_PI;
+    float slope=0;
+    float frame_stationarity;
+    float relativeE;
+    float frame_prob;
+    float alpha, alphaE, alphaE2;
+    float frame_loudness;
+    float bandwidth_mask;
+    int bandwidth=0;
+    float bandE[NB_TBANDS];
+    celt_encoder_ctl(celt_enc, CELT_GET_MODE(&mode));
+
+    tonal->last_transition++;
+    alpha = 1.f/IMIN(20, 1+tonal->count);
+    alphaE = 1.f/IMIN(50, 1+tonal->count);
+    alphaE2 = 1.f/IMIN(6000, 1+tonal->count);
+
+    if (tonal->count<4)
+       tonal->music_prob = .5;
+    kfft = mode->mdct.kfft[0];
+    if (C==1)
+    {
+       for (i=0;i<N2;i++)
+       {
+          float w = .5-.5*cos(M_PI*(i+1)/N2);
+          in[i].r = MULT16_16(w, x[i]);
+          in[i].i = MULT16_16(w, x[N-N2+i]);
+          in[N-i-1].r = MULT16_16(w, x[N-i-1]);
+          in[N-i-1].i = MULT16_16(w, x[2*N-N2-i-1]);
+       }
+    } else {
+       for (i=0;i<N2;i++)
+       {
+          float w = .5-.5*cos(M_PI*(i+1)/N2);
+          in[i].r = MULT16_16(w, x[2*i]+x[2*i+1]);
+          in[i].i = MULT16_16(w, x[2*(N-N2+i)]+x[2*(N-N2+i)+1]);
+          in[N-i-1].r = MULT16_16(w, x[2*(N-i-1)]+x[2*(N-i-1)+1]);
+          in[N-i-1].i = MULT16_16(w, x[2*(2*N-N2-i-1)]+x[2*(2*N-N2-i-1)+1]);
+       }
+    }
+    opus_fft(kfft, in, out);
+
+    for (i=1;i<N2;i++)
+    {
+       float X1r, X2r, X1i, X2i;
+       float angle, d_angle, d2_angle;
+       float angle2, d_angle2, d2_angle2;
+       float mod1, mod2, avg_mod;
+       X1r = out[i].r+out[N-i].r;
+       X1i = out[i].i-out[N-i].i;
+       X2r = out[i].i+out[N-i].i;
+       X2i = out[N-i].r-out[i].r;
+
+       angle = (.5/M_PI)*atan2(X1i, X1r);
+       d_angle = angle - A[i];
+       d2_angle = d_angle - dA[i];
+
+       angle2 = (.5/M_PI)*atan2(X2i, X2r);
+       d_angle2 = angle2 - angle;
+       d2_angle2 = d_angle2 - d_angle;
+
+       mod1 = d2_angle - floor(.5+d2_angle);
+       noisiness[i] = fabs(mod1);
+       mod1 *= mod1;
+       mod1 *= mod1;
+
+       mod2 = d2_angle2 - floor(.5+d2_angle2);
+       noisiness[i] += fabs(mod2);
+       mod2 *= mod2;
+       mod2 *= mod2;
+
+       avg_mod = .25*(d2A[i]+2*mod1+mod2);
+       tonality[i] = 1./(1+40*16*pi4*avg_mod)-.015;
+
+       A[i] = angle2;
+       dA[i] = d_angle2;
+       d2A[i] = mod2;
+    }
+
+    frame_tonality = 0;
+    max_frame_tonality = 0;
+    tw_sum = 0;
+    info->activity = 0;
+    frame_noisiness = 0;
+    frame_stationarity = 0;
+    if (!tonal->count)
+    {
+       for (b=0;b<NB_TBANDS;b++)
+       {
+          tonal->lowE[b] = 1e10;
+          tonal->highE[b] = -1e10;
+       }
+    }
+    relativeE = 0;
+    info->boost_amount[0]=info->boost_amount[1]=0;
+    info->boost_band[0]=info->boost_band[1]=0;
+    frame_loudness = 0;
+    bandwidth_mask = 0;
+    for (b=0;b<NB_TBANDS;b++)
+    {
+       float E=0, tE=0, nE=0;
+       float L1, L2;
+       float stationarity;
+       for (i=tbands[b];i<tbands[b+1];i++)
+       {
+          float binE = out[i].r*out[i].r + out[N-i].r*out[N-i].r
+                     + out[i].i*out[i].i + out[N-i].i*out[N-i].i;
+          E += binE;
+          tE += binE*tonality[i];
+          nE += binE*2*(.5-noisiness[i]);
+       }
+       bandE[b] = E;
+       tonal->E[tonal->E_count][b] = E;
+       frame_noisiness += nE/(1e-15+E);
+
+       frame_loudness += sqrt(E+1e-10);
+       /* Add a reasonable noise floor */
+       tonal->meanE[b] = (1-alphaE2)*tonal->meanE[b] + alphaE2*E;
+       tonal->meanRE[b] = (1-alphaE2)*tonal->meanRE[b] + alphaE2*sqrt(E);
+       /* 13 dB slope for spreading function */
+       bandwidth_mask = MAX32(.05*bandwidth_mask, E);
+       /* Checks if band looks like stationary noise or if it's below a (trivial) masking curve */
+       if (tonal->meanRE[b]*tonal->meanRE[b] < tonal->meanE[b]*.95 && E>.1*bandwidth_mask)
+          bandwidth = b;
+       logE[b] = log(E+1e-10);
+       tonal->lowE[b] = MIN32(logE[b], tonal->lowE[b]+.01);
+       tonal->highE[b] = MAX32(logE[b], tonal->highE[b]-.1);
+       if (tonal->highE[b] < tonal->lowE[b]+1)
+       {
+          tonal->highE[b]+=.5;
+          tonal->lowE[b]-=.5;
+       }
+       relativeE += (logE[b]-tonal->lowE[b])/(EPSILON+tonal->highE[b]-tonal->lowE[b]);
+
+       L1=L2=0;
+       for (i=0;i<NB_FRAMES;i++)
+       {
+          L1 += sqrt(tonal->E[i][b]);
+          L2 += tonal->E[i][b];
+       }
+
+       stationarity = MIN16(0.99,L1/sqrt(EPSILON+NB_FRAMES*L2));
+       stationarity *= stationarity;
+       stationarity *= stationarity;
+       frame_stationarity += stationarity;
+       /*band_tonality[b] = tE/(1e-15+E)*/;
+       band_tonality[b] = MAX16(tE/(EPSILON+E), stationarity*tonal->prev_band_tonality[b]);
+#if 0
+       if (b>=NB_TONAL_SKIP_BANDS)
+       {
+          frame_tonality += tweight[b]*band_tonality[b];
+          tw_sum += tweight[b];
+       }
+#else
+       frame_tonality += band_tonality[b];
+       if (b>=NB_TBANDS-NB_TONAL_SKIP_BANDS)
+          frame_tonality -= band_tonality[b-NB_TBANDS+NB_TONAL_SKIP_BANDS];
+#endif
+       max_frame_tonality = MAX16(max_frame_tonality, (1+.03*(b-NB_TBANDS))*frame_tonality);
+       slope += band_tonality[b]*(b-8);
+       /*printf("%f %f ", band_tonality[b], stationarity);*/
+       if (band_tonality[b] > info->boost_amount[1] && b>=7 && b < NB_TBANDS-1)
+       {
+          if (band_tonality[b] > info->boost_amount[0])
+          {
+             info->boost_amount[1] = info->boost_amount[0];
+             info->boost_band[1] = info->boost_band[0];
+             info->boost_amount[0] = band_tonality[b];
+             info->boost_band[0] = b;
+          } else {
+             info->boost_amount[1] = band_tonality[b];
+             info->boost_band[1] = b;
+          }
+       }
+       tonal->prev_band_tonality[b] = band_tonality[b];
+    }
+
+    frame_loudness = 20*log10(frame_loudness);
+    tonal->Etracker = MAX32(tonal->Etracker-.03, frame_loudness);
+    tonal->lowECount *= (1-alphaE);
+    if (frame_loudness < tonal->Etracker-30)
+       tonal->lowECount += alphaE;
+
+    for (i=0;i<8;i++)
+    {
+       float sum=0;
+       for (b=0;b<16;b++)
+          sum += dct_table[i*16+b]*logE[b];
+       BFCC[i] = sum;
+    }
+
+    frame_stationarity /= NB_TBANDS;
+    relativeE /= NB_TBANDS;
+    if (tonal->count<10)
+       relativeE = .5;
+    frame_noisiness /= NB_TBANDS;
+#if 1
+    info->activity = frame_noisiness + (1-frame_noisiness)*relativeE;
+#else
+    info->activity = .5*(1+frame_noisiness-frame_stationarity);
+#endif
+    frame_tonality = (max_frame_tonality/(NB_TBANDS-NB_TONAL_SKIP_BANDS));
+    frame_tonality = MAX16(frame_tonality, tonal->prev_tonality*.8);
+    tonal->prev_tonality = frame_tonality;
+    info->boost_amount[0] -= frame_tonality+.2;
+    info->boost_amount[1] -= frame_tonality+.2;
+    if (band_tonality[info->boost_band[0]] < band_tonality[info->boost_band[0]+1]+.15
+        || band_tonality[info->boost_band[0]] < band_tonality[info->boost_band[0]-1]+.15)
+       info->boost_amount[0]=0;
+    if (band_tonality[info->boost_band[1]] < band_tonality[info->boost_band[1]+1]+.15
+        || band_tonality[info->boost_band[1]] < band_tonality[info->boost_band[1]-1]+.15)
+       info->boost_amount[1]=0;
+
+    slope /= 8*8;
+    info->tonality_slope = slope;
+
+    tonal->E_count = (tonal->E_count+1)%NB_FRAMES;
+    tonal->count++;
+    info->tonality = frame_tonality;
+
+    for (i=0;i<4;i++)
+       features[i] = -0.12299*(BFCC[i]+tonal->mem[i+24]) + 0.49195*(tonal->mem[i]+tonal->mem[i+16]) + 0.69693*tonal->mem[i+8] - 1.4349*tonal->cmean[i];
+
+    for (i=0;i<4;i++)
+       tonal->cmean[i] = (1-alpha)*tonal->cmean[i] + alpha*BFCC[i];
+
+    for (i=0;i<4;i++)
+        features[4+i] = 0.63246*(BFCC[i]-tonal->mem[i+24]) + 0.31623*(tonal->mem[i]-tonal->mem[i+16]);
+    for (i=0;i<3;i++)
+        features[8+i] = 0.53452*(BFCC[i]+tonal->mem[i+24]) - 0.26726*(tonal->mem[i]+tonal->mem[i+16]) -0.53452*tonal->mem[i+8];
+
+    if (tonal->count > 5)
+    {
+       for (i=0;i<9;i++)
+          tonal->std[i] = (1-alpha)*tonal->std[i] + alpha*features[i]*features[i];
+    }
+
+    for (i=0;i<8;i++)
+    {
+       tonal->mem[i+24] = tonal->mem[i+16];
+       tonal->mem[i+16] = tonal->mem[i+8];
+       tonal->mem[i+8] = tonal->mem[i];
+       tonal->mem[i] = BFCC[i];
+    }
+    for (i=0;i<9;i++)
+       features[11+i] = sqrt(tonal->std[i]);
+    features[20] = info->tonality;
+    features[21] = info->activity;
+    features[22] = frame_stationarity;
+    features[23] = info->tonality_slope;
+    features[24] = tonal->lowECount;
+
+#ifndef FIXED_POINT
+    mlp_process(&net, features, &frame_prob);
+    /* Adds a "probability dead zone", with a cap on certainty */
+    frame_prob = .90*frame_prob*frame_prob*frame_prob;
+
+    frame_prob = .5*(frame_prob+1);
+
+    /*printf("%f\n", frame_prob);*/
+    {
+       float tau, beta;
+       float p0, p1;
+       float max_certainty;
+       /* One transition every 3 minutes */
+       tau = .00005;
+       beta = .1;
+       max_certainty = 1.f/(10+1*tonal->last_transition);
+       p0 = (1-tonal->music_prob)*(1-tau) +    tonal->music_prob *tau;
+       p1 =    tonal->music_prob *(1-tau) + (1-tonal->music_prob)*tau;
+       p0 *= pow(1-frame_prob, beta);
+       p1 *= pow(frame_prob, beta);
+       tonal->music_prob = MAX16(max_certainty, MIN16(1-max_certainty, p1/(p0+p1)));
+       info->music_prob = tonal->music_prob;
+       /*printf("%f %f\n", frame_prob, info->music_prob);*/
+    }
+    if (tonal->last_music != (tonal->music_prob>.5))
+       tonal->last_transition=0;
+    tonal->last_music = tonal->music_prob>.5;
+#else
+    info->music_prob = 0;
+#endif
+    /*for (i=0;i<25;i++)
+       printf("%f ", features[i]);
+    printf("\n");*/
+
+    /* FIXME: Can't detect SWB for now because the last band ends at 12 kHz */
+    if (bandwidth == NB_TBANDS-1 || tonal->count<100)
+    {
+       tonal->opus_bandwidth = OPUS_BANDWIDTH_FULLBAND;
+    } else {
+       int close_enough = 0;
+       if (bandE[bandwidth-1] < 3000*bandE[NB_TBANDS-1] && bandwidth < NB_TBANDS-1)
+          close_enough=1;
+       if (bandwidth<=11 || (bandwidth==12 && close_enough))
+          tonal->opus_bandwidth = OPUS_BANDWIDTH_NARROWBAND;
+       else if (bandwidth<=13)
+          tonal->opus_bandwidth = OPUS_BANDWIDTH_MEDIUMBAND;
+       else if (bandwidth<=15 || (bandwidth==16 && close_enough))
+          tonal->opus_bandwidth = OPUS_BANDWIDTH_WIDEBAND;
+    }
+    info->noisiness = frame_noisiness;
+    info->valid = 1;
+}

src/mlp.c

+/* Copyright (c) 2008-2011 Octasic Inc.
+   Written by Jean-Marc Valin */
+/*
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions
+   are met:
+
+   - Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+
+   - Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
+   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+
+#include <math.h>
+#include "mlp.h"
+#include "arch.h"
+#include "tansig_table.h"
+#define MAX_NEURONS 100
+
+#ifdef FIXED_POINT
+extern const opus_val16 tansig_table[501];
+static inline opus_val16 tansig_approx(opus_val32 _x) /* Q19 */
+{
+	int i;
+	opus_val16 xx; /* Q11 */
+	/*double x, y;*/
+	opus_val16 dy, yy; /* Q14 */
+	/*x = 1.9073e-06*_x;*/
+	if (_x>=QCONST32(10,19))
+		return QCONST32(1.,14);
+	if (_x<=-QCONST32(10,19))
+		return -QCONST32(1.,14);
+	xx = EXTRACT16(SHR32(_x, 8));
+	/*i = lrint(25*x);*/
+	i = SHR32(ADD32(1024,MULT16_16(25, xx)),11);
+	/*x -= .04*i;*/
+	xx -= EXTRACT16(SHR32(MULT16_16(20972,i),8));
+	/*x = xx*(1./2048);*/
+	/*y = tansig_table[250+i];*/
+	yy = tansig_table[250+i];
+	/*y = yy*(1./16384);*/
+	dy = 16384-MULT16_16_Q14(yy,yy);
+	yy = yy + MULT16_16_Q14(MULT16_16_Q11(xx,dy),(16384 - MULT16_16_Q11(yy,xx)));
+	return yy;
+}
+#else
+/*extern const float tansig_table[501];*/
+static inline double tansig_approx(double x)
+{
+	int i;
+	double y, dy;
+	if (x>=10)
+		return 1;
+	if (x<=-10)
+		return -1;
+	i = lrint(25*x);
+	x -= .04*i;
+	y = tansig_table[250+i];
+	dy = 1-y*y;
+	y = y + x*dy*(1 - y*x);
+	return y;
+}
+#endif
+
+void mlp_process(const MLP *m, const opus_val16 *in, opus_val16 *out)
+{
+	int j;
+	opus_val16 hidden[MAX_NEURONS];
+	const opus_val16 *W = m->weights;
+	/* Copy to tmp_in */
+	for (j=0;j<m->topo[1];j++)
+	{
+		int k;
+		opus_val32 sum = SHL32(EXTEND32(*W++),8);
+		for (k=0;k<m->topo[0];k++)
+			sum = MAC16_16(sum, in[k],*W++);
+		hidden[j] = tansig_approx(sum);
+	}
+	for (j=0;j<m->topo[2];j++)
+	{
+		int k;
+		opus_val32 sum = SHL32(EXTEND32(*W++),14);
+		for (k=0;k<m->topo[1];k++)
+			sum = MAC16_16(sum, hidden[k], *W++);
+		out[j] = tansig_approx(EXTRACT16(PSHR32(sum,17)));
+	}
+}
+

src/mlp.h

+/* Copyright (c) 2008-2011 Octasic Inc.
+   Written by Jean-Marc Valin */
+/*
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions
+   are met:
+
+   - Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+
+   - Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
+   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+#ifndef _MLP_H_
+#define _MLP_H_
+
+#include "arch.h"
+
+typedef struct {
+	int layers;
+	const int *topo;
+	const opus_val16 *weights;
+} MLP;
+
+void mlp_process(const MLP *m, const opus_val16 *in, opus_val16 *out);
+
+#endif /* _MLP_H_ */

src/mlp_data.c

+#include "mlp.h"
+
+/* RMS error was 0.179835, seed was 1322103961 */
+
+static const float weights[271] = {
+
+/* hidden layer */
+1.55597, -0.0739792, -0.0646761, -0.099531, -0.0794943, 
+0.0180174, -0.0391354, 0.0508224, -0.0160169, -0.0773263, 
+-0.0300002, -0.0865361, 0.124477, -0.28648, -0.0860702, 
+-0.518949, -0.0873341, -0.235393, -0.907833, -0.383573, 
+0.535388, -0.57944, 0.98116, 0.8482, 1.12426, 
+-3.23721, -0.647072, -0.0265139, 0.0711052, -0.00125666, 
+-0.0396181, -0.44282, -0.510495, -0.201865, 0.0134336, 
+-0.167205, -0.155406, 0.00041678, -0.00468705, -0.0233224, 
+0.264279, -0.301375, 0.00234895, 0.0144741, -0.137535, 
+0.200323, 0.0192027, 3.19818, 2.03495, 0.705517, 
+-4.6025, -0.11485, -0.792716, 0.150714, 0.10608, 
+0.240633, 0.0690698, 0.0695297, 0.124819, 0.0501433, 
+0.0460952, 0.147639, 0.10327, 0.158007, 0.113714, 
+0.0276191, 0.0680749, -0.130012, 0.0796126, 0.133067, 
+0.51495, 0.747578, -0.128742, 5.98112, -1.16698, 
+-0.276492, -1.73549, -3.90234, 2.01489, -0.040118, 
+-0.113002, -0.146751, -0.113569, 0.0534873, 0.0989832, 
+0.0872875, 0.049266, 0.0367557, -0.00889148, -0.0648461, 
+-0.00190352, 0.0143773, 0.0259364, -0.0592133, -0.0672924, 
+0.1399, -0.0987886, -0.347402, 0.101326, -0.0680876, 
+0.469186, 0.246922, 10.4017, 3.44846, -0.662725, 
+-0.0328208, -0.0561274, -0.0167744, 0.00044282, -0.0457645, 
+-0.0408314, -0.013113, -0.0373873, -0.0474122, -0.0273745, 
+-0.0308505, 0.000582959, -0.0421135, 0.464859, 0.196842, 
+0.320538, 0.0435528, -0.200168, 0.266475, -0.0853727, 
+1.20397, 0.711542, -1.04397, -1.47759, 1.26768, 
+0.446958, 0.266477, -0.30802, 0.28431, -0.118541, 
+0.00836345, 0.0689026, -0.0137996, -0.0395417, 0.26982, 
+-0.206255, 0.16066, 0.114757, 0.359587, -0.106503, 
+-0.0948534, 0.175358, -0.122966, -0.0056675, 0.483848, 
+-0.134916, -0.427567, -0.140172, -1.0866, -2.73921, 
+0.549843, 0.17685, 0.0010675, -0.00137386, 0.0884424, 
+-0.0698736, -0.00174136, 0.0718775, -0.0396849, 0.0448056, 
+0.0577853, -0.0372353, 0.134599, 0.0260656, 0.140322, 
+0.22704, -0.020568, -0.0142424, -0.21723, -0.997704, 
+-0.884573, -0.163495, 2.33617, 0.224142, 0.19635, 
+-0.957387, 0.144678, 1.47035, -0.00700498, -0.0472309, 
+-0.0137848, -0.0189145, 0.00856479, 0.0316965, 0.00613373, 
+0.00209807, 0.00270964, -0.0490206, 0.0105712, -0.0465045, 
+-0.0381532, -0.0985268, -0.108297, 0.0146409, -0.0040718, 
+-0.0698572, -0.380568, -0.230479, 3.98917, 0.457652, 
+-1.02355, -7.4435, -0.475314, 1.61743, 0.0254017, 
+-0.00791293, 0.047217, 0.0220995, -0.0304311, 0.0052168, 
+-0.0404054, -0.0230293, 0.00169229, -0.0138178, 0.0043137, 
+-0.0598088, -0.133601, 0.0555138, -0.177358, -0.159856, 
+-0.137281, 0.108051, -0.305973, 0.393775, 0.0747287, 
+0.783993, -0.875086, 1.06862, 0.340519, -0.352681, 
+-0.0830912, -0.100017, 0.0729085, -0.00829403, 0.027489, 
+-0.0779597, 0.082286, -0.164181, -0.41519, 0.00282335, 
+-0.29573, 0.125571, 0.726935, 0.392137, 0.491348, 
+0.0723196, -0.0259758, -0.0636332, -0.452384, -0.000225974, 
+-2.34001, 2.45211, -0.544628, 5.62944, -3.44507, 
+
+/* output layer */
+-3.13835, 0.994751, 0.444901, 1.59518, 1.23665, 
+3.37012, -1.34606, 1.99131, 1.33476, 1.3885, 
+1.12559, };
+
+static const int topo[3] = {25, 10, 1};
+
+const MLP net = {
+	3,
+	topo,
+	weights
+};
+

src/mlp_train.c

+/* Copyright (c) 2008-2011 Octasic Inc.
+   Written by Jean-Marc Valin */
+/*
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions
+   are met:
+
+   - Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+
+   - Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
+   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+
+#include "mlp_train.h"
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <semaphore.h>
+#include <pthread.h>
+#include <time.h>
+#include <signal.h>
+
+int stopped = 0;
+
+void handler(int sig)
+{
+	stopped = 1;
+	signal(sig, handler);
+}
+
+MLPTrain * mlp_init(int *topo, int nbLayers, float *inputs, float *outputs, int nbSamples)
+{
+	int i, j, k;
+	MLPTrain *net;
+	int inDim, outDim;
+	net = malloc(sizeof(*net));
+	net->topo = malloc(nbLayers*sizeof(net->topo[0]));
+	for (i=0;i<nbLayers;i++)
+		net->topo[i] = topo[i];
+	inDim = topo[0];
+	outDim = topo[nbLayers-1];
+	net->in_rate = malloc((inDim+1)*sizeof(net->in_rate[0]));
+	net->weights = malloc((nbLayers-1)*sizeof(net->weights));
+	net->best_weights = malloc((nbLayers-1)*sizeof(net->weights));
+	for (i=0;i<nbLayers-1;i++)
+	{
+		net->weights[i] = malloc((topo[i]+1)*topo[i+1]*sizeof(net->weights[0][0]));
+		net->best_weights[i] = malloc((topo[i]+1)*topo[i+1]*sizeof(net->weights[0][0]));
+	}
+	double inMean[inDim];
+	for (j=0;j<inDim;j++)
+	{
+		double std=0;
+		inMean[j] = 0;
+		for (i=0;i<nbSamples;i++)
+		{
+			inMean[j] += inputs[i*inDim+j];
+			std += inputs[i*inDim+j]*inputs[i*inDim+j];
+		}
+		inMean[j] /= nbSamples;
+		std /= nbSamples;
+		net->in_rate[1+j] = .5/(.0001+std);
+		std = std-inMean[j]*inMean[j];
+		if (std<.001)
+			std = .001;
+		std = 1/sqrt(inDim*std);
+		for (k=0;k<topo[1];k++)
+			net->weights[0][k*(topo[0]+1)+j+1] = randn(std);
+	}
+	net->in_rate[0] = 1;
+	for (j=0;j<topo[1];j++)
+	{
+		double sum = 0;
+		for (k=0;k<inDim;k++)
+			sum += inMean[k]*net->weights[0][j*(topo[0]+1)+k+1];
+		net->weights[0][j*(topo[0]+1)] = -sum;
+	}
+	for (j=0;j<outDim;j++)
+	{
+		double mean = 0;
+		double std;
+		for (i=0;i<nbSamples;i++)
+			mean += outputs[i*outDim+j];
+		mean /= nbSamples;
+		std = 1/sqrt(topo[nbLayers-2]);
+		net->weights[nbLayers-2][j*(topo[nbLayers-2]+1)] = mean;
+		for (k=0;k<topo[nbLayers-2];k++)
+			net->weights[nbLayers-2][j*(topo[nbLayers-2]+1)+k+1] = randn(std);
+	}
+	return net;
+}
+
+#define MAX_NEURONS 100
+
+double compute_gradient(MLPTrain *net, float *inputs, float *outputs, int nbSamples, double *W0_grad, double *W1_grad, double *error_rate)
+{
+	int i,j;
+	int s;
+	int inDim, outDim, hiddenDim;
+	int *topo;
+	double *W0, *W1;
+	double rms=0;
+	int W0_size, W1_size;
+	double hidden[MAX_NEURONS];
+	double netOut[MAX_NEURONS];
+	double error[MAX_NEURONS];
+
+        *error_rate = 0;
+	topo = net->topo;
+	inDim = net->topo[0];
+	hiddenDim = net->topo[1];
+	outDim = net->topo[2];
+	W0_size = (topo[0]+1)*topo[1];
+	W1_size = (topo[1]+1)*topo[2];
+	W0 = net->weights[0];
+	W1 = net->weights[1];
+	memset(W0_grad, 0, W0_size*sizeof(double));
+	memset(W1_grad, 0, W1_size*sizeof(double));
+	for (i=0;i<outDim;i++)
+		netOut[i] = outputs[i];
+	for (s=0;s<nbSamples;s++)
+	{
+		float *in, *out;
+		in = inputs+s*inDim;
+		out = outputs + s*outDim;
+		for (i=0;i<hiddenDim;i++)
+		{
+			double sum = W0[i*(inDim+1)];
+			for (j=0;j<inDim;j++)
+				sum += W0[i*(inDim+1)+j+1]*in[j];
+			hidden[i] = tansig_approx(sum);
+		}
+		for (i=0;i<outDim;i++)
+		{
+			double sum = W1[i*(hiddenDim+1)];
+			for (j=0;j<hiddenDim;j++)
+				sum += W1[i*(hiddenDim+1)+j+1]*hidden[j];
+			netOut[i] = tansig_approx(sum);
+			error[i] = out[i] - netOut[i];
+			rms += error[i]*error[i];
+			*error_rate += fabs(error[i])>1;
+                        //error[i] = error[i]/(1+fabs(error[i]));
+		}
+		/* Back-propagate error */
+		for (i=0;i<outDim;i++)
+		{
+                        float grad = 1-netOut[i]*netOut[i];
+			W1_grad[i*(hiddenDim+1)] += error[i]*grad;
+			for (j=0;j<hiddenDim;j++)
+				W1_grad[i*(hiddenDim+1)+j+1] += grad*error[i]*hidden[j];
+		}
+		for (i=0;i<hiddenDim;i++)
+		{
+			double grad;
+			grad = 0;
+			for (j=0;j<outDim;j++)
+				grad += error[j]*W1[j*(hiddenDim+1)+i+1];
+			grad *= 1-hidden[i]*hidden[i];
+			W0_grad[i*(inDim+1)] += grad;
+			for (j=0;j<inDim;j++)
+				W0_grad[i*(inDim+1)+j+1] += grad*in[j];
+		}
+	}
+	return rms;
+}
+
+#define NB_THREADS 8
+
+sem_t sem_begin[NB_THREADS];
+sem_t sem_end[NB_THREADS];
+
+struct GradientArg {
+	int id;
+	int done;
+	MLPTrain *net;
+	float *inputs;
+	float *outputs;
+	int nbSamples;
+	double *W0_grad;
+	double *W1_grad;
+	double rms;
+	double error_rate;
+};
+
+void *gradient_thread_process(void *_arg)
+{
+	int W0_size, W1_size;
+	struct GradientArg *arg = _arg;
+	int *topo = arg->net->topo;
+	W0_size = (topo[0]+1)*topo[1];
+	W1_size = (topo[1]+1)*topo[2];
+	double W0_grad[W0_size];
+	double W1_grad[W1_size];
+	arg->W0_grad = W0_grad;
+	arg->W1_grad = W1_grad;
+	while (1)
+	{
+		sem_wait(&sem_begin[arg->id]);
+		if (arg->done)
+			break;
+		arg->rms = compute_gradient(arg->net, arg->inputs, arg->outputs, arg->nbSamples, arg->W0_grad, arg->W1_grad, &arg->error_rate);
+		sem_post(&sem_end[arg->id]);
+	}
+	fprintf(stderr, "done\n");
+	return NULL;
+}
+
+float mlp_train_backprop(MLPTrain *net, float *inputs, float *outputs, int nbSamples, int nbEpoch, float rate)
+{
+	int i, j;
+	int e;
+	float best_rms = 1e10;
+	int inDim, outDim, hiddenDim;
+	int *topo;
+	double *W0, *W1, *best_W0, *best_W1;
+	double *W0_old, *W1_old;
+	double *W0_old2, *W1_old2;
+	double *W0_grad, *W1_grad;
+	double *W0_oldgrad, *W1_oldgrad;
+	double *W0_rate, *W1_rate;
+	double *best_W0_rate, *best_W1_rate;
+	int W0_size, W1_size;
+	topo = net->topo;
+	W0_size = (topo[0]+1)*topo[1];
+	W1_size = (topo[1]+1)*topo[2];
+	struct GradientArg args[NB_THREADS];
+	pthread_t thread[NB_THREADS];
+	int samplePerPart = nbSamples/NB_THREADS;
+	int count_worse=0;
+	int count_retries=0;
+
+	topo = net->topo;
+	inDim = net->topo[0];
+	hiddenDim = net->topo[1];
+	outDim = net->topo[2];
+	W0 = net->weights[0];
+	W1 = net->weights[1];
+	best_W0 = net->best_weights[0];
+	best_W1 = net->best_weights[1];
+	W0_old = malloc(W0_size*sizeof(double));
+	W1_old = malloc(W1_size*sizeof(double));
+	W0_old2 = malloc(W0_size*sizeof(double));
+	W1_old2 = malloc(W1_size*sizeof(double));
+	W0_grad = malloc(W0_size*sizeof(double));
+	W1_grad = malloc(W1_size*sizeof(double));
+	W0_oldgrad = malloc(W0_size*sizeof(double));
+	W1_oldgrad = malloc(W1_size*sizeof(double));
+	W0_rate = malloc(W0_size*sizeof(double));
+	W1_rate = malloc(W1_size*sizeof(double));
+	best_W0_rate = malloc(W0_size*sizeof(double));
+	best_W1_rate = malloc(W1_size*sizeof(double));
+	memcpy(W0_old, W0, W0_size*sizeof(double));
+	memcpy(W0_old2, W0, W0_size*sizeof(double));
+	memset(W0_grad, 0, W0_size*sizeof(double));
+	memset(W0_oldgrad, 0, W0_size*sizeof(double));
+	memcpy(W1_old, W1, W1_size*sizeof(double));
+	memcpy(W1_old2, W1, W1_size*sizeof(double));
+	memset(W1_grad, 0, W1_size*sizeof(double));
+	memset(W1_oldgrad, 0, W1_size*sizeof(double));
+	
+	rate /= nbSamples;
+	for (i=0;i<hiddenDim;i++)
+		for (j=0;j<inDim+1;j++)
+			W0_rate[i*(inDim+1)+j] = rate*net->in_rate[j];
+	for (i=0;i<W1_size;i++)
+		W1_rate[i] = rate;
+	
+	for (i=0;i<NB_THREADS;i++)
+	{
+		args[i].net = net;
+		args[i].inputs = inputs+i*samplePerPart*inDim;
+		args[i].outputs = outputs+i*samplePerPart*outDim;
+		args[i].nbSamples = samplePerPart;
+		args[i].id = i;
+		args[i].done = 0;
+		sem_init(&sem_begin[i], 0, 0);
+		sem_init(&sem_end[i], 0, 0);
+		pthread_create(&thread[i], NULL, gradient_thread_process, &args[i]);
+	}
+	for (e=0;e<nbEpoch;e++)
+	{
+		double rms=0;
+                double error_rate = 0;
+		for (i=0;i<NB_THREADS;i++)
+		{
+			sem_post(&sem_begin[i]);
+		}
+		memset(W0_grad, 0, W0_size*sizeof(double));
+		memset(W1_grad, 0, W1_size*sizeof(double));
+		for (i=0;i<NB_THREADS;i++)
+		{
+			sem_wait(&sem_end[i]);
+			rms += args[i].rms;
+			error_rate += args[i].error_rate;
+			for (j=0;j<W0_size;j++)
+				W0_grad[j] += args[i].W0_grad[j];
+			for (j=0;j<W1_size;j++)
+				W1_grad[j] += args[i].W1_grad[j];
+		}
+
+		float mean_rate = 0, min_rate = 1e10;
+		rms = (rms/(outDim*nbSamples));
+		error_rate = (error_rate/(outDim*nbSamples));
+		fprintf (stderr, "%f (%f %f) ", error_rate, rms, best_rms);
+		if (rms < best_rms)
+		{
+			best_rms = rms;
+			for (i=0;i<W0_size;i++)
+			{
+				best_W0[i] = W0[i];
+				best_W0_rate[i] = W0_rate[i];
+			}
+			for (i=0;i<W1_size;i++)
+			{
+				best_W1[i] = W1[i];
+				best_W1_rate[i] = W1_rate[i];
+			}
+			count_worse=0;
+			count_retries=0;
+		} else {
+			count_worse++;
+			if (count_worse>30)
+			{
+			    count_retries++;
+				count_worse=0;
+				for (i=0;i<W0_size;i++)
+				{
+					W0[i] = best_W0[i];
+					best_W0_rate[i] *= .7;
+					if (best_W0_rate[i]<1e-15) best_W0_rate[i]=1e-15;
+					W0_rate[i] = best_W0_rate[i];
+					W0_grad[i] = 0;
+				}
+				for (i=0;i<W1_size;i++)
+				{
+					W1[i] = best_W1[i];
+					best_W1_rate[i] *= .8;
+					if (best_W1_rate[i]<1e-15) best_W1_rate[i]=1e-15;
+					W1_rate[i] = best_W1_rate[i];
+					W1_grad[i] = 0;
+				}
+			}
+		}
+		if (count_retries>10)
+		    break;
+		for (i=0;i<W0_size;i++)
+		{
+			if (W0_oldgrad[i]*W0_grad[i] > 0)
+				W0_rate[i] *= 1.01;
+			else if (W0_oldgrad[i]*W0_grad[i] < 0)
+				W0_rate[i] *= .9;
+			mean_rate += W0_rate[i];
+			if (W0_rate[i] < min_rate)
+				min_rate = W0_rate[i];
+			if (W0_rate[i] < 1e-15)
+				W0_rate[i] = 1e-15;
+			/*if (W0_rate[i] > .01)
+				W0_rate[i] = .01;*/
+			W0_oldgrad[i] = W0_grad[i];
+			W0_old2[i] = W0_old[i];
+			W0_old[i] = W0[i];
+			W0[i] += W0_grad[i]*W0_rate[i];
+		}
+		for (i=0;i<W1_size;i++)
+		{
+			if (W1_oldgrad[i]*W1_grad[i] > 0)
+				W1_rate[i] *= 1.01;
+			else if (W1_oldgrad[i]*W1_grad[i] < 0)
+				W1_rate[i] *= .9;
+			mean_rate += W1_rate[i];
+			if (W1_rate[i] < min_rate)
+				min_rate = W1_rate[i];
+			if (W1_rate[i] < 1e-15)
+				W1_rate[i] = 1e-15;
+			W1_oldgrad[i] = W1_grad[i];
+			W1_old2[i] = W1_old[i];
+			W1_old[i] = W1[i];
+			W1[i] += W1_grad[i]*W1_rate[i];
+		}
+		mean_rate /= (topo[0]+1)*topo[1] + (topo[1]+1)*topo[2];
+		fprintf (stderr, "%g %d", mean_rate, e);
+		if (count_retries)
+		    fprintf(stderr, " %d", count_retries);
+		fprintf(stderr, "\n");
+		if (stopped)
+			break;
+	}
+	for (i=0;i<NB_THREADS;i++)
+	{
+		args[i].done = 1;
+		sem_post(&sem_begin[i]);
+		pthread_join(thread[i], NULL);
+		fprintf (stderr, "joined %d\n", i);
+	}
+	free(W0_old);
+	free(W1_old);
+	free(W0_grad);
+	free(W1_grad);
+	free(W0_rate);
+	free(W1_rate);
+	return best_rms;
+}
+
+int main(int argc, char **argv)
+{
+	int i, j;
+	int nbInputs;
+	int nbOutputs;
+	int nbHidden;
+	int nbSamples;
+	int nbEpoch;
+	int nbRealInputs;
+	unsigned int seed;
+	int ret;
+	float rms;
+	float *inputs;
+	float *outputs;
+	if (argc!=6)
+	{
+		fprintf (stderr, "usage: mlp_train <inputs> <hidden> <outputs> <nb samples> <nb epoch>\n");
+		return 1;
+	}
+	nbInputs = atoi(argv[1]);
+	nbHidden = atoi(argv[2]);
+	nbOutputs = atoi(argv[3]);
+	nbSamples = atoi(argv[4]);
+	nbEpoch = atoi(argv[5]);
+	nbRealInputs = nbInputs;
+	inputs = malloc(nbInputs*nbSamples*sizeof(*inputs));
+	outputs = malloc(nbOutputs*nbSamples*sizeof(*outputs));
+	
+	seed = time(NULL);
+	fprintf (stderr, "Seed is %u\n", seed);
+	srand(seed);
+	build_tansig_table();
+	signal(SIGTERM, handler);
+	signal(SIGINT, handler);
+	signal(SIGHUP, handler);
+	for (i=0;i<nbSamples;i++)
+	{
+		for (j=0;j<nbRealInputs;j++)
+			ret = scanf(" %f", &inputs[i*nbInputs+j]);
+		for (j=0;j<nbOutputs;j++)
+			ret = scanf(" %f", &outputs[i*nbOutputs+j]);
+		if (feof(stdin))
+		{
+			nbSamples = i;
+			break;
+		}
+	}
+	int topo[3] = {nbInputs, nbHidden, nbOutputs};
+	MLPTrain *net;
+
+	fprintf (stderr, "Got %d samples\n", nbSamples);
+	net = mlp_init(topo, 3, inputs, outputs, nbSamples);
+	rms = mlp_train_backprop(net, inputs, outputs, nbSamples, nbEpoch, 1);
+	printf ("#include \"mlp.h\"\n\n");
+	printf ("/* RMS error was %f, seed was %u */\n\n", rms, seed);
+	printf ("static const float weights[%d] = {\n", (topo[0]+1)*topo[1] + (topo[1]+1)*topo[2]);
+	printf ("\n/* hidden layer */\n");
+	for (i=0;i<(topo[0]+1)*topo[1];i++)
+	{
+		printf ("%g, ", net->weights[0][i]);
+		if (i%5==4)
+			printf("\n");
+	}
+	printf ("\n/* output layer */\n");
+	for (i=0;i<(topo[1]+1)*topo[2];i++)
+	{
+		printf ("%g, ", net->weights[1][i]);
+		if (i%5==4)
+			printf("\n");
+	}
+	printf ("};\n\n");
+	printf ("static const int topo[3] = {%d, %d, %d};\n\n", topo[0], topo[1], topo[2]);
+	printf ("const MLP net = {\n");
+	printf ("\t3,\n");
+	printf ("\ttopo,\n");
+	printf ("\tweights\n};\n");
+	return 0;
+}

src/mlp_train.h

+/* Copyright (c) 2008-2011 Octasic Inc.
+   Written by Jean-Marc Valin */
+/*
+   Redistribution and use in source and binary forms, with or without
+   modification, are permitted provided that the following conditions
+   are met:
+
+   - Redistributions of source code must retain the above copyright
+   notice, this list of conditions and the following disclaimer.
+
+   - Redistributions in binary form must reproduce the above copyright
+   notice, this list of conditions and the following disclaimer in the
+   documentation and/or other materials provided with the distribution.
+
+   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
+   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+#ifndef _MLP_TRAIN_H_
+#define _MLP_TRAIN_H_
+
+#include <math.h>
+#include <stdlib.h>
+
+double tansig_table[501];
+static inline double tansig_double(double x) 
+{
+	return 2./(1.+exp(-2.*x)) - 1.;
+}
+static inline void build_tansig_table()
+{
+	int i;
+	for (i=0;i<501;i++)
+		tansig_table[i] = tansig_double(.04*(i-250));
+}
+
+static inline double tansig_approx(double x)
+{
+	int i;
+	double y, dy;
+	if (x>=10)
+		return 1;
+	if (x<=-10)
+		return -1;
+	i = lrint(25*x);
+	x -= .04*i;
+	y = tansig_table[250+i];
+	dy = 1-y*y;
+	y = y + x*dy*(1 - y*x);
+	return y;
+}
+
+inline float randn(float sd)
+{
+   float U1, U2, S, x;
+   do {
+      U1 = ((float)rand())/RAND_MAX;
+      U2 = ((float)rand())/RAND_MAX;
+      U1 = 2*U1-1;
+      U2 = 2*U2-1;
+      S = U1*U1 + U2*U2;
+   } while (S >= 1 || S == 0.0f);
+   x = sd*sqrt(-2 * log(S) / S) * U1;
+   return x;
+}
+
+
+typedef struct {
+	int layers;
+	int *topo;
+	double **weights;
+	double **best_weights;
+	double *in_rate;
+} MLPTrain;
+
+
+#endif /* _MLP_TRAIN_H_ */

src/opus_encoder.c

@@ -40,6 +40,7 @@
 #include "arch.h"
 #include "opus_private.h"
 #include "os_support.h"
+#include "analysis.c"
 
 #include "tuning_parameters.h"
 #ifdef FIXED_POINT
@@ -84,7 +85,9 @@ struct OpusEncoder {
     /* Sampling rate (at the API level) */
     int          first;
     opus_val16   delay_buffer[MAX_ENCODER_BUFFER*2];
-
+#ifndef FIXED_POINT
+    TonalityAnalysisState analysis;
+#endif
     opus_uint32  rangeFinal;
 };
 
@@ -365,6 +368,56 @@ static void hp_cutoff(const opus_val16 *in, opus_int32 cutoff_Hz, opus_val16 *ou
 #endif
 }
 
+#ifdef FIXED_POINT
+static void dc_reject(const opus_val16 *in, opus_int32 cutoff_Hz, opus_val16 *out, opus_val32 *hp_mem, int len, int channels, opus_int32 Fs)
+{
+   int c, i;
+   int shift;
+
+   /* Approximates -round(log2(4.*cutoff_Hz/Fs)) */
+   shift=celt_ilog2(Fs/(cutoff_Hz*3));
+   for (c=0;c<channels;c++)
+   {
+      for (i=0;i<len;i++)
+      {
+         opus_val32 x, tmp, y;
+         x = SHL32(EXTEND32(in[channels*i+c]), 15);
+         /* First stage */
+         tmp = x-hp_mem[2*c];
+         hp_mem[2*c] = hp_mem[2*c] + PSHR32(x - hp_mem[2*c], shift);
+         /* Second stage */
+         y = tmp - hp_mem[2*c+1];
+         hp_mem[2*c+1] = hp_mem[2*c+1] + PSHR32(tmp - hp_mem[2*c+1], shift);
+         out[channels*i+c] = EXTRACT16(SATURATE(PSHR32(y, 15), 32767));
+      }
+   }
+}
+
+#else
+static void dc_reject(const opus_val16 *in, opus_int32 cutoff_Hz, opus_val16 *out, opus_val32 *hp_mem, int len, int channels, opus_int32 Fs)
+{
+   int c, i;
+   float coef;
+
+   coef = 4.*cutoff_Hz/Fs;
+   for (c=0;c<channels;c++)
+   {
+      for (i=0;i<len;i++)
+      {
+         opus_val32 x, tmp, y;
+         x = in[channels*i+c];
+         /* First stage */
+         tmp = x-hp_mem[2*c];
+         hp_mem[2*c] = hp_mem[2*c] + coef*(x - hp_mem[2*c]);
+         /* Second stage */
+         y = tmp - hp_mem[2*c+1];
+         hp_mem[2*c+1] = hp_mem[2*c+1] + coef*(tmp - hp_mem[2*c+1]);
+         out[channels*i+c] = y;
+      }
+   }
+}
+#endif
+
 static void stereo_fade(const opus_val16 *in, opus_val16 *out, opus_val16 g1, opus_val16 g2,
         int overlap, int frame_size, int channels, const opus_val16 *window, opus_int32 Fs)
 {
@@ -469,6 +522,11 @@ opus_int32 opus_encode_float(OpusEncoder *st, const opus_val16 *pcm, int frame_s
     opus_int32 max_rate;
     int curr_bandwidth;
     opus_int32 max_data_bytes;
+    int extra_buffer, total_buffer;
+    int perform_analysis=0;
+#ifndef FIXED_POINT
+    AnalysisInfo analysis_info;
+#endif
     VARDECL(opus_val16, tmp_prefill);
 
     ALLOC_STACK;
@@ -490,11 +548,20 @@ opus_int32 opus_encode_float(OpusEncoder *st, const opus_val16 *pcm, int frame_s
     silk_enc = (char*)st+st->silk_enc_offset;
     celt_enc = (CELTEncoder*)((char*)st+st->celt_enc_offset);
 
+#ifndef FIXED_POINT
+    perform_analysis = st->silk_mode.complexity >= 7 && frame_size >= st->Fs/100 && st->Fs==48000;
+#endif
     if (st->application == OPUS_APPLICATION_RESTRICTED_LOWDELAY)
        delay_compensation = 0;
     else
        delay_compensation = st->delay_compensation;
-
+    if (perform_analysis)
+    {
+       total_buffer = IMAX(st->Fs/200, delay_compensation);
+    } else {
+       total_buffer = delay_compensation;
+    }
+    extra_buffer = total_buffer-delay_compensation;
     st->bitrate_bps = user_bitrate_to_bitrate(st, frame_size, max_data_bytes);
 
     frame_rate = st->Fs/frame_size;
@@ -820,9 +887,9 @@ opus_int32 opus_encode_float(OpusEncoder *st, const opus_val16 *pcm, int frame_s
 
     ec_enc_init(&enc, data, max_data_bytes-1);
 
-    ALLOC(pcm_buf, (delay_compensation+frame_size)*st->channels, opus_val16);
-    for (i=0;i<delay_compensation*st->channels;i++)
-       pcm_buf[i] = st->delay_buffer[(st->encoder_buffer-delay_compensation)*st->channels+i];
+    ALLOC(pcm_buf, (total_buffer+frame_size)*st->channels, opus_val16);
+    for (i=0;i<total_buffer*st->channels;i++)
+       pcm_buf[i] = st->delay_buffer[(st->encoder_buffer-total_buffer)*st->channels+i];
 
     if (st->mode == MODE_CELT_ONLY)
        hp_freq_smth1 = silk_LSHIFT( silk_lin2log( VARIABLE_HP_MIN_CUTOFF_HZ ), 8 );
@@ -837,12 +904,26 @@ opus_int32 opus_encode_float(OpusEncoder *st, const opus_val16 *pcm, int frame_s
 
     if (st->application == OPUS_APPLICATION_VOIP)
     {
-       hp_cutoff(pcm, cutoff_Hz, &pcm_buf[delay_compensation*st->channels], st->hp_mem, frame_size, st->channels, st->Fs);
+       hp_cutoff(pcm, cutoff_Hz, &pcm_buf[total_buffer*st->channels], st->hp_mem, frame_size, st->channels, st->Fs);
     } else {
-       for (i=0;i<frame_size*st->channels;i++)
-          pcm_buf[delay_compensation*st->channels + i] = pcm[i];
+       dc_reject(pcm, 3, &pcm_buf[total_buffer*st->channels], st->hp_mem, frame_size, st->channels, st->Fs);
     }
 
+#ifndef FIXED_POINT
+    if (perform_analysis)
+    {
+       int nb_analysis_frames;
+       nb_analysis_frames = frame_size/(st->Fs/100);
+       for (i=0;i<nb_analysis_frames;i++)
+          tonality_analysis(&st->analysis, &analysis_info, celt_enc, pcm_buf+i*(st->Fs/100)*st->channels, st->channels);
+       if (st->signal_type == OPUS_AUTO)
+          st->voice_ratio = floor(.5+100*(1-analysis_info.music_prob));
+    } else {
+       analysis_info.valid = 0;
+       st->voice_ratio = -1;
+    }
+#endif
+
     /* SILK processing */
     if (st->mode != MODE_CELT_ONLY)
     {
@@ -948,10 +1029,10 @@ opus_int32 opus_encode_float(OpusEncoder *st, const opus_val16 *pcm, int frame_s
         }
 
 #ifdef FIXED_POINT
-        pcm_silk = pcm_buf+delay_compensation*st->channels;
+        pcm_silk = pcm_buf+total_buffer*st->channels;
 #else
         for (i=0;i<frame_size*st->channels;i++)
-            pcm_silk[i] = FLOAT2INT16(pcm_buf[delay_compensation*st->channels + i]);
+            pcm_silk[i] = FLOAT2INT16(pcm_buf[total_buffer*st->channels + i]);
 #endif
         ret = silk_Encode( silk_enc, &st->silk_mode, pcm_silk, frame_size, &enc, &nBytes, 0 );
         if( ret ) {
@@ -1052,13 +1133,13 @@ opus_int32 opus_encode_float(OpusEncoder *st, const opus_val16 *pcm, int frame_s
     if (st->mode != MODE_SILK_ONLY && st->mode != st->prev_mode && st->prev_mode > 0)
     {
        for (i=0;i<st->channels*st->Fs/400;i++)
-          tmp_prefill[i] = st->delay_buffer[(st->encoder_buffer-st->delay_compensation-st->Fs/400)*st->channels + i];
+          tmp_prefill[i] = st->delay_buffer[(st->encoder_buffer-total_buffer-st->Fs/400)*st->channels + i];
     }
 
-    for (i=0;i<st->channels*(st->encoder_buffer-(frame_size+delay_compensation));i++)
+    for (i=0;i<st->channels*(st->encoder_buffer-(frame_size+total_buffer));i++)
         st->delay_buffer[i] = st->delay_buffer[i+st->channels*frame_size];
     for (;i<st->encoder_buffer*st->channels;i++)
-        st->delay_buffer[i] = pcm_buf[(frame_size+delay_compensation-st->encoder_buffer)*st->channels+i];
+        st->delay_buffer[i] = pcm_buf[(frame_size+total_buffer-st->encoder_buffer)*st->channels+i];
 
 
     if (st->mode != MODE_HYBRID || st->stream_channels==1)
@@ -1079,7 +1160,7 @@ opus_int32 opus_encode_float(OpusEncoder *st, const opus_val16 *pcm, int frame_s
             g1 *= (1./16384);
             g2 *= (1./16384);
 #endif
-            stereo_fade(pcm_buf, pcm_buf, g1, g2, celt_mode->overlap,
+            stereo_fade(pcm_buf+extra_buffer*st->channels, pcm_buf+extra_buffer*st->channels, g1, g2, celt_mode->overlap,
                   frame_size, st->channels, celt_mode->window, st->Fs);
             st->hybrid_stereo_width_Q14 = st->silk_mode.stereoWidth_Q14;
         }
@@ -1131,7 +1212,7 @@ opus_int32 opus_encode_float(OpusEncoder *st, const opus_val16 *pcm, int frame_s
         int err;
         celt_encoder_ctl(celt_enc, CELT_SET_START_BAND(0));
         celt_encoder_ctl(celt_enc, OPUS_SET_VBR(0));
-        err = celt_encode_with_ec(celt_enc, pcm_buf, st->Fs/200, data+nb_compr_bytes, redundancy_bytes, NULL);
+        err = celt_encode_with_ec(celt_enc, pcm_buf+extra_buffer*st->channels, st->Fs/200, data+nb_compr_bytes, redundancy_bytes, NULL);
         if (err < 0)
         {
            RESTORE_STACK;
@@ -1157,7 +1238,11 @@ opus_int32 opus_encode_float(OpusEncoder *st, const opus_val16 *pcm, int frame_s
         /* If false, we already busted the budget and we'll end up with a "PLC packet" */
         if (ec_tell(&enc) <= 8*nb_compr_bytes)
         {
-           ret = celt_encode_with_ec(celt_enc, pcm_buf, frame_size, NULL, nb_compr_bytes, &enc);
+#ifndef FIXED_POINT
+           if (perform_analysis)
+              celt_encoder_ctl(celt_enc, CELT_SET_ANALYSIS(&analysis_info));
+#endif
+           ret = celt_encode_with_ec(celt_enc, pcm_buf+extra_buffer*st->channels, frame_size, NULL, nb_compr_bytes, &enc);
            if (ret < 0)
            {
               RESTORE_STACK;
@@ -1180,9 +1265,9 @@ opus_int32 opus_encode_float(OpusEncoder *st, const opus_val16 *pcm, int frame_s
         celt_encoder_ctl(celt_enc, CELT_SET_PREDICTION(0));
 
         /* NOTE: We could speed this up slightly (at the expense of code size) by just adding a function that prefills the buffer */
-        celt_encode_with_ec(celt_enc, pcm_buf+st->channels*(frame_size-N2-N4), N4, dummy, 2, NULL);
+        celt_encode_with_ec(celt_enc, pcm_buf+st->channels*(extra_buffer+frame_size-N2-N4), N4, dummy, 2, NULL);
 
-        err = celt_encode_with_ec(celt_enc, pcm_buf+st->channels*(frame_size-N2), N2, data+nb_compr_bytes, redundancy_bytes, NULL);
+        err = celt_encode_with_ec(celt_enc, pcm_buf+st->channels*(extra_buffer+frame_size-N2), N2, data+nb_compr_bytes, redundancy_bytes, NULL);
         if (err < 0)
         {
            RESTORE_STACK;

src/tansig_table.h

+/* This file is auto-generated by gen_tables */
+
+static const opus_val16 tansig_table[501] = {
+-1.000000, -1.000000, -1.000000, -1.000000, -1.000000, 
+-1.000000, -1.000000, -1.000000, -1.000000, -1.000000, 
+-1.000000, -1.000000, -1.000000, -1.000000, -1.000000, 
+-1.000000, -1.000000, -1.000000, -1.000000, -1.000000, 
+-1.000000, -1.000000, -1.000000, -1.000000, -1.000000, 
+-1.000000, -1.000000, -1.000000, -1.000000, -1.000000, 
+-1.000000, -1.000000, -1.000000, -1.000000, -1.000000, 
+-1.000000, -1.000000, -1.000000, -1.000000, -1.000000, 
+-1.000000, -1.000000, -1.000000, -1.000000, -1.000000, 
+-1.000000, -1.000000, -1.000000, -1.000000, -1.000000, 
+-1.000000, -1.000000, -1.000000, -1.000000, -1.000000, 
+-1.000000, -1.000000, -1.000000, -1.000000, -1.000000, 
+-1.000000, -0.999999, -0.999999, -0.999999, -0.999999, 
+-0.999999, -0.999999, -0.999999, -0.999999, -0.999999, 
+-0.999999, -0.999999, -0.999999, -0.999999, -0.999998, 
+-0.999998, -0.999998, -0.999998, -0.999998, -0.999998, 
+-0.999997, -0.999997, -0.999997, -0.999997, -0.999997, 
+-0.999996, -0.999996, -0.999996, -0.999995, -0.999995, 
+-0.999994, -0.999994, -0.999994, -0.999993, -0.999992, 
+-0.999992, -0.999991, -0.999990, -0.999990, -0.999989, 
+-0.999988, -0.999987, -0.999986, -0.999984, -0.999983, 
+-0.999982, -0.999980, -0.999978, -0.999977, -0.999975, 
+-0.999973, -0.999970, -0.999968, -0.999965, -0.999962, 
+-0.999959, -0.999956, -0.999952, -0.999948, -0.999944, 
+-0.999939, -0.999934, -0.999929, -0.999923, -0.999916, 
+-0.999909, -0.999902, -0.999893, -0.999885, -0.999875, 
+-0.999865, -0.999853, -0.999841, -0.999828, -0.999813, 
+-0.999798, -0.999781, -0.999763, -0.999743, -0.999722, 
+-0.999699, -0.999673, -0.999646, -0.999617, -0.999585, 
+-0.999550, -0.999513, -0.999472, -0.999428, -0.999381, 
+-0.999329, -0.999273, -0.999213, -0.999147, -0.999076, 
+-0.999000, -0.998916, -0.998826, -0.998728, -0.998623, 
+-0.998508, -0.998384, -0.998249, -0.998104, -0.997946, 
+-0.997775, -0.997590, -0.997389, -0.997172, -0.996937, 
+-0.996682, -0.996407, -0.996108, -0.995784, -0.995434, 
+-0.995055, -0.994644, -0.994199, -0.993718, -0.993196, 
+-0.992631, -0.992020, -0.991359, -0.990642, -0.989867, 
+-0.989027, -0.988119, -0.987136, -0.986072, -0.984921, 
+-0.983675, -0.982327, -0.980869, -0.979293, -0.977587, 
+-0.975743, -0.973749, -0.971594, -0.969265, -0.966747, 
+-0.964028, -0.961090, -0.957917, -0.954492, -0.950795, 
+-0.946806, -0.942503, -0.937863, -0.932862, -0.927473, 
+-0.921669, -0.915420, -0.908698, -0.901468, -0.893698, 
+-0.885352, -0.876393, -0.866784, -0.856485, -0.845456, 
+-0.833655, -0.821040, -0.807569, -0.793199, -0.777888, 
+-0.761594, -0.744277, -0.725897, -0.706419, -0.685809, 
+-0.664037, -0.641077, -0.616909, -0.591519, -0.564900, 
+-0.537050, -0.507977, -0.477700, -0.446244, -0.413644, 
+-0.379949, -0.345214, -0.309507, -0.272905, -0.235496, 
+-0.197375, -0.158649, -0.119427, -0.079830, -0.039979, 
+0.000000, 0.039979, 0.079830, 0.119427, 0.158649, 
+0.197375, 0.235496, 0.272905, 0.309507, 0.345214, 
+0.379949, 0.413644, 0.446244, 0.477700, 0.507977, 
+0.537050, 0.564900, 0.591519, 0.616909, 0.641077, 
+0.664037, 0.685809, 0.706419, 0.725897, 0.744277, 
+0.761594, 0.777888, 0.793199, 0.807569, 0.821040, 
+0.833655, 0.845456, 0.856485, 0.866784, 0.876393, 
+0.885352, 0.893698, 0.901468, 0.908698, 0.915420, 
+0.921669, 0.927473, 0.932862, 0.937863, 0.942503, 
+0.946806, 0.950795, 0.954492, 0.957917, 0.961090, 
+0.964028, 0.966747, 0.969265, 0.971594, 0.973749, 
+0.975743, 0.977587, 0.979293, 0.980869, 0.982327, 
+0.983675, 0.984921, 0.986072, 0.987136, 0.988119, 
+0.989027, 0.989867, 0.990642, 0.991359, 0.992020, 
+0.992631, 0.993196, 0.993718, 0.994199, 0.994644, 
+0.995055, 0.995434, 0.995784, 0.996108, 0.996407, 
+0.996682, 0.996937, 0.997172, 0.997389, 0.997590, 
+0.997775, 0.997946, 0.998104, 0.998249, 0.998384, 
+0.998508, 0.998623, 0.998728, 0.998826, 0.998916, 
+0.999000, 0.999076, 0.999147, 0.999213, 0.999273, 
+0.999329, 0.999381, 0.999428, 0.999472, 0.999513, 
+0.999550, 0.999585, 0.999617, 0.999646, 0.999673, 
+0.999699, 0.999722, 0.999743, 0.999763, 0.999781, 
+0.999798, 0.999813, 0.999828, 0.999841, 0.999853, 
+0.999865, 0.999875, 0.999885, 0.999893, 0.999902, 
+0.999909, 0.999916, 0.999923, 0.999929, 0.999934, 
+0.999939, 0.999944, 0.999948, 0.999952, 0.999956, 
+0.999959, 0.999962, 0.999965, 0.999968, 0.999970, 
+0.999973, 0.999975, 0.999977, 0.999978, 0.999980, 
+0.999982, 0.999983, 0.999984, 0.999986, 0.999987, 
+0.999988, 0.999989, 0.999990, 0.999990, 0.999991, 
+0.999992, 0.999992, 0.999993, 0.999994, 0.999994, 
+0.999994, 0.999995, 0.999995, 0.999996, 0.999996, 
+0.999996, 0.999997, 0.999997, 0.999997, 0.999997, 
+0.999997, 0.999998, 0.999998, 0.999998, 0.999998, 
+0.999998, 0.999998, 0.999999, 0.999999, 0.999999, 
+0.999999, 0.999999, 0.999999, 0.999999, 0.999999, 
+0.999999, 0.999999, 0.999999, 0.999999, 0.999999, 
+1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+1.000000, 1.000000, 1.000000, 1.000000, 1.000000, 
+1.000000, 
+};