Adds functions for multistream padding/unpadding and single-stream unpadding
[opus.git] / src / analysis.c
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 OPUS_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 = (kiss_fft_scalar)(w*tonal->inmem[i]);
257        in[i].i = (kiss_fft_scalar)(w*tonal->inmem[N2+i]);
258        in[N-i-1].r = (kiss_fft_scalar)(w*tonal->inmem[N-i-1]);
259        in[N-i-1].i = (kiss_fft_scalar)(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 = (float)out[i].r+out[N-i].r;
274        X1i = (float)out[i].i-out[N-i].i;
275        X2r = (float)out[i].i+out[N-i].i;
276        X2i = (float)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     for (b=0;b<NB_TBANDS;b++)
321     {
322        float E=0, tE=0, nE=0;
323        float L1, L2;
324        float stationarity;
325        for (i=tbands[b];i<tbands[b+1];i++)
326        {
327           float binE = out[i].r*(float)out[i].r + out[N-i].r*(float)out[N-i].r
328                      + out[i].i*(float)out[i].i + out[N-i].i*(float)out[N-i].i;
329 #ifdef FIXED_POINT
330           /* FIXME: It's probably best to change the BFCC filter initial state instead */
331           binE *= 5.55e-17f;
332 #endif
333           E += binE;
334           tE += binE*tonality[i];
335           nE += binE*2.f*(.5f-noisiness[i]);
336        }
337        tonal->E[tonal->E_count][b] = E;
338        frame_noisiness += nE/(1e-15f+E);
339
340        frame_loudness += (float)sqrt(E+1e-10f);
341        logE[b] = (float)log(E+1e-10f);
342        tonal->lowE[b] = MIN32(logE[b], tonal->lowE[b]+.01f);
343        tonal->highE[b] = MAX32(logE[b], tonal->highE[b]-.1f);
344        if (tonal->highE[b] < tonal->lowE[b]+1.f)
345        {
346           tonal->highE[b]+=.5f;
347           tonal->lowE[b]-=.5f;
348        }
349        relativeE += (logE[b]-tonal->lowE[b])/(1e-15f+tonal->highE[b]-tonal->lowE[b]);
350
351        L1=L2=0;
352        for (i=0;i<NB_FRAMES;i++)
353        {
354           L1 += (float)sqrt(tonal->E[i][b]);
355           L2 += tonal->E[i][b];
356        }
357
358        stationarity = MIN16(0.99f,L1/(float)sqrt(1e-15+NB_FRAMES*L2));
359        stationarity *= stationarity;
360        stationarity *= stationarity;
361        frame_stationarity += stationarity;
362        /*band_tonality[b] = tE/(1e-15+E)*/;
363        band_tonality[b] = MAX16(tE/(1e-15f+E), stationarity*tonal->prev_band_tonality[b]);
364 #if 0
365        if (b>=NB_TONAL_SKIP_BANDS)
366        {
367           frame_tonality += tweight[b]*band_tonality[b];
368           tw_sum += tweight[b];
369        }
370 #else
371        frame_tonality += band_tonality[b];
372        if (b>=NB_TBANDS-NB_TONAL_SKIP_BANDS)
373           frame_tonality -= band_tonality[b-NB_TBANDS+NB_TONAL_SKIP_BANDS];
374 #endif
375        max_frame_tonality = MAX16(max_frame_tonality, (1.f+.03f*(b-NB_TBANDS))*frame_tonality);
376        slope += band_tonality[b]*(b-8);
377        /*printf("%f %f ", band_tonality[b], stationarity);*/
378        tonal->prev_band_tonality[b] = band_tonality[b];
379     }
380
381     bandwidth_mask = 0;
382     bandwidth = 0;
383     maxE = 0;
384     noise_floor = 5.7e-4f/(1<<(IMAX(0,lsb_depth-8)));
385 #ifdef FIXED_POINT
386     noise_floor *= 1<<(15+SIG_SHIFT);
387 #endif
388     noise_floor *= noise_floor;
389     for (b=0;b<NB_TOT_BANDS;b++)
390     {
391        float E=0;
392        int band_start, band_end;
393        /* Keep a margin of 300 Hz for aliasing */
394        band_start = extra_bands[b];
395        band_end = extra_bands[b+1];
396        for (i=band_start;i<band_end;i++)
397        {
398           float binE = out[i].r*(float)out[i].r + out[N-i].r*(float)out[N-i].r
399                      + out[i].i*(float)out[i].i + out[N-i].i*(float)out[N-i].i;
400           E += binE;
401        }
402        maxE = MAX32(maxE, E);
403        tonal->meanE[b] = MAX32((1-alphaE2)*tonal->meanE[b], E);
404        E = MAX32(E, tonal->meanE[b]);
405        /* Use a simple follower with 13 dB/Bark slope for spreading function */
406        bandwidth_mask = MAX32(.05f*bandwidth_mask, E);
407        /* Consider the band "active" only if all these conditions are met:
408           1) less than 10 dB below the simple follower
409           2) less than 90 dB below the peak band (maximal masking possible considering
410              both the ATH and the loudness-dependent slope of the spreading function)
411           3) above the PCM quantization noise floor
412        */
413        if (E>.1*bandwidth_mask && E*1e9f > maxE && E > noise_floor*(band_end-band_start))
414           bandwidth = b;
415     }
416     if (tonal->count<=2)
417        bandwidth = 20;
418     frame_loudness = 20*(float)log10(frame_loudness);
419     tonal->Etracker = MAX32(tonal->Etracker-.03f, frame_loudness);
420     tonal->lowECount *= (1-alphaE);
421     if (frame_loudness < tonal->Etracker-30)
422        tonal->lowECount += alphaE;
423
424     for (i=0;i<8;i++)
425     {
426        float sum=0;
427        for (b=0;b<16;b++)
428           sum += dct_table[i*16+b]*logE[b];
429        BFCC[i] = sum;
430     }
431
432     frame_stationarity /= NB_TBANDS;
433     relativeE /= NB_TBANDS;
434     if (tonal->count<10)
435        relativeE = .5;
436     frame_noisiness /= NB_TBANDS;
437 #if 1
438     info->activity = frame_noisiness + (1-frame_noisiness)*relativeE;
439 #else
440     info->activity = .5*(1+frame_noisiness-frame_stationarity);
441 #endif
442     frame_tonality = (max_frame_tonality/(NB_TBANDS-NB_TONAL_SKIP_BANDS));
443     frame_tonality = MAX16(frame_tonality, tonal->prev_tonality*.8f);
444     tonal->prev_tonality = frame_tonality;
445
446     slope /= 8*8;
447     info->tonality_slope = slope;
448
449     tonal->E_count = (tonal->E_count+1)%NB_FRAMES;
450     tonal->count++;
451     info->tonality = frame_tonality;
452
453     for (i=0;i<4;i++)
454        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];
455
456     for (i=0;i<4;i++)
457        tonal->cmean[i] = (1-alpha)*tonal->cmean[i] + alpha*BFCC[i];
458
459     for (i=0;i<4;i++)
460         features[4+i] = 0.63246f*(BFCC[i]-tonal->mem[i+24]) + 0.31623f*(tonal->mem[i]-tonal->mem[i+16]);
461     for (i=0;i<3;i++)
462         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];
463
464     if (tonal->count > 5)
465     {
466        for (i=0;i<9;i++)
467           tonal->std[i] = (1-alpha)*tonal->std[i] + alpha*features[i]*features[i];
468     }
469
470     for (i=0;i<8;i++)
471     {
472        tonal->mem[i+24] = tonal->mem[i+16];
473        tonal->mem[i+16] = tonal->mem[i+8];
474        tonal->mem[i+8] = tonal->mem[i];
475        tonal->mem[i] = BFCC[i];
476     }
477     for (i=0;i<9;i++)
478        features[11+i] = (float)sqrt(tonal->std[i]);
479     features[20] = info->tonality;
480     features[21] = info->activity;
481     features[22] = frame_stationarity;
482     features[23] = info->tonality_slope;
483     features[24] = tonal->lowECount;
484
485 #ifndef DISABLE_FLOAT_API
486     mlp_process(&net, features, frame_probs);
487     frame_probs[0] = .5f*(frame_probs[0]+1);
488     /* Curve fitting between the MLP probability and the actual probability */
489     frame_probs[0] = .01f + 1.21f*frame_probs[0]*frame_probs[0] - .23f*(float)pow(frame_probs[0], 10);
490     /* Probability of active audio (as opposed to silence) */
491     frame_probs[1] = .5f*frame_probs[1]+.5f;
492     /* Consider that silence has a 50-50 probability. */
493     frame_probs[0] = frame_probs[1]*frame_probs[0] + (1-frame_probs[1])*.5f;
494
495     /*printf("%f %f ", frame_probs[0], frame_probs[1]);*/
496     {
497        /* Probability of state transition */
498        float tau;
499        /* Represents independence of the MLP probabilities, where
500           beta=1 means fully independent. */
501        float beta;
502        /* Denormalized probability of speech (p0) and music (p1) after update */
503        float p0, p1;
504        /* Probabilities for "all speech" and "all music" */
505        float s0, m0;
506        /* Probability sum for renormalisation */
507        float psum;
508        /* Instantaneous probability of speech and music, with beta pre-applied. */
509        float speech0;
510        float music0;
511
512        /* One transition every 3 minutes of active audio */
513        tau = .00005f*frame_probs[1];
514        beta = .05f;
515        if (1) {
516           /* Adapt beta based on how "unexpected" the new prob is */
517           float p, q;
518           p = MAX16(.05f,MIN16(.95f,frame_probs[0]));
519           q = MAX16(.05f,MIN16(.95f,tonal->music_prob));
520           beta = .01f+.05f*ABS16(p-q)/(p*(1-q)+q*(1-p));
521        }
522        /* p0 and p1 are the probabilities of speech and music at this frame
523           using only information from previous frame and applying the
524           state transition model */
525        p0 = (1-tonal->music_prob)*(1-tau) +    tonal->music_prob *tau;
526        p1 =    tonal->music_prob *(1-tau) + (1-tonal->music_prob)*tau;
527        /* We apply the current probability with exponent beta to work around
528           the fact that the probability estimates aren't independent. */
529        p0 *= (float)pow(1-frame_probs[0], beta);
530        p1 *= (float)pow(frame_probs[0], beta);
531        /* Normalise the probabilities to get the Marokv probability of music. */
532        tonal->music_prob = p1/(p0+p1);
533        info->music_prob = tonal->music_prob;
534
535        /* This chunk of code deals with delayed decision. */
536        psum=1e-20f;
537        /* Instantaneous probability of speech and music, with beta pre-applied. */
538        speech0 = (float)pow(1-frame_probs[0], beta);
539        music0  = (float)pow(frame_probs[0], beta);
540        if (tonal->count==1)
541        {
542           tonal->pspeech[0]=.5;
543           tonal->pmusic [0]=.5;
544        }
545        /* Updated probability of having only speech (s0) or only music (m0),
546           before considering the new observation. */
547        s0 = tonal->pspeech[0] + tonal->pspeech[1];
548        m0 = tonal->pmusic [0] + tonal->pmusic [1];
549        /* Updates s0 and m0 with instantaneous probability. */
550        tonal->pspeech[0] = s0*(1-tau)*speech0;
551        tonal->pmusic [0] = m0*(1-tau)*music0;
552        /* Propagate the transition probabilities */
553        for (i=1;i<DETECT_SIZE-1;i++)
554        {
555           tonal->pspeech[i] = tonal->pspeech[i+1]*speech0;
556           tonal->pmusic [i] = tonal->pmusic [i+1]*music0;
557        }
558        /* Probability that the latest frame is speech, when all the previous ones were music. */
559        tonal->pspeech[DETECT_SIZE-1] = m0*tau*speech0;
560        /* Probability that the latest frame is music, when all the previous ones were speech. */
561        tonal->pmusic [DETECT_SIZE-1] = s0*tau*music0;
562
563        /* Renormalise probabilities to 1 */
564        for (i=0;i<DETECT_SIZE;i++)
565           psum += tonal->pspeech[i] + tonal->pmusic[i];
566        psum = 1.f/psum;
567        for (i=0;i<DETECT_SIZE;i++)
568        {
569           tonal->pspeech[i] *= psum;
570           tonal->pmusic [i] *= psum;
571        }
572        psum = tonal->pmusic[0];
573        for (i=1;i<DETECT_SIZE;i++)
574           psum += tonal->pspeech[i];
575
576        /* Estimate our confidence in the speech/music decisions */
577        if (frame_probs[1]>.75)
578        {
579           if (tonal->music_prob>.9)
580           {
581              float adapt;
582              adapt = 1.f/(++tonal->music_confidence_count);
583              tonal->music_confidence_count = IMIN(tonal->music_confidence_count, 500);
584              tonal->music_confidence += adapt*MAX16(-.2f,frame_probs[0]-tonal->music_confidence);
585           }
586           if (tonal->music_prob<.1)
587           {
588              float adapt;
589              adapt = 1.f/(++tonal->speech_confidence_count);
590              tonal->speech_confidence_count = IMIN(tonal->speech_confidence_count, 500);
591              tonal->speech_confidence += adapt*MIN16(.2f,frame_probs[0]-tonal->speech_confidence);
592           }
593        } else {
594           if (tonal->music_confidence_count==0)
595              tonal->music_confidence = .9f;
596           if (tonal->speech_confidence_count==0)
597              tonal->speech_confidence = .1f;
598        }
599     }
600     if (tonal->last_music != (tonal->music_prob>.5f))
601        tonal->last_transition=0;
602     tonal->last_music = tonal->music_prob>.5f;
603 #else
604     info->music_prob = 0;
605 #endif
606     /*for (i=0;i<25;i++)
607        printf("%f ", features[i]);
608     printf("\n");*/
609
610     info->bandwidth = bandwidth;
611     /*printf("%d %d\n", info->bandwidth, info->opus_bandwidth);*/
612     info->noisiness = frame_noisiness;
613     info->valid = 1;
614     if (info_out!=NULL)
615        OPUS_COPY(info_out, info, 1);
616     RESTORE_STACK;
617 }
618
619 void run_analysis(TonalityAnalysisState *analysis, const CELTMode *celt_mode, const void *analysis_pcm,
620                  int analysis_frame_size, int frame_size, int c1, int c2, int C, opus_int32 Fs,
621                  int lsb_depth, downmix_func downmix, AnalysisInfo *analysis_info)
622 {
623    int offset;
624    int pcm_len;
625
626    if (analysis_pcm != NULL)
627    {
628       /* Avoid overflow/wrap-around of the analysis buffer */
629       analysis_frame_size = IMIN((DETECT_SIZE-5)*Fs/100, analysis_frame_size);
630
631       pcm_len = analysis_frame_size - analysis->analysis_offset;
632       offset = analysis->analysis_offset;
633       do {
634          tonality_analysis(analysis, NULL, celt_mode, analysis_pcm, IMIN(480, pcm_len), offset, c1, c2, C, lsb_depth, downmix);
635          offset += 480;
636          pcm_len -= 480;
637       } while (pcm_len>0);
638       analysis->analysis_offset = analysis_frame_size;
639
640       analysis->analysis_offset -= frame_size;
641    }
642
643    analysis_info->valid = 0;
644    tonality_get_info(analysis, analysis_info, frame_size);
645 }