Home | History | Annotate | Download | only in src
      1 /* Copyright (c) 2011 Xiph.Org Foundation
      2    Written by Jean-Marc Valin */
      3 /*
      4    Redistribution and use in source and binary forms, with or without
      5    modification, are permitted provided that the following conditions
      6    are met:
      7 
      8    - Redistributions of source code must retain the above copyright
      9    notice, this list of conditions and the following disclaimer.
     10 
     11    - Redistributions in binary form must reproduce the above copyright
     12    notice, this list of conditions and the following disclaimer in the
     13    documentation and/or other materials provided with the distribution.
     14 
     15    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
     16    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
     17    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
     18    A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
     19    CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
     20    EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
     21    PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
     22    PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
     23    LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
     24    NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
     25    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
     26 */
     27 
     28 #ifdef HAVE_CONFIG_H
     29 #include "config.h"
     30 #endif
     31 
     32 #include "kiss_fft.h"
     33 #include "celt.h"
     34 #include "modes.h"
     35 #include "arch.h"
     36 #include "quant_bands.h"
     37 #include <stdio.h>
     38 #include "analysis.h"
     39 #include "mlp.h"
     40 #include "stack_alloc.h"
     41 
     42 extern const MLP net;
     43 
     44 #ifndef M_PI
     45 #define M_PI 3.141592653
     46 #endif
     47 
     48 static const float dct_table[128] = {
     49         0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f,
     50         0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f,
     51         0.351851f, 0.338330f, 0.311806f, 0.273300f, 0.224292f, 0.166664f, 0.102631f, 0.034654f,
     52        -0.034654f,-0.102631f,-0.166664f,-0.224292f,-0.273300f,-0.311806f,-0.338330f,-0.351851f,
     53         0.346760f, 0.293969f, 0.196424f, 0.068975f,-0.068975f,-0.196424f,-0.293969f,-0.346760f,
     54        -0.346760f,-0.293969f,-0.196424f,-0.068975f, 0.068975f, 0.196424f, 0.293969f, 0.346760f,
     55         0.338330f, 0.224292f, 0.034654f,-0.166664f,-0.311806f,-0.351851f,-0.273300f,-0.102631f,
     56         0.102631f, 0.273300f, 0.351851f, 0.311806f, 0.166664f,-0.034654f,-0.224292f,-0.338330f,
     57         0.326641f, 0.135299f,-0.135299f,-0.326641f,-0.326641f,-0.135299f, 0.135299f, 0.326641f,
     58         0.326641f, 0.135299f,-0.135299f,-0.326641f,-0.326641f,-0.135299f, 0.135299f, 0.326641f,
     59         0.311806f, 0.034654f,-0.273300f,-0.338330f,-0.102631f, 0.224292f, 0.351851f, 0.166664f,
     60        -0.166664f,-0.351851f,-0.224292f, 0.102631f, 0.338330f, 0.273300f,-0.034654f,-0.311806f,
     61         0.293969f,-0.068975f,-0.346760f,-0.196424f, 0.196424f, 0.346760f, 0.068975f,-0.293969f,
     62        -0.293969f, 0.068975f, 0.346760f, 0.196424f,-0.196424f,-0.346760f,-0.068975f, 0.293969f,
     63         0.273300f,-0.166664f,-0.338330f, 0.034654f, 0.351851f, 0.102631f,-0.311806f,-0.224292f,
     64         0.224292f, 0.311806f,-0.102631f,-0.351851f,-0.034654f, 0.338330f, 0.166664f,-0.273300f,
     65 };
     66 
     67 static const float analysis_window[240] = {
     68       0.000043f, 0.000171f, 0.000385f, 0.000685f, 0.001071f, 0.001541f, 0.002098f, 0.002739f,
     69       0.003466f, 0.004278f, 0.005174f, 0.006156f, 0.007222f, 0.008373f, 0.009607f, 0.010926f,
     70       0.012329f, 0.013815f, 0.015385f, 0.017037f, 0.018772f, 0.020590f, 0.022490f, 0.024472f,
     71       0.026535f, 0.028679f, 0.030904f, 0.033210f, 0.035595f, 0.038060f, 0.040604f, 0.043227f,
     72       0.045928f, 0.048707f, 0.051564f, 0.054497f, 0.057506f, 0.060591f, 0.063752f, 0.066987f,
     73       0.070297f, 0.073680f, 0.077136f, 0.080665f, 0.084265f, 0.087937f, 0.091679f, 0.095492f,
     74       0.099373f, 0.103323f, 0.107342f, 0.111427f, 0.115579f, 0.119797f, 0.124080f, 0.128428f,
     75       0.132839f, 0.137313f, 0.141849f, 0.146447f, 0.151105f, 0.155823f, 0.160600f, 0.165435f,
     76       0.170327f, 0.175276f, 0.180280f, 0.185340f, 0.190453f, 0.195619f, 0.200838f, 0.206107f,
     77       0.211427f, 0.216797f, 0.222215f, 0.227680f, 0.233193f, 0.238751f, 0.244353f, 0.250000f,
     78       0.255689f, 0.261421f, 0.267193f, 0.273005f, 0.278856f, 0.284744f, 0.290670f, 0.296632f,
     79       0.302628f, 0.308658f, 0.314721f, 0.320816f, 0.326941f, 0.333097f, 0.339280f, 0.345492f,
     80       0.351729f, 0.357992f, 0.364280f, 0.370590f, 0.376923f, 0.383277f, 0.389651f, 0.396044f,
     81       0.402455f, 0.408882f, 0.415325f, 0.421783f, 0.428254f, 0.434737f, 0.441231f, 0.447736f,
     82       0.454249f, 0.460770f, 0.467298f, 0.473832f, 0.480370f, 0.486912f, 0.493455f, 0.500000f,
     83       0.506545f, 0.513088f, 0.519630f, 0.526168f, 0.532702f, 0.539230f, 0.545751f, 0.552264f,
     84       0.558769f, 0.565263f, 0.571746f, 0.578217f, 0.584675f, 0.591118f, 0.597545f, 0.603956f,
     85       0.610349f, 0.616723f, 0.623077f, 0.629410f, 0.635720f, 0.642008f, 0.648271f, 0.654508f,
     86       0.660720f, 0.666903f, 0.673059f, 0.679184f, 0.685279f, 0.691342f, 0.697372f, 0.703368f,
     87       0.709330f, 0.715256f, 0.721144f, 0.726995f, 0.732807f, 0.738579f, 0.744311f, 0.750000f,
     88       0.755647f, 0.761249f, 0.766807f, 0.772320f, 0.777785f, 0.783203f, 0.788573f, 0.793893f,
     89       0.799162f, 0.804381f, 0.809547f, 0.814660f, 0.819720f, 0.824724f, 0.829673f, 0.834565f,
     90       0.839400f, 0.844177f, 0.848895f, 0.853553f, 0.858151f, 0.862687f, 0.867161f, 0.871572f,
     91       0.875920f, 0.880203f, 0.884421f, 0.888573f, 0.892658f, 0.896677f, 0.900627f, 0.904508f,
     92       0.908321f, 0.912063f, 0.915735f, 0.919335f, 0.922864f, 0.926320f, 0.929703f, 0.933013f,
     93       0.936248f, 0.939409f, 0.942494f, 0.945503f, 0.948436f, 0.951293f, 0.954072f, 0.956773f,
     94       0.959396f, 0.961940f, 0.964405f, 0.966790f, 0.969096f, 0.971321f, 0.973465f, 0.975528f,
     95       0.977510f, 0.979410f, 0.981228f, 0.982963f, 0.984615f, 0.986185f, 0.987671f, 0.989074f,
     96       0.990393f, 0.991627f, 0.992778f, 0.993844f, 0.994826f, 0.995722f, 0.996534f, 0.997261f,
     97       0.997902f, 0.998459f, 0.998929f, 0.999315f, 0.999615f, 0.999829f, 0.999957f, 1.000000f,
     98 };
     99 
    100 static const int tbands[NB_TBANDS+1] = {
    101        2,  4,  6,  8, 10, 12, 14, 16, 20, 24, 28, 32, 40, 48, 56, 68, 80, 96, 120
    102 };
    103 
    104 static const int extra_bands[NB_TOT_BANDS+1] = {
    105       1, 2,  4,  6,  8, 10, 12, 14, 16, 20, 24, 28, 32, 40, 48, 56, 68, 80, 96, 120, 160, 200
    106 };
    107 
    108 /*static const float tweight[NB_TBANDS+1] = {
    109       .3, .4, .5, .6, .7, .8, .9, 1., 1., 1., 1., 1., 1., 1., .8, .7, .6, .5
    110 };*/
    111 
    112 #define NB_TONAL_SKIP_BANDS 9
    113 
    114 #define cA 0.43157974f
    115 #define cB 0.67848403f
    116 #define cC 0.08595542f
    117 #define cE ((float)M_PI/2)
    118 static inline float fast_atan2f(float y, float x) {
    119    float x2, y2;
    120    /* Should avoid underflow on the values we'll get */
    121    if (ABS16(x)+ABS16(y)<1e-9f)
    122    {
    123       x*=1e12f;
    124       y*=1e12f;
    125    }
    126    x2 = x*x;
    127    y2 = y*y;
    128    if(x2<y2){
    129       float den = (y2 + cB*x2) * (y2 + cC*x2);
    130       if (den!=0)
    131          return -x*y*(y2 + cA*x2) / den + (y<0 ? -cE : cE);
    132       else
    133          return (y<0 ? -cE : cE);
    134    }else{
    135       float den = (x2 + cB*y2) * (x2 + cC*y2);
    136       if (den!=0)
    137          return  x*y*(x2 + cA*y2) / den + (y<0 ? -cE : cE) - (x*y<0 ? -cE : cE);
    138       else
    139          return (y<0 ? -cE : cE) - (x*y<0 ? -cE : cE);
    140    }
    141 }
    142 
    143 void tonality_get_info(TonalityAnalysisState *tonal, AnalysisInfo *info_out, int len)
    144 {
    145    int pos;
    146    int curr_lookahead;
    147    float psum;
    148    int i;
    149 
    150    pos = tonal->read_pos;
    151    curr_lookahead = tonal->write_pos-tonal->read_pos;
    152    if (curr_lookahead<0)
    153       curr_lookahead += DETECT_SIZE;
    154 
    155    if (len > 480 && pos != tonal->write_pos)
    156    {
    157       pos++;
    158       if (pos==DETECT_SIZE)
    159          pos=0;
    160    }
    161    if (pos == tonal->write_pos)
    162       pos--;
    163    if (pos<0)
    164       pos = DETECT_SIZE-1;
    165    OPUS_COPY(info_out, &tonal->info[pos], 1);
    166    tonal->read_subframe += len/120;
    167    while (tonal->read_subframe>=4)
    168    {
    169       tonal->read_subframe -= 4;
    170       tonal->read_pos++;
    171    }
    172    if (tonal->read_pos>=DETECT_SIZE)
    173       tonal->read_pos-=DETECT_SIZE;
    174 
    175    /* Compensate for the delay in the features themselves.
    176       FIXME: Need a better estimate the 10 I just made up */
    177    curr_lookahead = IMAX(curr_lookahead-10, 0);
    178 
    179    psum=0;
    180    /* Summing the probability of transition patterns that involve music at
    181       time (DETECT_SIZE-curr_lookahead-1) */
    182    for (i=0;i<DETECT_SIZE-curr_lookahead;i++)
    183       psum += tonal->pmusic[i];
    184    for (;i<DETECT_SIZE;i++)
    185       psum += tonal->pspeech[i];
    186    psum = psum*tonal->music_confidence + (1-psum)*tonal->speech_confidence;
    187    /*printf("%f %f %f\n", psum, info_out->music_prob, info_out->tonality);*/
    188 
    189    info_out->music_prob = psum;
    190 }
    191 
    192 void tonality_analysis(TonalityAnalysisState *tonal, AnalysisInfo *info_out, const CELTMode *celt_mode, const void *x, int len, int offset, int c1, int c2, int C, int lsb_depth, downmix_func downmix)
    193 {
    194     int i, b;
    195     const kiss_fft_state *kfft;
    196     VARDECL(kiss_fft_cpx, in);
    197     VARDECL(kiss_fft_cpx, out);
    198     int N = 480, N2=240;
    199     float * OPUS_RESTRICT A = tonal->angle;
    200     float * OPUS_RESTRICT dA = tonal->d_angle;
    201     float * OPUS_RESTRICT d2A = tonal->d2_angle;
    202     VARDECL(float, tonality);
    203     VARDECL(float, noisiness);
    204     float band_tonality[NB_TBANDS];
    205     float logE[NB_TBANDS];
    206     float BFCC[8];
    207     float features[25];
    208     float frame_tonality;
    209     float max_frame_tonality;
    210     /*float tw_sum=0;*/
    211     float frame_noisiness;
    212     const float pi4 = (float)(M_PI*M_PI*M_PI*M_PI);
    213     float slope=0;
    214     float frame_stationarity;
    215     float relativeE;
    216     float frame_probs[2];
    217     float alpha, alphaE, alphaE2;
    218     float frame_loudness;
    219     float bandwidth_mask;
    220     int bandwidth=0;
    221     float maxE = 0;
    222     float noise_floor;
    223     int remaining;
    224     AnalysisInfo *info;
    225     SAVE_STACK;
    226 
    227     tonal->last_transition++;
    228     alpha = 1.f/IMIN(20, 1+tonal->count);
    229     alphaE = 1.f/IMIN(50, 1+tonal->count);
    230     alphaE2 = 1.f/IMIN(1000, 1+tonal->count);
    231 
    232     if (tonal->count<4)
    233        tonal->music_prob = .5;
    234     kfft = celt_mode->mdct.kfft[0];
    235     if (tonal->count==0)
    236        tonal->mem_fill = 240;
    237     downmix(x, &tonal->inmem[tonal->mem_fill], IMIN(len, ANALYSIS_BUF_SIZE-tonal->mem_fill), offset, c1, c2, C);
    238     if (tonal->mem_fill+len < ANALYSIS_BUF_SIZE)
    239     {
    240        tonal->mem_fill += len;
    241        /* Don't have enough to update the analysis */
    242        RESTORE_STACK;
    243        return;
    244     }
    245     info = &tonal->info[tonal->write_pos++];
    246     if (tonal->write_pos>=DETECT_SIZE)
    247        tonal->write_pos-=DETECT_SIZE;
    248 
    249     ALLOC(in, 480, kiss_fft_cpx);
    250     ALLOC(out, 480, kiss_fft_cpx);
    251     ALLOC(tonality, 240, float);
    252     ALLOC(noisiness, 240, float);
    253     for (i=0;i<N2;i++)
    254     {
    255        float w = analysis_window[i];
    256        in[i].r = w*tonal->inmem[i];
    257        in[i].i = w*tonal->inmem[N2+i];
    258        in[N-i-1].r = w*tonal->inmem[N-i-1];
    259        in[N-i-1].i = w*tonal->inmem[N+N2-i-1];
    260     }
    261     OPUS_MOVE(tonal->inmem, tonal->inmem+ANALYSIS_BUF_SIZE-240, 240);
    262     remaining = len - (ANALYSIS_BUF_SIZE-tonal->mem_fill);
    263     downmix(x, &tonal->inmem[240], remaining, offset+ANALYSIS_BUF_SIZE-tonal->mem_fill, c1, c2, C);
    264     tonal->mem_fill = 240 + remaining;
    265     opus_fft(kfft, in, out);
    266 
    267     for (i=1;i<N2;i++)
    268     {
    269        float X1r, X2r, X1i, X2i;
    270        float angle, d_angle, d2_angle;
    271        float angle2, d_angle2, d2_angle2;
    272        float mod1, mod2, avg_mod;
    273        X1r = out[i].r+out[N-i].r;
    274        X1i = out[i].i-out[N-i].i;
    275        X2r = out[i].i+out[N-i].i;
    276        X2i = out[N-i].r-out[i].r;
    277 
    278        angle = (float)(.5f/M_PI)*fast_atan2f(X1i, X1r);
    279        d_angle = angle - A[i];
    280        d2_angle = d_angle - dA[i];
    281 
    282        angle2 = (float)(.5f/M_PI)*fast_atan2f(X2i, X2r);
    283        d_angle2 = angle2 - angle;
    284        d2_angle2 = d_angle2 - d_angle;
    285 
    286        mod1 = d2_angle - (float)floor(.5+d2_angle);
    287        noisiness[i] = ABS16(mod1);
    288        mod1 *= mod1;
    289        mod1 *= mod1;
    290 
    291        mod2 = d2_angle2 - (float)floor(.5+d2_angle2);
    292        noisiness[i] += ABS16(mod2);
    293        mod2 *= mod2;
    294        mod2 *= mod2;
    295 
    296        avg_mod = .25f*(d2A[i]+2.f*mod1+mod2);
    297        tonality[i] = 1.f/(1.f+40.f*16.f*pi4*avg_mod)-.015f;
    298 
    299        A[i] = angle2;
    300        dA[i] = d_angle2;
    301        d2A[i] = mod2;
    302     }
    303 
    304     frame_tonality = 0;
    305     max_frame_tonality = 0;
    306     /*tw_sum = 0;*/
    307     info->activity = 0;
    308     frame_noisiness = 0;
    309     frame_stationarity = 0;
    310     if (!tonal->count)
    311     {
    312        for (b=0;b<NB_TBANDS;b++)
    313        {
    314           tonal->lowE[b] = 1e10;
    315           tonal->highE[b] = -1e10;
    316        }
    317     }
    318     relativeE = 0;
    319     frame_loudness = 0;
    320     bandwidth_mask = 0;
    321     for (b=0;b<NB_TBANDS;b++)
    322     {
    323        float E=0, tE=0, nE=0;
    324        float L1, L2;
    325        float stationarity;
    326        for (i=tbands[b];i<tbands[b+1];i++)
    327        {
    328           float binE = out[i].r*(float)out[i].r + out[N-i].r*(float)out[N-i].r
    329                      + out[i].i*(float)out[i].i + out[N-i].i*(float)out[N-i].i;
    330 #ifdef FIXED_POINT
    331           /* FIXME: It's probably best to change the BFCC filter initial state instead */
    332           binE *= 5.55e-17f;
    333 #endif
    334           E += binE;
    335           tE += binE*tonality[i];
    336           nE += binE*2.f*(.5f-noisiness[i]);
    337        }
    338        tonal->E[tonal->E_count][b] = E;
    339        frame_noisiness += nE/(1e-15f+E);
    340 
    341        frame_loudness += sqrt(E+1e-10f);
    342        logE[b] = (float)log(E+1e-10f);
    343        tonal->lowE[b] = MIN32(logE[b], tonal->lowE[b]+.01f);
    344        tonal->highE[b] = MAX32(logE[b], tonal->highE[b]-.1f);
    345        if (tonal->highE[b] < tonal->lowE[b]+1.f)
    346        {
    347           tonal->highE[b]+=.5f;
    348           tonal->lowE[b]-=.5f;
    349        }
    350        relativeE += (logE[b]-tonal->lowE[b])/(1e-15+tonal->highE[b]-tonal->lowE[b]);
    351 
    352        L1=L2=0;
    353        for (i=0;i<NB_FRAMES;i++)
    354        {
    355           L1 += sqrt(tonal->E[i][b]);
    356           L2 += tonal->E[i][b];
    357        }
    358 
    359        stationarity = MIN16(0.99f,L1/sqrt(1e-15+NB_FRAMES*L2));
    360        stationarity *= stationarity;
    361        stationarity *= stationarity;
    362        frame_stationarity += stationarity;
    363        /*band_tonality[b] = tE/(1e-15+E)*/;
    364        band_tonality[b] = MAX16(tE/(1e-15+E), stationarity*tonal->prev_band_tonality[b]);
    365 #if 0
    366        if (b>=NB_TONAL_SKIP_BANDS)
    367        {
    368           frame_tonality += tweight[b]*band_tonality[b];
    369           tw_sum += tweight[b];
    370        }
    371 #else
    372        frame_tonality += band_tonality[b];
    373        if (b>=NB_TBANDS-NB_TONAL_SKIP_BANDS)
    374           frame_tonality -= band_tonality[b-NB_TBANDS+NB_TONAL_SKIP_BANDS];
    375 #endif
    376        max_frame_tonality = MAX16(max_frame_tonality, (1.f+.03f*(b-NB_TBANDS))*frame_tonality);
    377        slope += band_tonality[b]*(b-8);
    378        /*printf("%f %f ", band_tonality[b], stationarity);*/
    379        tonal->prev_band_tonality[b] = band_tonality[b];
    380     }
    381 
    382     bandwidth_mask = 0;
    383     bandwidth = 0;
    384     maxE = 0;
    385     noise_floor = 5.7e-4f/(1<<(IMAX(0,lsb_depth-8)));
    386 #ifdef FIXED_POINT
    387     noise_floor *= 1<<(15+SIG_SHIFT);
    388 #endif
    389     noise_floor *= noise_floor;
    390     for (b=0;b<NB_TOT_BANDS;b++)
    391     {
    392        float E=0;
    393        int band_start, band_end;
    394        /* Keep a margin of 300 Hz for aliasing */
    395        band_start = extra_bands[b];
    396        band_end = extra_bands[b+1];
    397        for (i=band_start;i<band_end;i++)
    398        {
    399           float binE = out[i].r*(float)out[i].r + out[N-i].r*(float)out[N-i].r
    400                      + out[i].i*(float)out[i].i + out[N-i].i*(float)out[N-i].i;
    401           E += binE;
    402        }
    403        maxE = MAX32(maxE, E);
    404        tonal->meanE[b] = MAX32((1-alphaE2)*tonal->meanE[b], E);
    405        E = MAX32(E, tonal->meanE[b]);
    406        /* Use a simple follower with 13 dB/Bark slope for spreading function */
    407        bandwidth_mask = MAX32(.05f*bandwidth_mask, E);
    408        /* Consider the band "active" only if all these conditions are met:
    409           1) less than 10 dB below the simple follower
    410           2) less than 90 dB below the peak band (maximal masking possible considering
    411              both the ATH and the loudness-dependent slope of the spreading function)
    412           3) above the PCM quantization noise floor
    413        */
    414        if (E>.1*bandwidth_mask && E*1e9f > maxE && E > noise_floor*(band_end-band_start))
    415           bandwidth = b;
    416     }
    417     if (tonal->count<=2)
    418        bandwidth = 20;
    419     frame_loudness = 20*(float)log10(frame_loudness);
    420     tonal->Etracker = MAX32(tonal->Etracker-.03f, frame_loudness);
    421     tonal->lowECount *= (1-alphaE);
    422     if (frame_loudness < tonal->Etracker-30)
    423        tonal->lowECount += alphaE;
    424 
    425     for (i=0;i<8;i++)
    426     {
    427        float sum=0;
    428        for (b=0;b<16;b++)
    429           sum += dct_table[i*16+b]*logE[b];
    430        BFCC[i] = sum;
    431     }
    432 
    433     frame_stationarity /= NB_TBANDS;
    434     relativeE /= NB_TBANDS;
    435     if (tonal->count<10)
    436        relativeE = .5;
    437     frame_noisiness /= NB_TBANDS;
    438 #if 1
    439     info->activity = frame_noisiness + (1-frame_noisiness)*relativeE;
    440 #else
    441     info->activity = .5*(1+frame_noisiness-frame_stationarity);
    442 #endif
    443     frame_tonality = (max_frame_tonality/(NB_TBANDS-NB_TONAL_SKIP_BANDS));
    444     frame_tonality = MAX16(frame_tonality, tonal->prev_tonality*.8f);
    445     tonal->prev_tonality = frame_tonality;
    446 
    447     slope /= 8*8;
    448     info->tonality_slope = slope;
    449 
    450     tonal->E_count = (tonal->E_count+1)%NB_FRAMES;
    451     tonal->count++;
    452     info->tonality = frame_tonality;
    453 
    454     for (i=0;i<4;i++)
    455        features[i] = -0.12299f*(BFCC[i]+tonal->mem[i+24]) + 0.49195f*(tonal->mem[i]+tonal->mem[i+16]) + 0.69693f*tonal->mem[i+8] - 1.4349f*tonal->cmean[i];
    456 
    457     for (i=0;i<4;i++)
    458        tonal->cmean[i] = (1-alpha)*tonal->cmean[i] + alpha*BFCC[i];
    459 
    460     for (i=0;i<4;i++)
    461         features[4+i] = 0.63246f*(BFCC[i]-tonal->mem[i+24]) + 0.31623f*(tonal->mem[i]-tonal->mem[i+16]);
    462     for (i=0;i<3;i++)
    463         features[8+i] = 0.53452f*(BFCC[i]+tonal->mem[i+24]) - 0.26726f*(tonal->mem[i]+tonal->mem[i+16]) -0.53452f*tonal->mem[i+8];
    464 
    465     if (tonal->count > 5)
    466     {
    467        for (i=0;i<9;i++)
    468           tonal->std[i] = (1-alpha)*tonal->std[i] + alpha*features[i]*features[i];
    469     }
    470 
    471     for (i=0;i<8;i++)
    472     {
    473        tonal->mem[i+24] = tonal->mem[i+16];
    474        tonal->mem[i+16] = tonal->mem[i+8];
    475        tonal->mem[i+8] = tonal->mem[i];
    476        tonal->mem[i] = BFCC[i];
    477     }
    478     for (i=0;i<9;i++)
    479        features[11+i] = sqrt(tonal->std[i]);
    480     features[20] = info->tonality;
    481     features[21] = info->activity;
    482     features[22] = frame_stationarity;
    483     features[23] = info->tonality_slope;
    484     features[24] = tonal->lowECount;
    485 
    486 #ifndef DISABLE_FLOAT_API
    487     mlp_process(&net, features, frame_probs);
    488     frame_probs[0] = .5f*(frame_probs[0]+1);
    489     /* Curve fitting between the MLP probability and the actual probability */
    490     frame_probs[0] = .01f + 1.21f*frame_probs[0]*frame_probs[0] - .23f*(float)pow(frame_probs[0], 10);
    491     /* Probability of active audio (as opposed to silence) */
    492     frame_probs[1] = .5f*frame_probs[1]+.5f;
    493     /* Consider that silence has a 50-50 probability. */
    494     frame_probs[0] = frame_probs[1]*frame_probs[0] + (1-frame_probs[1])*.5f;
    495 
    496     /*printf("%f %f ", frame_probs[0], frame_probs[1]);*/
    497     {
    498        /* Probability of state transition */
    499        float tau;
    500        /* Represents independence of the MLP probabilities, where
    501           beta=1 means fully independent. */
    502        float beta;
    503        /* Denormalized probability of speech (p0) and music (p1) after update */
    504        float p0, p1;
    505        /* Probabilities for "all speech" and "all music" */
    506        float s0, m0;
    507        /* Probability sum for renormalisation */
    508        float psum;
    509        /* Instantaneous probability of speech and music, with beta pre-applied. */
    510        float speech0;
    511        float music0;
    512 
    513        /* One transition every 3 minutes of active audio */
    514        tau = .00005f*frame_probs[1];
    515        beta = .05f;
    516        if (1) {
    517           /* Adapt beta based on how "unexpected" the new prob is */
    518           float p, q;
    519           p = MAX16(.05f,MIN16(.95f,frame_probs[0]));
    520           q = MAX16(.05f,MIN16(.95f,tonal->music_prob));
    521           beta = .01f+.05f*ABS16(p-q)/(p*(1-q)+q*(1-p));
    522        }
    523        /* p0 and p1 are the probabilities of speech and music at this frame
    524           using only information from previous frame and applying the
    525           state transition model */
    526        p0 = (1-tonal->music_prob)*(1-tau) +    tonal->music_prob *tau;
    527        p1 =    tonal->music_prob *(1-tau) + (1-tonal->music_prob)*tau;
    528        /* We apply the current probability with exponent beta to work around
    529           the fact that the probability estimates aren't independent. */
    530        p0 *= (float)pow(1-frame_probs[0], beta);
    531        p1 *= (float)pow(frame_probs[0], beta);
    532        /* Normalise the probabilities to get the Marokv probability of music. */
    533        tonal->music_prob = p1/(p0+p1);
    534        info->music_prob = tonal->music_prob;
    535 
    536        /* This chunk of code deals with delayed decision. */
    537        psum=1e-20f;
    538        /* Instantaneous probability of speech and music, with beta pre-applied. */
    539        speech0 = (float)pow(1-frame_probs[0], beta);
    540        music0  = (float)pow(frame_probs[0], beta);
    541        if (tonal->count==1)
    542        {
    543           tonal->pspeech[0]=.5;
    544           tonal->pmusic [0]=.5;
    545        }
    546        /* Updated probability of having only speech (s0) or only music (m0),
    547           before considering the new observation. */
    548        s0 = tonal->pspeech[0] + tonal->pspeech[1];
    549        m0 = tonal->pmusic [0] + tonal->pmusic [1];
    550        /* Updates s0 and m0 with instantaneous probability. */
    551        tonal->pspeech[0] = s0*(1-tau)*speech0;
    552        tonal->pmusic [0] = m0*(1-tau)*music0;
    553        /* Propagate the transition probabilities */
    554        for (i=1;i<DETECT_SIZE-1;i++)
    555        {
    556           tonal->pspeech[i] = tonal->pspeech[i+1]*speech0;
    557           tonal->pmusic [i] = tonal->pmusic [i+1]*music0;
    558        }
    559        /* Probability that the latest frame is speech, when all the previous ones were music. */
    560        tonal->pspeech[DETECT_SIZE-1] = m0*tau*speech0;
    561        /* Probability that the latest frame is music, when all the previous ones were speech. */
    562        tonal->pmusic [DETECT_SIZE-1] = s0*tau*music0;
    563 
    564        /* Renormalise probabilities to 1 */
    565        for (i=0;i<DETECT_SIZE;i++)
    566           psum += tonal->pspeech[i] + tonal->pmusic[i];
    567        psum = 1.f/psum;
    568        for (i=0;i<DETECT_SIZE;i++)
    569        {
    570           tonal->pspeech[i] *= psum;
    571           tonal->pmusic [i] *= psum;
    572        }
    573        psum = tonal->pmusic[0];
    574        for (i=1;i<DETECT_SIZE;i++)
    575           psum += tonal->pspeech[i];
    576 
    577        /* Estimate our confidence in the speech/music decisions */
    578        if (frame_probs[1]>.75)
    579        {
    580           if (tonal->music_prob>.9)
    581           {
    582              float adapt;
    583              adapt = 1.f/(++tonal->music_confidence_count);
    584              tonal->music_confidence_count = IMIN(tonal->music_confidence_count, 500);
    585              tonal->music_confidence += adapt*MAX16(-.2f,frame_probs[0]-tonal->music_confidence);
    586           }
    587           if (tonal->music_prob<.1)
    588           {
    589              float adapt;
    590              adapt = 1.f/(++tonal->speech_confidence_count);
    591              tonal->speech_confidence_count = IMIN(tonal->speech_confidence_count, 500);
    592              tonal->speech_confidence += adapt*MIN16(.2f,frame_probs[0]-tonal->speech_confidence);
    593           }
    594        } else {
    595           if (tonal->music_confidence_count==0)
    596              tonal->music_confidence = .9f;
    597           if (tonal->speech_confidence_count==0)
    598              tonal->speech_confidence = .1f;
    599        }
    600        psum = MAX16(tonal->speech_confidence, MIN16(tonal->music_confidence, psum));
    601     }
    602     if (tonal->last_music != (tonal->music_prob>.5f))
    603        tonal->last_transition=0;
    604     tonal->last_music = tonal->music_prob>.5f;
    605 #else
    606     info->music_prob = 0;
    607 #endif
    608     /*for (i=0;i<25;i++)
    609        printf("%f ", features[i]);
    610     printf("\n");*/
    611 
    612     info->bandwidth = bandwidth;
    613     /*printf("%d %d\n", info->bandwidth, info->opus_bandwidth);*/
    614     info->noisiness = frame_noisiness;
    615     info->valid = 1;
    616     if (info_out!=NULL)
    617        OPUS_COPY(info_out, info, 1);
    618     RESTORE_STACK;
    619 }
    620 
    621 void run_analysis(TonalityAnalysisState *analysis, const CELTMode *celt_mode, const void *analysis_pcm,
    622                  int analysis_frame_size, int frame_size, int c1, int c2, int C, opus_int32 Fs,
    623                  int lsb_depth, downmix_func downmix, AnalysisInfo *analysis_info)
    624 {
    625    int offset;
    626    int pcm_len;
    627 
    628    if (analysis_pcm != NULL)
    629    {
    630       /* Avoid overflow/wrap-around of the analysis buffer */
    631       analysis_frame_size = IMIN((DETECT_SIZE-5)*Fs/100, analysis_frame_size);
    632 
    633       pcm_len = analysis_frame_size - analysis->analysis_offset;
    634       offset = analysis->analysis_offset;
    635       do {
    636          tonality_analysis(analysis, NULL, celt_mode, analysis_pcm, IMIN(480, pcm_len), offset, c1, c2, C, lsb_depth, downmix);
    637          offset += 480;
    638          pcm_len -= 480;
    639       } while (pcm_len>0);
    640       analysis->analysis_offset = analysis_frame_size;
    641 
    642       analysis->analysis_offset -= frame_size;
    643    }
    644 
    645    analysis_info->valid = 0;
    646    tonality_get_info(analysis, analysis_info, frame_size);
    647 }
    648