Initialize speech/music prob based on application
[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 #define ANALYSIS_C
33
34 #include <stdio.h>
35
36 #include "mathops.h"
37 #include "kiss_fft.h"
38 #include "celt.h"
39 #include "modes.h"
40 #include "arch.h"
41 #include "quant_bands.h"
42 #include "analysis.h"
43 #include "mlp.h"
44 #include "stack_alloc.h"
45 #include "float_cast.h"
46
47 #ifndef M_PI
48 #define M_PI 3.141592653
49 #endif
50
51 #ifndef DISABLE_FLOAT_API
52
53 static const float dct_table[128] = {
54         0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f,
55         0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f,
56         0.351851f, 0.338330f, 0.311806f, 0.273300f, 0.224292f, 0.166664f, 0.102631f, 0.034654f,
57        -0.034654f,-0.102631f,-0.166664f,-0.224292f,-0.273300f,-0.311806f,-0.338330f,-0.351851f,
58         0.346760f, 0.293969f, 0.196424f, 0.068975f,-0.068975f,-0.196424f,-0.293969f,-0.346760f,
59        -0.346760f,-0.293969f,-0.196424f,-0.068975f, 0.068975f, 0.196424f, 0.293969f, 0.346760f,
60         0.338330f, 0.224292f, 0.034654f,-0.166664f,-0.311806f,-0.351851f,-0.273300f,-0.102631f,
61         0.102631f, 0.273300f, 0.351851f, 0.311806f, 0.166664f,-0.034654f,-0.224292f,-0.338330f,
62         0.326641f, 0.135299f,-0.135299f,-0.326641f,-0.326641f,-0.135299f, 0.135299f, 0.326641f,
63         0.326641f, 0.135299f,-0.135299f,-0.326641f,-0.326641f,-0.135299f, 0.135299f, 0.326641f,
64         0.311806f, 0.034654f,-0.273300f,-0.338330f,-0.102631f, 0.224292f, 0.351851f, 0.166664f,
65        -0.166664f,-0.351851f,-0.224292f, 0.102631f, 0.338330f, 0.273300f,-0.034654f,-0.311806f,
66         0.293969f,-0.068975f,-0.346760f,-0.196424f, 0.196424f, 0.346760f, 0.068975f,-0.293969f,
67        -0.293969f, 0.068975f, 0.346760f, 0.196424f,-0.196424f,-0.346760f,-0.068975f, 0.293969f,
68         0.273300f,-0.166664f,-0.338330f, 0.034654f, 0.351851f, 0.102631f,-0.311806f,-0.224292f,
69         0.224292f, 0.311806f,-0.102631f,-0.351851f,-0.034654f, 0.338330f, 0.166664f,-0.273300f,
70 };
71
72 static const float analysis_window[240] = {
73       0.000043f, 0.000171f, 0.000385f, 0.000685f, 0.001071f, 0.001541f, 0.002098f, 0.002739f,
74       0.003466f, 0.004278f, 0.005174f, 0.006156f, 0.007222f, 0.008373f, 0.009607f, 0.010926f,
75       0.012329f, 0.013815f, 0.015385f, 0.017037f, 0.018772f, 0.020590f, 0.022490f, 0.024472f,
76       0.026535f, 0.028679f, 0.030904f, 0.033210f, 0.035595f, 0.038060f, 0.040604f, 0.043227f,
77       0.045928f, 0.048707f, 0.051564f, 0.054497f, 0.057506f, 0.060591f, 0.063752f, 0.066987f,
78       0.070297f, 0.073680f, 0.077136f, 0.080665f, 0.084265f, 0.087937f, 0.091679f, 0.095492f,
79       0.099373f, 0.103323f, 0.107342f, 0.111427f, 0.115579f, 0.119797f, 0.124080f, 0.128428f,
80       0.132839f, 0.137313f, 0.141849f, 0.146447f, 0.151105f, 0.155823f, 0.160600f, 0.165435f,
81       0.170327f, 0.175276f, 0.180280f, 0.185340f, 0.190453f, 0.195619f, 0.200838f, 0.206107f,
82       0.211427f, 0.216797f, 0.222215f, 0.227680f, 0.233193f, 0.238751f, 0.244353f, 0.250000f,
83       0.255689f, 0.261421f, 0.267193f, 0.273005f, 0.278856f, 0.284744f, 0.290670f, 0.296632f,
84       0.302628f, 0.308658f, 0.314721f, 0.320816f, 0.326941f, 0.333097f, 0.339280f, 0.345492f,
85       0.351729f, 0.357992f, 0.364280f, 0.370590f, 0.376923f, 0.383277f, 0.389651f, 0.396044f,
86       0.402455f, 0.408882f, 0.415325f, 0.421783f, 0.428254f, 0.434737f, 0.441231f, 0.447736f,
87       0.454249f, 0.460770f, 0.467298f, 0.473832f, 0.480370f, 0.486912f, 0.493455f, 0.500000f,
88       0.506545f, 0.513088f, 0.519630f, 0.526168f, 0.532702f, 0.539230f, 0.545751f, 0.552264f,
89       0.558769f, 0.565263f, 0.571746f, 0.578217f, 0.584675f, 0.591118f, 0.597545f, 0.603956f,
90       0.610349f, 0.616723f, 0.623077f, 0.629410f, 0.635720f, 0.642008f, 0.648271f, 0.654508f,
91       0.660720f, 0.666903f, 0.673059f, 0.679184f, 0.685279f, 0.691342f, 0.697372f, 0.703368f,
92       0.709330f, 0.715256f, 0.721144f, 0.726995f, 0.732807f, 0.738579f, 0.744311f, 0.750000f,
93       0.755647f, 0.761249f, 0.766807f, 0.772320f, 0.777785f, 0.783203f, 0.788573f, 0.793893f,
94       0.799162f, 0.804381f, 0.809547f, 0.814660f, 0.819720f, 0.824724f, 0.829673f, 0.834565f,
95       0.839400f, 0.844177f, 0.848895f, 0.853553f, 0.858151f, 0.862687f, 0.867161f, 0.871572f,
96       0.875920f, 0.880203f, 0.884421f, 0.888573f, 0.892658f, 0.896677f, 0.900627f, 0.904508f,
97       0.908321f, 0.912063f, 0.915735f, 0.919335f, 0.922864f, 0.926320f, 0.929703f, 0.933013f,
98       0.936248f, 0.939409f, 0.942494f, 0.945503f, 0.948436f, 0.951293f, 0.954072f, 0.956773f,
99       0.959396f, 0.961940f, 0.964405f, 0.966790f, 0.969096f, 0.971321f, 0.973465f, 0.975528f,
100       0.977510f, 0.979410f, 0.981228f, 0.982963f, 0.984615f, 0.986185f, 0.987671f, 0.989074f,
101       0.990393f, 0.991627f, 0.992778f, 0.993844f, 0.994826f, 0.995722f, 0.996534f, 0.997261f,
102       0.997902f, 0.998459f, 0.998929f, 0.999315f, 0.999615f, 0.999829f, 0.999957f, 1.000000f,
103 };
104
105 static const int tbands[NB_TBANDS+1] = {
106       4, 8, 12, 16, 20, 24, 28, 32, 40, 48, 56, 64, 80, 96, 112, 136, 160, 192, 240
107 };
108
109 #define NB_TONAL_SKIP_BANDS 9
110
111 static opus_val32 silk_resampler_down2_hp(
112     opus_val32                  *S,                 /* I/O  State vector [ 2 ]                                          */
113     opus_val32                  *out,               /* O    Output signal [ floor(len/2) ]                              */
114     const opus_val32            *in,                /* I    Input signal [ len ]                                        */
115     int                         inLen               /* I    Number of input samples                                     */
116 )
117 {
118     int k, len2 = inLen/2;
119     opus_val32 in32, out32, out32_hp, Y, X;
120     opus_val64 hp_ener = 0;
121     /* Internal variables and state are in Q10 format */
122     for( k = 0; k < len2; k++ ) {
123         /* Convert to Q10 */
124         in32 = in[ 2 * k ];
125
126         /* All-pass section for even input sample */
127         Y      = SUB32( in32, S[ 0 ] );
128         X      = MULT16_32_Q15(QCONST16(0.6074371f, 15), Y);
129         out32  = ADD32( S[ 0 ], X );
130         S[ 0 ] = ADD32( in32, X );
131         out32_hp = out32;
132         /* Convert to Q10 */
133         in32 = in[ 2 * k + 1 ];
134
135         /* All-pass section for odd input sample, and add to output of previous section */
136         Y      = SUB32( in32, S[ 1 ] );
137         X      = MULT16_32_Q15(QCONST16(0.15063f, 15), Y);
138         out32  = ADD32( out32, S[ 1 ] );
139         out32  = ADD32( out32, X );
140         S[ 1 ] = ADD32( in32, X );
141
142         Y      = SUB32( -in32, S[ 2 ] );
143         X      = MULT16_32_Q15(QCONST16(0.15063f, 15), Y);
144         out32_hp  = ADD32( out32_hp, S[ 2 ] );
145         out32_hp  = ADD32( out32_hp, X );
146         S[ 2 ] = ADD32( -in32, X );
147
148         hp_ener += out32_hp*(opus_val64)out32_hp;
149         /* Add, convert back to int16 and store to output */
150         out[ k ] = HALF32(out32);
151     }
152 #ifdef FIXED_POINT
153     /* len2 can be up to 480, so we shift by 8 more to make it fit. */
154     hp_ener = hp_ener >> (2*SIG_SHIFT + 8);
155 #endif
156     return (opus_val32)hp_ener;
157 }
158
159 static opus_val32 downmix_and_resample(downmix_func downmix, const void *_x, opus_val32 *y, opus_val32 S[3], int subframe, int offset, int c1, int c2, int C, int Fs)
160 {
161    VARDECL(opus_val32, tmp);
162    opus_val32 scale;
163    int j;
164    opus_val32 ret = 0;
165    SAVE_STACK;
166
167    if (subframe==0) return 0;
168    if (Fs == 48000)
169    {
170       subframe *= 2;
171       offset *= 2;
172    } else if (Fs == 16000) {
173       subframe = subframe*2/3;
174       offset = offset*2/3;
175    }
176    ALLOC(tmp, subframe, opus_val32);
177
178    downmix(_x, tmp, subframe, offset, c1, c2, C);
179 #ifdef FIXED_POINT
180    scale = (1<<SIG_SHIFT);
181 #else
182    scale = 1.f/32768;
183 #endif
184    if (c2==-2)
185       scale /= C;
186    else if (c2>-1)
187       scale /= 2;
188    for (j=0;j<subframe;j++)
189       tmp[j] *= scale;
190    if (Fs == 48000)
191    {
192       ret = silk_resampler_down2_hp(S, y, tmp, subframe);
193    } else if (Fs == 24000) {
194       OPUS_COPY(y, tmp, subframe);
195    } else if (Fs == 16000) {
196       VARDECL(opus_val32, tmp3x);
197       ALLOC(tmp3x, 3*subframe, opus_val32);
198       /* Don't do this at home! This resampler is horrible and it's only (barely)
199          usable for the purpose of the analysis because we don't care about all
200          the aliasing between 8 kHz and 12 kHz. */
201       for (j=0;j<subframe;j++)
202       {
203          tmp3x[3*j] = tmp[j];
204          tmp3x[3*j+1] = tmp[j];
205          tmp3x[3*j+2] = tmp[j];
206       }
207       silk_resampler_down2_hp(S, y, tmp3x, 3*subframe);
208    }
209    RESTORE_STACK;
210    return ret;
211 }
212
213 void tonality_analysis_init(TonalityAnalysisState *tonal, opus_int32 Fs)
214 {
215   /* Initialize reusable fields. */
216   tonal->arch = opus_select_arch();
217   tonal->Fs = Fs;
218   /* Clear remaining fields. */
219   tonality_analysis_reset(tonal);
220 }
221
222 void tonality_analysis_reset(TonalityAnalysisState *tonal)
223 {
224   /* Clear non-reusable fields. */
225   char *start = (char*)&tonal->TONALITY_ANALYSIS_RESET_START;
226   OPUS_CLEAR(start, sizeof(TonalityAnalysisState) - (start - (char*)tonal));
227   tonal->music_confidence = .9f;
228   tonal->speech_confidence = .1f;
229 }
230
231 void tonality_get_info(TonalityAnalysisState *tonal, AnalysisInfo *info_out, int len)
232 {
233    int pos;
234    int curr_lookahead;
235    float psum;
236    int i;
237
238    pos = tonal->read_pos;
239    curr_lookahead = tonal->write_pos-tonal->read_pos;
240    if (curr_lookahead<0)
241       curr_lookahead += DETECT_SIZE;
242
243    /* On long frames, look at the second analysis window rather than the first. */
244    if (len > tonal->Fs/50 && pos != tonal->write_pos)
245    {
246       pos++;
247       if (pos==DETECT_SIZE)
248          pos=0;
249    }
250    if (pos == tonal->write_pos)
251       pos--;
252    if (pos<0)
253       pos = DETECT_SIZE-1;
254    OPUS_COPY(info_out, &tonal->info[pos], 1);
255    /* If possible, look ahead for a tone to compensate for the delay in the tone detector. */
256    for (i=0;i<3;i++)
257    {
258       pos++;
259       if (pos==DETECT_SIZE)
260          pos = 0;
261       if (pos == tonal->write_pos)
262          break;
263       info_out->tonality = MAX32(0, -.03f + MAX32(info_out->tonality, tonal->info[pos].tonality-.05f));
264    }
265    tonal->read_subframe += len/(tonal->Fs/400);
266    while (tonal->read_subframe>=8)
267    {
268       tonal->read_subframe -= 8;
269       tonal->read_pos++;
270    }
271    if (tonal->read_pos>=DETECT_SIZE)
272       tonal->read_pos-=DETECT_SIZE;
273
274    /* The -1 is to compensate for the delay in the features themselves. */
275    curr_lookahead = IMAX(curr_lookahead-1, 0);
276
277    psum=0;
278    /* Summing the probability of transition patterns that involve music at
279       time (DETECT_SIZE-curr_lookahead-1) */
280    for (i=0;i<DETECT_SIZE-curr_lookahead;i++)
281       psum += tonal->pmusic[i];
282    for (;i<DETECT_SIZE;i++)
283       psum += tonal->pspeech[i];
284    psum = psum*tonal->music_confidence + (1-psum)*tonal->speech_confidence;
285    /*printf("%f %f %f %f %f\n", psum, info_out->music_prob, info_out->vad_prob, info_out->activity_probability, info_out->tonality);*/
286
287    info_out->music_prob = psum;
288 }
289
290 static const float std_feature_bias[9] = {
291       5.684947f, 3.475288f, 1.770634f, 1.599784f, 3.773215f,
292       2.163313f, 1.260756f, 1.116868f, 1.918795f
293 };
294
295 #define LEAKAGE_OFFSET 2.5f
296 #define LEAKAGE_SLOPE 2.f
297
298 #ifdef FIXED_POINT
299 /* For fixed-point, the input is +/-2^15 shifted up by SIG_SHIFT, so we need to
300    compensate for that in the energy. */
301 #define SCALE_COMPENS (1.f/((opus_int32)1<<(15+SIG_SHIFT)))
302 #define SCALE_ENER(e) ((SCALE_COMPENS*SCALE_COMPENS)*(e))
303 #else
304 #define SCALE_ENER(e) (e)
305 #endif
306
307 static void tonality_analysis(TonalityAnalysisState *tonal, const CELTMode *celt_mode, const void *x, int len, int offset, int c1, int c2, int C, int lsb_depth, downmix_func downmix)
308 {
309     int i, b;
310     const kiss_fft_state *kfft;
311     VARDECL(kiss_fft_cpx, in);
312     VARDECL(kiss_fft_cpx, out);
313     int N = 480, N2=240;
314     float * OPUS_RESTRICT A = tonal->angle;
315     float * OPUS_RESTRICT dA = tonal->d_angle;
316     float * OPUS_RESTRICT d2A = tonal->d2_angle;
317     VARDECL(float, tonality);
318     VARDECL(float, noisiness);
319     float band_tonality[NB_TBANDS];
320     float logE[NB_TBANDS];
321     float BFCC[8];
322     float features[25];
323     float frame_tonality;
324     float max_frame_tonality;
325     /*float tw_sum=0;*/
326     float frame_noisiness;
327     const float pi4 = (float)(M_PI*M_PI*M_PI*M_PI);
328     float slope=0;
329     float frame_stationarity;
330     float relativeE;
331     float frame_probs[2];
332     float alpha, alphaE, alphaE2;
333     float frame_loudness;
334     float bandwidth_mask;
335     int bandwidth=0;
336     float maxE = 0;
337     float noise_floor;
338     int remaining;
339     AnalysisInfo *info;
340     float hp_ener;
341     float tonality2[240];
342     float midE[8];
343     float spec_variability=0;
344     float band_log2[NB_TBANDS+1];
345     float leakage_from[NB_TBANDS+1];
346     float leakage_to[NB_TBANDS+1];
347     SAVE_STACK;
348
349     alpha = 1.f/IMIN(10, 1+tonal->count);
350     alphaE = 1.f/IMIN(25, 1+tonal->count);
351     alphaE2 = 1.f/IMIN(500, 1+tonal->count);
352
353     if (tonal->Fs == 48000)
354     {
355        /* len and offset are now at 24 kHz. */
356        len/= 2;
357        offset /= 2;
358     } else if (tonal->Fs == 16000) {
359        len = 3*len/2;
360        offset = 3*offset/2;
361     }
362
363     if (tonal->count<4) {
364        if (tonal->application == OPUS_APPLICATION_VOIP)
365           tonal->music_prob = .1;
366        else
367           tonal->music_prob = .625;
368     }
369     kfft = celt_mode->mdct.kfft[0];
370     if (tonal->count==0)
371        tonal->mem_fill = 240;
372     tonal->hp_ener_accum += (float)downmix_and_resample(downmix, x,
373           &tonal->inmem[tonal->mem_fill], tonal->downmix_state,
374           IMIN(len, ANALYSIS_BUF_SIZE-tonal->mem_fill), offset, c1, c2, C, tonal->Fs);
375     if (tonal->mem_fill+len < ANALYSIS_BUF_SIZE)
376     {
377        tonal->mem_fill += len;
378        /* Don't have enough to update the analysis */
379        RESTORE_STACK;
380        return;
381     }
382     hp_ener = tonal->hp_ener_accum;
383     info = &tonal->info[tonal->write_pos++];
384     if (tonal->write_pos>=DETECT_SIZE)
385        tonal->write_pos-=DETECT_SIZE;
386
387     ALLOC(in, 480, kiss_fft_cpx);
388     ALLOC(out, 480, kiss_fft_cpx);
389     ALLOC(tonality, 240, float);
390     ALLOC(noisiness, 240, float);
391     for (i=0;i<N2;i++)
392     {
393        float w = analysis_window[i];
394        in[i].r = (kiss_fft_scalar)(w*tonal->inmem[i]);
395        in[i].i = (kiss_fft_scalar)(w*tonal->inmem[N2+i]);
396        in[N-i-1].r = (kiss_fft_scalar)(w*tonal->inmem[N-i-1]);
397        in[N-i-1].i = (kiss_fft_scalar)(w*tonal->inmem[N+N2-i-1]);
398     }
399     OPUS_MOVE(tonal->inmem, tonal->inmem+ANALYSIS_BUF_SIZE-240, 240);
400     remaining = len - (ANALYSIS_BUF_SIZE-tonal->mem_fill);
401     tonal->hp_ener_accum = (float)downmix_and_resample(downmix, x,
402           &tonal->inmem[240], tonal->downmix_state, remaining,
403           offset+ANALYSIS_BUF_SIZE-tonal->mem_fill, c1, c2, C, tonal->Fs);
404     tonal->mem_fill = 240 + remaining;
405     opus_fft(kfft, in, out, tonal->arch);
406 #ifndef FIXED_POINT
407     /* If there's any NaN on the input, the entire output will be NaN, so we only need to check one value. */
408     if (celt_isnan(out[0].r))
409     {
410        info->valid = 0;
411        RESTORE_STACK;
412        return;
413     }
414 #endif
415
416     for (i=1;i<N2;i++)
417     {
418        float X1r, X2r, X1i, X2i;
419        float angle, d_angle, d2_angle;
420        float angle2, d_angle2, d2_angle2;
421        float mod1, mod2, avg_mod;
422        X1r = (float)out[i].r+out[N-i].r;
423        X1i = (float)out[i].i-out[N-i].i;
424        X2r = (float)out[i].i+out[N-i].i;
425        X2i = (float)out[N-i].r-out[i].r;
426
427        angle = (float)(.5f/M_PI)*fast_atan2f(X1i, X1r);
428        d_angle = angle - A[i];
429        d2_angle = d_angle - dA[i];
430
431        angle2 = (float)(.5f/M_PI)*fast_atan2f(X2i, X2r);
432        d_angle2 = angle2 - angle;
433        d2_angle2 = d_angle2 - d_angle;
434
435        mod1 = d2_angle - (float)float2int(d2_angle);
436        noisiness[i] = ABS16(mod1);
437        mod1 *= mod1;
438        mod1 *= mod1;
439
440        mod2 = d2_angle2 - (float)float2int(d2_angle2);
441        noisiness[i] += ABS16(mod2);
442        mod2 *= mod2;
443        mod2 *= mod2;
444
445        avg_mod = .25f*(d2A[i]+mod1+2*mod2);
446        /* This introduces an extra delay of 2 frames in the detection. */
447        tonality[i] = 1.f/(1.f+40.f*16.f*pi4*avg_mod)-.015f;
448        /* No delay on this detection, but it's less reliable. */
449        tonality2[i] = 1.f/(1.f+40.f*16.f*pi4*mod2)-.015f;
450
451        A[i] = angle2;
452        dA[i] = d_angle2;
453        d2A[i] = mod2;
454     }
455     for (i=2;i<N2-1;i++)
456     {
457        float tt = MIN32(tonality2[i], MAX32(tonality2[i-1], tonality2[i+1]));
458        tonality[i] = .9f*MAX32(tonality[i], tt-.1f);
459     }
460     frame_tonality = 0;
461     max_frame_tonality = 0;
462     /*tw_sum = 0;*/
463     info->activity = 0;
464     frame_noisiness = 0;
465     frame_stationarity = 0;
466     if (!tonal->count)
467     {
468        for (b=0;b<NB_TBANDS;b++)
469        {
470           tonal->lowE[b] = 1e10;
471           tonal->highE[b] = -1e10;
472        }
473     }
474     relativeE = 0;
475     frame_loudness = 0;
476     /* The energy of the very first band is special because of DC. */
477     {
478        float E = 0;
479        float X1r, X2r;
480        X1r = 2*(float)out[0].r;
481        X2r = 2*(float)out[0].i;
482        E = X1r*X1r + X2r*X2r;
483        for (i=1;i<4;i++)
484        {
485           float binE = out[i].r*(float)out[i].r + out[N-i].r*(float)out[N-i].r
486                      + out[i].i*(float)out[i].i + out[N-i].i*(float)out[N-i].i;
487           E += binE;
488        }
489        E = SCALE_ENER(E);
490        band_log2[0] = .5f*1.442695f*(float)log(E+1e-10f);
491     }
492     for (b=0;b<NB_TBANDS;b++)
493     {
494        float E=0, tE=0, nE=0;
495        float L1, L2;
496        float stationarity;
497        for (i=tbands[b];i<tbands[b+1];i++)
498        {
499           float binE = out[i].r*(float)out[i].r + out[N-i].r*(float)out[N-i].r
500                      + out[i].i*(float)out[i].i + out[N-i].i*(float)out[N-i].i;
501           binE = SCALE_ENER(binE);
502           E += binE;
503           tE += binE*MAX32(0, tonality[i]);
504           nE += binE*2.f*(.5f-noisiness[i]);
505        }
506 #ifndef FIXED_POINT
507        /* Check for extreme band energies that could cause NaNs later. */
508        if (!(E<1e9f) || celt_isnan(E))
509        {
510           info->valid = 0;
511           RESTORE_STACK;
512           return;
513        }
514 #endif
515
516        tonal->E[tonal->E_count][b] = E;
517        frame_noisiness += nE/(1e-15f+E);
518
519        frame_loudness += (float)sqrt(E+1e-10f);
520        logE[b] = (float)log(E+1e-10f);
521        band_log2[b+1] = .5f*1.442695f*(float)log(E+1e-10f);
522        tonal->logE[tonal->E_count][b] = logE[b];
523        if (tonal->count==0)
524           tonal->highE[b] = tonal->lowE[b] = logE[b];
525        if (tonal->highE[b] > tonal->lowE[b] + 7.5)
526        {
527           if (tonal->highE[b] - logE[b] > logE[b] - tonal->lowE[b])
528              tonal->highE[b] -= .01f;
529           else
530              tonal->lowE[b] += .01f;
531        }
532        if (logE[b] > tonal->highE[b])
533        {
534           tonal->highE[b] = logE[b];
535           tonal->lowE[b] = MAX32(tonal->highE[b]-15, tonal->lowE[b]);
536        } else if (logE[b] < tonal->lowE[b])
537        {
538           tonal->lowE[b] = logE[b];
539           tonal->highE[b] = MIN32(tonal->lowE[b]+15, tonal->highE[b]);
540        }
541        relativeE += (logE[b]-tonal->lowE[b])/(1e-15f + (tonal->highE[b]-tonal->lowE[b]));
542
543        L1=L2=0;
544        for (i=0;i<NB_FRAMES;i++)
545        {
546           L1 += (float)sqrt(tonal->E[i][b]);
547           L2 += tonal->E[i][b];
548        }
549
550        stationarity = MIN16(0.99f,L1/(float)sqrt(1e-15+NB_FRAMES*L2));
551        stationarity *= stationarity;
552        stationarity *= stationarity;
553        frame_stationarity += stationarity;
554        /*band_tonality[b] = tE/(1e-15+E)*/;
555        band_tonality[b] = MAX16(tE/(1e-15f+E), stationarity*tonal->prev_band_tonality[b]);
556 #if 0
557        if (b>=NB_TONAL_SKIP_BANDS)
558        {
559           frame_tonality += tweight[b]*band_tonality[b];
560           tw_sum += tweight[b];
561        }
562 #else
563        frame_tonality += band_tonality[b];
564        if (b>=NB_TBANDS-NB_TONAL_SKIP_BANDS)
565           frame_tonality -= band_tonality[b-NB_TBANDS+NB_TONAL_SKIP_BANDS];
566 #endif
567        max_frame_tonality = MAX16(max_frame_tonality, (1.f+.03f*(b-NB_TBANDS))*frame_tonality);
568        slope += band_tonality[b]*(b-8);
569        /*printf("%f %f ", band_tonality[b], stationarity);*/
570        tonal->prev_band_tonality[b] = band_tonality[b];
571     }
572
573     leakage_from[0] = band_log2[0];
574     leakage_to[0] = band_log2[0] - LEAKAGE_OFFSET;
575     for (b=1;b<NB_TBANDS+1;b++)
576     {
577        float leak_slope = LEAKAGE_SLOPE*(tbands[b]-tbands[b-1])/4;
578        leakage_from[b] = MIN16(leakage_from[b-1]+leak_slope, band_log2[b]);
579        leakage_to[b] = MAX16(leakage_to[b-1]-leak_slope, band_log2[b]-LEAKAGE_OFFSET);
580     }
581     for (b=NB_TBANDS-2;b>=0;b--)
582     {
583        float leak_slope = LEAKAGE_SLOPE*(tbands[b+1]-tbands[b])/4;
584        leakage_from[b] = MIN16(leakage_from[b+1]+leak_slope, leakage_from[b]);
585        leakage_to[b] = MAX16(leakage_to[b+1]-leak_slope, leakage_to[b]);
586     }
587     celt_assert(NB_TBANDS+1 <= LEAK_BANDS);
588     for (b=0;b<NB_TBANDS+1;b++)
589     {
590        /* leak_boost[] is made up of two terms. The first, based on leakage_to[],
591           represents the boost needed to overcome the amount of analysis leakage
592           cause in a weaker band b by louder neighbouring bands.
593           The second, based on leakage_from[], applies to a loud band b for
594           which the quantization noise causes synthesis leakage to the weaker
595           neighbouring bands. */
596        float boost = MAX16(0, leakage_to[b] - band_log2[b]) +
597              MAX16(0, band_log2[b] - (leakage_from[b]+LEAKAGE_OFFSET));
598        info->leak_boost[b] = IMIN(255, (int)floor(.5 + 64.f*boost));
599     }
600     for (;b<LEAK_BANDS;b++) info->leak_boost[b] = 0;
601
602     for (i=0;i<NB_FRAMES;i++)
603     {
604        int j;
605        float mindist = 1e15f;
606        for (j=0;j<NB_FRAMES;j++)
607        {
608           int k;
609           float dist=0;
610           for (k=0;k<NB_TBANDS;k++)
611           {
612              float tmp;
613              tmp = tonal->logE[i][k] - tonal->logE[j][k];
614              dist += tmp*tmp;
615           }
616           if (j!=i)
617              mindist = MIN32(mindist, dist);
618        }
619        spec_variability += mindist;
620     }
621     spec_variability = (float)sqrt(spec_variability/NB_FRAMES/NB_TBANDS);
622     bandwidth_mask = 0;
623     bandwidth = 0;
624     maxE = 0;
625     noise_floor = 5.7e-4f/(1<<(IMAX(0,lsb_depth-8)));
626     noise_floor *= noise_floor;
627     for (b=0;b<NB_TBANDS;b++)
628     {
629        float E=0;
630        int band_start, band_end;
631        /* Keep a margin of 300 Hz for aliasing */
632        band_start = tbands[b];
633        band_end = tbands[b+1];
634        for (i=band_start;i<band_end;i++)
635        {
636           float binE = out[i].r*(float)out[i].r + out[N-i].r*(float)out[N-i].r
637                      + out[i].i*(float)out[i].i + out[N-i].i*(float)out[N-i].i;
638           E += binE;
639        }
640        E = SCALE_ENER(E);
641        maxE = MAX32(maxE, E);
642        tonal->meanE[b] = MAX32((1-alphaE2)*tonal->meanE[b], E);
643        E = MAX32(E, tonal->meanE[b]);
644        /* Use a simple follower with 13 dB/Bark slope for spreading function */
645        bandwidth_mask = MAX32(.05f*bandwidth_mask, E);
646        /* Consider the band "active" only if all these conditions are met:
647           1) less than 10 dB below the simple follower
648           2) less than 90 dB below the peak band (maximal masking possible considering
649              both the ATH and the loudness-dependent slope of the spreading function)
650           3) above the PCM quantization noise floor
651           We use b+1 because the first CELT band isn't included in tbands[]
652        */
653        if (E>.1*bandwidth_mask && E*1e9f > maxE && E > noise_floor*(band_end-band_start))
654           bandwidth = b+1;
655     }
656     /* Special case for the last two bands, for which we don't have spectrum but only
657        the energy above 12 kHz. */
658     {
659        float E = hp_ener*(1.f/(240*240));
660 #ifdef FIXED_POINT
661        /* silk_resampler_down2_hp() shifted right by an extra 8 bits. */
662        E *= 256.f*(1.f/Q15ONE)*(1.f/Q15ONE);
663 #endif
664        maxE = MAX32(maxE, E);
665        tonal->meanE[b] = MAX32((1-alphaE2)*tonal->meanE[b], E);
666        E = MAX32(E, tonal->meanE[b]);
667        /* Use a simple follower with 13 dB/Bark slope for spreading function */
668        bandwidth_mask = MAX32(.05f*bandwidth_mask, E);
669        if (E>.1*bandwidth_mask && E*1e9f > maxE && E > noise_floor*160)
670           bandwidth = 20;
671     }
672     if (tonal->count<=2)
673        bandwidth = 20;
674     frame_loudness = 20*(float)log10(frame_loudness);
675     tonal->Etracker = MAX32(tonal->Etracker-.003f, frame_loudness);
676     tonal->lowECount *= (1-alphaE);
677     if (frame_loudness < tonal->Etracker-30)
678        tonal->lowECount += alphaE;
679
680     for (i=0;i<8;i++)
681     {
682        float sum=0;
683        for (b=0;b<16;b++)
684           sum += dct_table[i*16+b]*logE[b];
685        BFCC[i] = sum;
686     }
687     for (i=0;i<8;i++)
688     {
689        float sum=0;
690        for (b=0;b<16;b++)
691           sum += dct_table[i*16+b]*.5f*(tonal->highE[b]+tonal->lowE[b]);
692        midE[i] = sum;
693     }
694
695     frame_stationarity /= NB_TBANDS;
696     relativeE /= NB_TBANDS;
697     if (tonal->count<10)
698        relativeE = .5;
699     frame_noisiness /= NB_TBANDS;
700 #if 1
701     info->activity = frame_noisiness + (1-frame_noisiness)*relativeE;
702 #else
703     info->activity = .5*(1+frame_noisiness-frame_stationarity);
704 #endif
705     frame_tonality = (max_frame_tonality/(NB_TBANDS-NB_TONAL_SKIP_BANDS));
706     frame_tonality = MAX16(frame_tonality, tonal->prev_tonality*.8f);
707     tonal->prev_tonality = frame_tonality;
708
709     slope /= 8*8;
710     info->tonality_slope = slope;
711
712     tonal->E_count = (tonal->E_count+1)%NB_FRAMES;
713     tonal->count = IMIN(tonal->count+1, ANALYSIS_COUNT_MAX);
714     info->tonality = frame_tonality;
715
716     for (i=0;i<4;i++)
717        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];
718
719     for (i=0;i<4;i++)
720        tonal->cmean[i] = (1-alpha)*tonal->cmean[i] + alpha*BFCC[i];
721
722     for (i=0;i<4;i++)
723         features[4+i] = 0.63246f*(BFCC[i]-tonal->mem[i+24]) + 0.31623f*(tonal->mem[i]-tonal->mem[i+16]);
724     for (i=0;i<3;i++)
725         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];
726
727     if (tonal->count > 5)
728     {
729        for (i=0;i<9;i++)
730           tonal->std[i] = (1-alpha)*tonal->std[i] + alpha*features[i]*features[i];
731     }
732     for (i=0;i<4;i++)
733        features[i] = BFCC[i]-midE[i];
734
735     for (i=0;i<8;i++)
736     {
737        tonal->mem[i+24] = tonal->mem[i+16];
738        tonal->mem[i+16] = tonal->mem[i+8];
739        tonal->mem[i+8] = tonal->mem[i];
740        tonal->mem[i] = BFCC[i];
741     }
742     for (i=0;i<9;i++)
743        features[11+i] = (float)sqrt(tonal->std[i]) - std_feature_bias[i];
744     features[18] = spec_variability - 0.78f;
745     features[20] = info->tonality - 0.154723f;
746     features[21] = info->activity - 0.724643f;
747     features[22] = frame_stationarity - 0.743717f;
748     features[23] = info->tonality_slope + 0.069216f;
749     features[24] = tonal->lowECount - 0.067930f;
750
751     mlp_process(&net, features, frame_probs);
752     frame_probs[0] = .5f*(frame_probs[0]+1);
753     /* Curve fitting between the MLP probability and the actual probability */
754     /*frame_probs[0] = .01f + 1.21f*frame_probs[0]*frame_probs[0] - .23f*(float)pow(frame_probs[0], 10);*/
755     /* Probability of active audio (as opposed to silence) */
756     frame_probs[1] = .5f*frame_probs[1]+.5f;
757     frame_probs[1] *= frame_probs[1];
758
759     /* Probability of speech or music vs noise */
760     info->activity_probability = frame_probs[1];
761
762     /*printf("%f %f\n", frame_probs[0], frame_probs[1]);*/
763     {
764        /* Probability of state transition */
765        float tau;
766        /* Represents independence of the MLP probabilities, where
767           beta=1 means fully independent. */
768        float beta;
769        /* Denormalized probability of speech (p0) and music (p1) after update */
770        float p0, p1;
771        /* Probabilities for "all speech" and "all music" */
772        float s0, m0;
773        /* Probability sum for renormalisation */
774        float psum;
775        /* Instantaneous probability of speech and music, with beta pre-applied. */
776        float speech0;
777        float music0;
778        float p, q;
779
780        /* More silence transitions for speech than for music. */
781        tau = .001f*tonal->music_prob + .01f*(1-tonal->music_prob);
782        p = MAX16(.05f,MIN16(.95f,frame_probs[1]));
783        q = MAX16(.05f,MIN16(.95f,tonal->vad_prob));
784        beta = .02f+.05f*ABS16(p-q)/(p*(1-q)+q*(1-p));
785        /* p0 and p1 are the probabilities of speech and music at this frame
786           using only information from previous frame and applying the
787           state transition model */
788        p0 = (1-tonal->vad_prob)*(1-tau) +    tonal->vad_prob *tau;
789        p1 =    tonal->vad_prob *(1-tau) + (1-tonal->vad_prob)*tau;
790        /* We apply the current probability with exponent beta to work around
791           the fact that the probability estimates aren't independent. */
792        p0 *= (float)pow(1-frame_probs[1], beta);
793        p1 *= (float)pow(frame_probs[1], beta);
794        /* Normalise the probabilities to get the Marokv probability of music. */
795        tonal->vad_prob = p1/(p0+p1);
796        info->vad_prob = tonal->vad_prob;
797        /* Consider that silence has a 50-50 probability of being speech or music. */
798        frame_probs[0] = tonal->vad_prob*frame_probs[0] + (1-tonal->vad_prob)*.5f;
799
800        /* One transition every 3 minutes of active audio */
801        tau = .0001f;
802        /* Adapt beta based on how "unexpected" the new prob is */
803        p = MAX16(.05f,MIN16(.95f,frame_probs[0]));
804        q = MAX16(.05f,MIN16(.95f,tonal->music_prob));
805        beta = .02f+.05f*ABS16(p-q)/(p*(1-q)+q*(1-p));
806        /* p0 and p1 are the probabilities of speech and music at this frame
807           using only information from previous frame and applying the
808           state transition model */
809        p0 = (1-tonal->music_prob)*(1-tau) +    tonal->music_prob *tau;
810        p1 =    tonal->music_prob *(1-tau) + (1-tonal->music_prob)*tau;
811        /* We apply the current probability with exponent beta to work around
812           the fact that the probability estimates aren't independent. */
813        p0 *= (float)pow(1-frame_probs[0], beta);
814        p1 *= (float)pow(frame_probs[0], beta);
815        /* Normalise the probabilities to get the Marokv probability of music. */
816        tonal->music_prob = p1/(p0+p1);
817        info->music_prob = tonal->music_prob;
818
819        /*printf("%f %f %f %f\n", frame_probs[0], frame_probs[1], tonal->music_prob, tonal->vad_prob);*/
820        /* This chunk of code deals with delayed decision. */
821        psum=1e-20f;
822        /* Instantaneous probability of speech and music, with beta pre-applied. */
823        speech0 = (float)pow(1-frame_probs[0], beta);
824        music0  = (float)pow(frame_probs[0], beta);
825        if (tonal->count==1)
826        {
827           if (tonal->application == OPUS_APPLICATION_VOIP)
828              tonal->pmusic[0] = .1;
829           else
830              tonal->pmusic[0] = .625;
831           tonal->pspeech[0] = 1-tonal->pmusic[0];
832        }
833        /* Updated probability of having only speech (s0) or only music (m0),
834           before considering the new observation. */
835        s0 = tonal->pspeech[0] + tonal->pspeech[1];
836        m0 = tonal->pmusic [0] + tonal->pmusic [1];
837        /* Updates s0 and m0 with instantaneous probability. */
838        tonal->pspeech[0] = s0*(1-tau)*speech0;
839        tonal->pmusic [0] = m0*(1-tau)*music0;
840        /* Propagate the transition probabilities */
841        for (i=1;i<DETECT_SIZE-1;i++)
842        {
843           tonal->pspeech[i] = tonal->pspeech[i+1]*speech0;
844           tonal->pmusic [i] = tonal->pmusic [i+1]*music0;
845        }
846        /* Probability that the latest frame is speech, when all the previous ones were music. */
847        tonal->pspeech[DETECT_SIZE-1] = m0*tau*speech0;
848        /* Probability that the latest frame is music, when all the previous ones were speech. */
849        tonal->pmusic [DETECT_SIZE-1] = s0*tau*music0;
850
851        /* Renormalise probabilities to 1 */
852        for (i=0;i<DETECT_SIZE;i++)
853           psum += tonal->pspeech[i] + tonal->pmusic[i];
854        psum = 1.f/psum;
855        for (i=0;i<DETECT_SIZE;i++)
856        {
857           tonal->pspeech[i] *= psum;
858           tonal->pmusic [i] *= psum;
859        }
860        psum = tonal->pmusic[0];
861        for (i=1;i<DETECT_SIZE;i++)
862           psum += tonal->pspeech[i];
863
864        /* Estimate our confidence in the speech/music decisions */
865        if (frame_probs[1]>.75)
866        {
867           if (tonal->music_prob>.9)
868           {
869              float adapt;
870              adapt = 1.f/(++tonal->music_confidence_count);
871              tonal->music_confidence_count = IMIN(tonal->music_confidence_count, 500);
872              tonal->music_confidence += adapt*MAX16(-.2f,frame_probs[0]-tonal->music_confidence);
873           }
874           if (tonal->music_prob<.1)
875           {
876              float adapt;
877              adapt = 1.f/(++tonal->speech_confidence_count);
878              tonal->speech_confidence_count = IMIN(tonal->speech_confidence_count, 500);
879              tonal->speech_confidence += adapt*MIN16(.2f,frame_probs[0]-tonal->speech_confidence);
880           }
881        }
882     }
883     tonal->last_music = tonal->music_prob>.5f;
884 #ifdef MLP_TRAINING
885     for (i=0;i<25;i++)
886        printf("%f ", features[i]);
887     printf("\n");
888 #endif
889
890     info->bandwidth = bandwidth;
891     /*printf("%d %d\n", info->bandwidth, info->opus_bandwidth);*/
892     info->noisiness = frame_noisiness;
893     info->valid = 1;
894     RESTORE_STACK;
895 }
896
897 void run_analysis(TonalityAnalysisState *analysis, const CELTMode *celt_mode, const void *analysis_pcm,
898                  int analysis_frame_size, int frame_size, int c1, int c2, int C, opus_int32 Fs,
899                  int lsb_depth, downmix_func downmix, AnalysisInfo *analysis_info)
900 {
901    int offset;
902    int pcm_len;
903
904    analysis_frame_size -= analysis_frame_size&1;
905    if (analysis_pcm != NULL)
906    {
907       /* Avoid overflow/wrap-around of the analysis buffer */
908       analysis_frame_size = IMIN((DETECT_SIZE-5)*Fs/50, analysis_frame_size);
909
910       pcm_len = analysis_frame_size - analysis->analysis_offset;
911       offset = analysis->analysis_offset;
912       while (pcm_len>0) {
913          tonality_analysis(analysis, celt_mode, analysis_pcm, IMIN(Fs/50, pcm_len), offset, c1, c2, C, lsb_depth, downmix);
914          offset += Fs/50;
915          pcm_len -= Fs/50;
916       }
917       analysis->analysis_offset = analysis_frame_size;
918
919       analysis->analysis_offset -= frame_size;
920    }
921
922    analysis_info->valid = 0;
923    tonality_get_info(analysis, analysis_info, frame_size);
924 }
925
926 #endif /* DISABLE_FLOAT_API */