Adding leakage modelling to boost bands
[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 }
228
229 void tonality_get_info(TonalityAnalysisState *tonal, AnalysisInfo *info_out, int len)
230 {
231    int pos;
232    int curr_lookahead;
233    float psum;
234    int i;
235
236    pos = tonal->read_pos;
237    curr_lookahead = tonal->write_pos-tonal->read_pos;
238    if (curr_lookahead<0)
239       curr_lookahead += DETECT_SIZE;
240
241    /* On long frames, look at the second analysis window rather than the first. */
242    if (len > tonal->Fs/50 && pos != tonal->write_pos)
243    {
244       pos++;
245       if (pos==DETECT_SIZE)
246          pos=0;
247    }
248    if (pos == tonal->write_pos)
249       pos--;
250    if (pos<0)
251       pos = DETECT_SIZE-1;
252    OPUS_COPY(info_out, &tonal->info[pos], 1);
253    /* If possible, look ahead for a tone to compensate for the delay in the tone detector. */
254    for (i=0;i<3;i++)
255    {
256       pos++;
257       if (pos==DETECT_SIZE)
258          pos = 0;
259       if (pos == tonal->write_pos)
260          break;
261       info_out->tonality = MAX32(0, -.03f + MAX32(info_out->tonality, tonal->info[pos].tonality-.05f));
262    }
263    tonal->read_subframe += len/(tonal->Fs/400);
264    while (tonal->read_subframe>=8)
265    {
266       tonal->read_subframe -= 8;
267       tonal->read_pos++;
268    }
269    if (tonal->read_pos>=DETECT_SIZE)
270       tonal->read_pos-=DETECT_SIZE;
271
272    /* The -1 is to compensate for the delay in the features themselves. */
273    curr_lookahead = IMAX(curr_lookahead-1, 0);
274
275    psum=0;
276    /* Summing the probability of transition patterns that involve music at
277       time (DETECT_SIZE-curr_lookahead-1) */
278    for (i=0;i<DETECT_SIZE-curr_lookahead;i++)
279       psum += tonal->pmusic[i];
280    for (;i<DETECT_SIZE;i++)
281       psum += tonal->pspeech[i];
282    psum = psum*tonal->music_confidence + (1-psum)*tonal->speech_confidence;
283    /*printf("%f %f %f %f %f\n", psum, info_out->music_prob, info_out->vad_prob, info_out->activity_probability, info_out->tonality);*/
284
285    info_out->music_prob = psum;
286 }
287
288 static const float std_feature_bias[9] = {
289       5.684947f, 3.475288f, 1.770634f, 1.599784f, 3.773215f,
290       2.163313f, 1.260756f, 1.116868f, 1.918795f
291 };
292
293 #define LEAKAGE_OFFSET 2.5f
294 #define LEAKAGE_SLOPE 2.f
295
296 #ifdef FIXED_POINT
297 /* For fixed-point, the input is +/-2^15 shifted up by SIG_SHIFT, so we need to
298    compensate for that in the energy. */
299 #define SCALE_COMPENS (1.f/((opus_int32)1<<(15+SIG_SHIFT)))
300 #define SCALE_ENER(e) ((SCALE_COMPENS*SCALE_COMPENS)*(e))
301 #else
302 #define SCALE_ENER(e) (e)
303 #endif
304
305 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)
306 {
307     int i, b;
308     const kiss_fft_state *kfft;
309     VARDECL(kiss_fft_cpx, in);
310     VARDECL(kiss_fft_cpx, out);
311     int N = 480, N2=240;
312     float * OPUS_RESTRICT A = tonal->angle;
313     float * OPUS_RESTRICT dA = tonal->d_angle;
314     float * OPUS_RESTRICT d2A = tonal->d2_angle;
315     VARDECL(float, tonality);
316     VARDECL(float, noisiness);
317     float band_tonality[NB_TBANDS];
318     float logE[NB_TBANDS];
319     float BFCC[8];
320     float features[25];
321     float frame_tonality;
322     float max_frame_tonality;
323     /*float tw_sum=0;*/
324     float frame_noisiness;
325     const float pi4 = (float)(M_PI*M_PI*M_PI*M_PI);
326     float slope=0;
327     float frame_stationarity;
328     float relativeE;
329     float frame_probs[2];
330     float alpha, alphaE, alphaE2;
331     float frame_loudness;
332     float bandwidth_mask;
333     int bandwidth=0;
334     float maxE = 0;
335     float noise_floor;
336     int remaining;
337     AnalysisInfo *info;
338     float hp_ener;
339     float tonality2[240];
340     float midE[8];
341     float spec_variability=0;
342     float band_log2[NB_TBANDS+1];
343     float leakage_from[NB_TBANDS+1];
344     float leakage_to[NB_TBANDS+1];
345     SAVE_STACK;
346
347     alpha = 1.f/IMIN(10, 1+tonal->count);
348     alphaE = 1.f/IMIN(25, 1+tonal->count);
349     alphaE2 = 1.f/IMIN(500, 1+tonal->count);
350
351     if (tonal->Fs == 48000)
352     {
353        /* len and offset are now at 24 kHz. */
354        len/= 2;
355        offset /= 2;
356     } else if (tonal->Fs == 16000) {
357        len = 3*len/2;
358        offset = 3*offset/2;
359     }
360
361     if (tonal->count<4)
362        tonal->music_prob = .5;
363     kfft = celt_mode->mdct.kfft[0];
364     if (tonal->count==0)
365        tonal->mem_fill = 240;
366     tonal->hp_ener_accum += (float)downmix_and_resample(downmix, x,
367           &tonal->inmem[tonal->mem_fill], tonal->downmix_state,
368           IMIN(len, ANALYSIS_BUF_SIZE-tonal->mem_fill), offset, c1, c2, C, tonal->Fs);
369     if (tonal->mem_fill+len < ANALYSIS_BUF_SIZE)
370     {
371        tonal->mem_fill += len;
372        /* Don't have enough to update the analysis */
373        RESTORE_STACK;
374        return;
375     }
376     hp_ener = tonal->hp_ener_accum;
377     info = &tonal->info[tonal->write_pos++];
378     if (tonal->write_pos>=DETECT_SIZE)
379        tonal->write_pos-=DETECT_SIZE;
380
381     ALLOC(in, 480, kiss_fft_cpx);
382     ALLOC(out, 480, kiss_fft_cpx);
383     ALLOC(tonality, 240, float);
384     ALLOC(noisiness, 240, float);
385     for (i=0;i<N2;i++)
386     {
387        float w = analysis_window[i];
388        in[i].r = (kiss_fft_scalar)(w*tonal->inmem[i]);
389        in[i].i = (kiss_fft_scalar)(w*tonal->inmem[N2+i]);
390        in[N-i-1].r = (kiss_fft_scalar)(w*tonal->inmem[N-i-1]);
391        in[N-i-1].i = (kiss_fft_scalar)(w*tonal->inmem[N+N2-i-1]);
392     }
393     OPUS_MOVE(tonal->inmem, tonal->inmem+ANALYSIS_BUF_SIZE-240, 240);
394     remaining = len - (ANALYSIS_BUF_SIZE-tonal->mem_fill);
395     tonal->hp_ener_accum = (float)downmix_and_resample(downmix, x,
396           &tonal->inmem[240], tonal->downmix_state, remaining,
397           offset+ANALYSIS_BUF_SIZE-tonal->mem_fill, c1, c2, C, tonal->Fs);
398     tonal->mem_fill = 240 + remaining;
399     opus_fft(kfft, in, out, tonal->arch);
400 #ifndef FIXED_POINT
401     /* If there's any NaN on the input, the entire output will be NaN, so we only need to check one value. */
402     if (celt_isnan(out[0].r))
403     {
404        info->valid = 0;
405        RESTORE_STACK;
406        return;
407     }
408 #endif
409
410     for (i=1;i<N2;i++)
411     {
412        float X1r, X2r, X1i, X2i;
413        float angle, d_angle, d2_angle;
414        float angle2, d_angle2, d2_angle2;
415        float mod1, mod2, avg_mod;
416        X1r = (float)out[i].r+out[N-i].r;
417        X1i = (float)out[i].i-out[N-i].i;
418        X2r = (float)out[i].i+out[N-i].i;
419        X2i = (float)out[N-i].r-out[i].r;
420
421        angle = (float)(.5f/M_PI)*fast_atan2f(X1i, X1r);
422        d_angle = angle - A[i];
423        d2_angle = d_angle - dA[i];
424
425        angle2 = (float)(.5f/M_PI)*fast_atan2f(X2i, X2r);
426        d_angle2 = angle2 - angle;
427        d2_angle2 = d_angle2 - d_angle;
428
429        mod1 = d2_angle - (float)float2int(d2_angle);
430        noisiness[i] = ABS16(mod1);
431        mod1 *= mod1;
432        mod1 *= mod1;
433
434        mod2 = d2_angle2 - (float)float2int(d2_angle2);
435        noisiness[i] += ABS16(mod2);
436        mod2 *= mod2;
437        mod2 *= mod2;
438
439        avg_mod = .25f*(d2A[i]+mod1+2*mod2);
440        /* This introduces an extra delay of 2 frames in the detection. */
441        tonality[i] = 1.f/(1.f+40.f*16.f*pi4*avg_mod)-.015f;
442        /* No delay on this detection, but it's less reliable. */
443        tonality2[i] = 1.f/(1.f+40.f*16.f*pi4*mod2)-.015f;
444
445        A[i] = angle2;
446        dA[i] = d_angle2;
447        d2A[i] = mod2;
448     }
449     for (i=2;i<N2-1;i++)
450     {
451        float tt = MIN32(tonality2[i], MAX32(tonality2[i-1], tonality2[i+1]));
452        tonality[i] = .9f*MAX32(tonality[i], tt-.1f);
453     }
454     frame_tonality = 0;
455     max_frame_tonality = 0;
456     /*tw_sum = 0;*/
457     info->activity = 0;
458     frame_noisiness = 0;
459     frame_stationarity = 0;
460     if (!tonal->count)
461     {
462        for (b=0;b<NB_TBANDS;b++)
463        {
464           tonal->lowE[b] = 1e10;
465           tonal->highE[b] = -1e10;
466        }
467     }
468     relativeE = 0;
469     frame_loudness = 0;
470     /* The energy of the very first band is special because of DC. */
471     {
472        float E = 0;
473        float X1r, X2r;
474        X1r = 2*(float)out[0].r;
475        X2r = 2*(float)out[0].i;
476        E = X1r*X1r + X2r*X2r;
477        for (i=1;i<4;i++)
478        {
479           float binE = out[i].r*(float)out[i].r + out[N-i].r*(float)out[N-i].r
480                      + out[i].i*(float)out[i].i + out[N-i].i*(float)out[N-i].i;
481           E += binE;
482        }
483        E = SCALE_ENER(E);
484        band_log2[0] = (float).5*log2(E+1e-10f);
485     }
486     for (b=0;b<NB_TBANDS;b++)
487     {
488        float E=0, tE=0, nE=0;
489        float L1, L2;
490        float stationarity;
491        for (i=tbands[b];i<tbands[b+1];i++)
492        {
493           float binE = out[i].r*(float)out[i].r + out[N-i].r*(float)out[N-i].r
494                      + out[i].i*(float)out[i].i + out[N-i].i*(float)out[N-i].i;
495           binE = SCALE_ENER(binE);
496           E += binE;
497           tE += binE*MAX32(0, tonality[i]);
498           nE += binE*2.f*(.5f-noisiness[i]);
499        }
500 #ifndef FIXED_POINT
501        /* Check for extreme band energies that could cause NaNs later. */
502        if (!(E<1e9f) || celt_isnan(E))
503        {
504           info->valid = 0;
505           RESTORE_STACK;
506           return;
507        }
508 #endif
509
510        tonal->E[tonal->E_count][b] = E;
511        frame_noisiness += nE/(1e-15f+E);
512
513        frame_loudness += (float)sqrt(E+1e-10f);
514        logE[b] = (float)log(E+1e-10f);
515        band_log2[b+1] = (float).5*log2(E+1e-10f);
516        tonal->logE[tonal->E_count][b] = logE[b];
517        if (tonal->count==0)
518           tonal->highE[b] = tonal->lowE[b] = logE[b];
519        if (tonal->highE[b] > tonal->lowE[b] + 7.5)
520        {
521           if (tonal->highE[b] - logE[b] > logE[b] - tonal->lowE[b])
522              tonal->highE[b] -= .01f;
523           else
524              tonal->lowE[b] += .01f;
525        }
526        if (logE[b] > tonal->highE[b])
527        {
528           tonal->highE[b] = logE[b];
529           tonal->lowE[b] = MAX32(tonal->highE[b]-15, tonal->lowE[b]);
530        } else if (logE[b] < tonal->lowE[b])
531        {
532           tonal->lowE[b] = logE[b];
533           tonal->highE[b] = MIN32(tonal->lowE[b]+15, tonal->highE[b]);
534        }
535        relativeE += (logE[b]-tonal->lowE[b])/(1e-15f + (tonal->highE[b]-tonal->lowE[b]));
536
537        L1=L2=0;
538        for (i=0;i<NB_FRAMES;i++)
539        {
540           L1 += (float)sqrt(tonal->E[i][b]);
541           L2 += tonal->E[i][b];
542        }
543
544        stationarity = MIN16(0.99f,L1/(float)sqrt(1e-15+NB_FRAMES*L2));
545        stationarity *= stationarity;
546        stationarity *= stationarity;
547        frame_stationarity += stationarity;
548        /*band_tonality[b] = tE/(1e-15+E)*/;
549        band_tonality[b] = MAX16(tE/(1e-15f+E), stationarity*tonal->prev_band_tonality[b]);
550 #if 0
551        if (b>=NB_TONAL_SKIP_BANDS)
552        {
553           frame_tonality += tweight[b]*band_tonality[b];
554           tw_sum += tweight[b];
555        }
556 #else
557        frame_tonality += band_tonality[b];
558        if (b>=NB_TBANDS-NB_TONAL_SKIP_BANDS)
559           frame_tonality -= band_tonality[b-NB_TBANDS+NB_TONAL_SKIP_BANDS];
560 #endif
561        max_frame_tonality = MAX16(max_frame_tonality, (1.f+.03f*(b-NB_TBANDS))*frame_tonality);
562        slope += band_tonality[b]*(b-8);
563        /*printf("%f %f ", band_tonality[b], stationarity);*/
564        tonal->prev_band_tonality[b] = band_tonality[b];
565     }
566
567     leakage_from[0] = band_log2[0];
568     leakage_to[0] = band_log2[0] - LEAKAGE_OFFSET;
569     for (b=1;b<NB_TBANDS+1;b++)
570     {
571        float leak_slope = LEAKAGE_SLOPE*(tbands[b]-tbands[b-1])/4;
572        leakage_from[b] = MIN16(leakage_from[b-1]+leak_slope, band_log2[b]);
573        leakage_to[b] = MAX16(leakage_to[b-1]-leak_slope, band_log2[b]-LEAKAGE_OFFSET);
574     }
575     for (b=NB_TBANDS-2;b>=0;b--)
576     {
577        float leak_slope = LEAKAGE_SLOPE*(tbands[b+1]-tbands[b])/4;
578        leakage_from[b] = MIN16(leakage_from[b+1]+leak_slope, leakage_from[b]);
579        leakage_to[b] = MAX16(leakage_to[b+1]-leak_slope, leakage_to[b]);
580     }
581     celt_assert(NB_TBANDS+1 <= LEAK_BANDS);
582     for (b=0;b<NB_TBANDS+1;b++)
583     {
584        /* leak_boost[] is made up of two terms. The first, based on leakage_high[],
585           represents the boost needed to overcome the amount of analysis leakage
586           cause in a weaker band b by louder neighroubing bands.
587           The second, based on leakage_low[], applies to a loud band b for
588           which the quantization noise causes synthesis leakage to the weaker
589           neighbouring bands. */
590        float boost = MAX16(0, leakage_to[b] - band_log2[b]) +
591              MAX16(0, band_log2[b] - (leakage_from[b]+LEAKAGE_OFFSET));
592        info->leak_boost[b] = IMIN(255, (int)floor(.5 + 64.f*boost));
593     }
594     for (;b<LEAK_BANDS;b++) info->leak_boost[b] = 0;
595
596     for (i=0;i<NB_FRAMES;i++)
597     {
598        int j;
599        float mindist = 1e15f;
600        for (j=0;j<NB_FRAMES;j++)
601        {
602           int k;
603           float dist=0;
604           for (k=0;k<NB_TBANDS;k++)
605           {
606              float tmp;
607              tmp = tonal->logE[i][k] - tonal->logE[j][k];
608              dist += tmp*tmp;
609           }
610           if (j!=i)
611              mindist = MIN32(mindist, dist);
612        }
613        spec_variability += mindist;
614     }
615     spec_variability = (float)sqrt(spec_variability/NB_FRAMES/NB_TBANDS);
616     bandwidth_mask = 0;
617     bandwidth = 0;
618     maxE = 0;
619     noise_floor = 5.7e-4f/(1<<(IMAX(0,lsb_depth-8)));
620     noise_floor *= noise_floor;
621     for (b=0;b<NB_TBANDS;b++)
622     {
623        float E=0;
624        int band_start, band_end;
625        /* Keep a margin of 300 Hz for aliasing */
626        band_start = tbands[b];
627        band_end = tbands[b+1];
628        for (i=band_start;i<band_end;i++)
629        {
630           float binE = out[i].r*(float)out[i].r + out[N-i].r*(float)out[N-i].r
631                      + out[i].i*(float)out[i].i + out[N-i].i*(float)out[N-i].i;
632           E += binE;
633        }
634        E = SCALE_ENER(E);
635        maxE = MAX32(maxE, E);
636        tonal->meanE[b] = MAX32((1-alphaE2)*tonal->meanE[b], E);
637        E = MAX32(E, tonal->meanE[b]);
638        /* Use a simple follower with 13 dB/Bark slope for spreading function */
639        bandwidth_mask = MAX32(.05f*bandwidth_mask, E);
640        /* Consider the band "active" only if all these conditions are met:
641           1) less than 10 dB below the simple follower
642           2) less than 90 dB below the peak band (maximal masking possible considering
643              both the ATH and the loudness-dependent slope of the spreading function)
644           3) above the PCM quantization noise floor
645           We use b+1 because the first CELT band isn't included in tbands[]
646        */
647        if (E>.1*bandwidth_mask && E*1e9f > maxE && E > noise_floor*(band_end-band_start))
648           bandwidth = b+1;
649     }
650     /* Special case for the last two bands, for which we don't have spectrum but only
651        the energy above 12 kHz. */
652     {
653        float E = hp_ener*(1.f/(240*240));
654 #ifdef FIXED_POINT
655        /* silk_resampler_down2_hp() shifted right by an extra 8 bits. */
656        E *= 256.f*(1.f/Q15ONE)*(1.f/Q15ONE);
657 #endif
658        maxE = MAX32(maxE, E);
659        tonal->meanE[b] = MAX32((1-alphaE2)*tonal->meanE[b], E);
660        E = MAX32(E, tonal->meanE[b]);
661        /* Use a simple follower with 13 dB/Bark slope for spreading function */
662        bandwidth_mask = MAX32(.05f*bandwidth_mask, E);
663        if (E>.1*bandwidth_mask && E*1e9f > maxE && E > noise_floor*160)
664           bandwidth = 20;
665     }
666     if (tonal->count<=2)
667        bandwidth = 20;
668     frame_loudness = 20*(float)log10(frame_loudness);
669     tonal->Etracker = MAX32(tonal->Etracker-.003f, frame_loudness);
670     tonal->lowECount *= (1-alphaE);
671     if (frame_loudness < tonal->Etracker-30)
672        tonal->lowECount += alphaE;
673
674     for (i=0;i<8;i++)
675     {
676        float sum=0;
677        for (b=0;b<16;b++)
678           sum += dct_table[i*16+b]*logE[b];
679        BFCC[i] = sum;
680     }
681     for (i=0;i<8;i++)
682     {
683        float sum=0;
684        for (b=0;b<16;b++)
685           sum += dct_table[i*16+b]*.5f*(tonal->highE[b]+tonal->lowE[b]);
686        midE[i] = sum;
687     }
688
689     frame_stationarity /= NB_TBANDS;
690     relativeE /= NB_TBANDS;
691     if (tonal->count<10)
692        relativeE = .5;
693     frame_noisiness /= NB_TBANDS;
694 #if 1
695     info->activity = frame_noisiness + (1-frame_noisiness)*relativeE;
696 #else
697     info->activity = .5*(1+frame_noisiness-frame_stationarity);
698 #endif
699     frame_tonality = (max_frame_tonality/(NB_TBANDS-NB_TONAL_SKIP_BANDS));
700     frame_tonality = MAX16(frame_tonality, tonal->prev_tonality*.8f);
701     tonal->prev_tonality = frame_tonality;
702
703     slope /= 8*8;
704     info->tonality_slope = slope;
705
706     tonal->E_count = (tonal->E_count+1)%NB_FRAMES;
707     tonal->count = IMIN(tonal->count+1, ANALYSIS_COUNT_MAX);
708     info->tonality = frame_tonality;
709
710     for (i=0;i<4;i++)
711        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];
712
713     for (i=0;i<4;i++)
714        tonal->cmean[i] = (1-alpha)*tonal->cmean[i] + alpha*BFCC[i];
715
716     for (i=0;i<4;i++)
717         features[4+i] = 0.63246f*(BFCC[i]-tonal->mem[i+24]) + 0.31623f*(tonal->mem[i]-tonal->mem[i+16]);
718     for (i=0;i<3;i++)
719         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];
720
721     if (tonal->count > 5)
722     {
723        for (i=0;i<9;i++)
724           tonal->std[i] = (1-alpha)*tonal->std[i] + alpha*features[i]*features[i];
725     }
726     for (i=0;i<4;i++)
727        features[i] = BFCC[i]-midE[i];
728
729     for (i=0;i<8;i++)
730     {
731        tonal->mem[i+24] = tonal->mem[i+16];
732        tonal->mem[i+16] = tonal->mem[i+8];
733        tonal->mem[i+8] = tonal->mem[i];
734        tonal->mem[i] = BFCC[i];
735     }
736     for (i=0;i<9;i++)
737        features[11+i] = (float)sqrt(tonal->std[i]) - std_feature_bias[i];
738     features[18] = spec_variability - 0.78f;
739     features[20] = info->tonality - 0.154723f;
740     features[21] = info->activity - 0.724643f;
741     features[22] = frame_stationarity - 0.743717f;
742     features[23] = info->tonality_slope + 0.069216f;
743     features[24] = tonal->lowECount - 0.067930f;
744
745     mlp_process(&net, features, frame_probs);
746     frame_probs[0] = .5f*(frame_probs[0]+1);
747     /* Curve fitting between the MLP probability and the actual probability */
748     /*frame_probs[0] = .01f + 1.21f*frame_probs[0]*frame_probs[0] - .23f*(float)pow(frame_probs[0], 10);*/
749     /* Probability of active audio (as opposed to silence) */
750     frame_probs[1] = .5f*frame_probs[1]+.5f;
751     frame_probs[1] *= frame_probs[1];
752
753     /* Probability of speech or music vs noise */
754     info->activity_probability = frame_probs[1];
755
756     /*printf("%f %f\n", frame_probs[0], frame_probs[1]);*/
757     {
758        /* Probability of state transition */
759        float tau;
760        /* Represents independence of the MLP probabilities, where
761           beta=1 means fully independent. */
762        float beta;
763        /* Denormalized probability of speech (p0) and music (p1) after update */
764        float p0, p1;
765        /* Probabilities for "all speech" and "all music" */
766        float s0, m0;
767        /* Probability sum for renormalisation */
768        float psum;
769        /* Instantaneous probability of speech and music, with beta pre-applied. */
770        float speech0;
771        float music0;
772        float p, q;
773
774        /* More silence transitions for speech than for music. */
775        tau = .001f*tonal->music_prob + .01f*(1-tonal->music_prob);
776        p = MAX16(.05f,MIN16(.95f,frame_probs[1]));
777        q = MAX16(.05f,MIN16(.95f,tonal->vad_prob));
778        beta = .02f+.05f*ABS16(p-q)/(p*(1-q)+q*(1-p));
779        /* p0 and p1 are the probabilities of speech and music at this frame
780           using only information from previous frame and applying the
781           state transition model */
782        p0 = (1-tonal->vad_prob)*(1-tau) +    tonal->vad_prob *tau;
783        p1 =    tonal->vad_prob *(1-tau) + (1-tonal->vad_prob)*tau;
784        /* We apply the current probability with exponent beta to work around
785           the fact that the probability estimates aren't independent. */
786        p0 *= (float)pow(1-frame_probs[1], beta);
787        p1 *= (float)pow(frame_probs[1], beta);
788        /* Normalise the probabilities to get the Marokv probability of music. */
789        tonal->vad_prob = p1/(p0+p1);
790        info->vad_prob = tonal->vad_prob;
791        /* Consider that silence has a 50-50 probability of being speech or music. */
792        frame_probs[0] = tonal->vad_prob*frame_probs[0] + (1-tonal->vad_prob)*.5f;
793
794        /* One transition every 3 minutes of active audio */
795        tau = .0001f;
796        /* Adapt beta based on how "unexpected" the new prob is */
797        p = MAX16(.05f,MIN16(.95f,frame_probs[0]));
798        q = MAX16(.05f,MIN16(.95f,tonal->music_prob));
799        beta = .02f+.05f*ABS16(p-q)/(p*(1-q)+q*(1-p));
800        /* p0 and p1 are the probabilities of speech and music at this frame
801           using only information from previous frame and applying the
802           state transition model */
803        p0 = (1-tonal->music_prob)*(1-tau) +    tonal->music_prob *tau;
804        p1 =    tonal->music_prob *(1-tau) + (1-tonal->music_prob)*tau;
805        /* We apply the current probability with exponent beta to work around
806           the fact that the probability estimates aren't independent. */
807        p0 *= (float)pow(1-frame_probs[0], beta);
808        p1 *= (float)pow(frame_probs[0], beta);
809        /* Normalise the probabilities to get the Marokv probability of music. */
810        tonal->music_prob = p1/(p0+p1);
811        info->music_prob = tonal->music_prob;
812
813        /*printf("%f %f %f %f\n", frame_probs[0], frame_probs[1], tonal->music_prob, tonal->vad_prob);*/
814        /* This chunk of code deals with delayed decision. */
815        psum=1e-20f;
816        /* Instantaneous probability of speech and music, with beta pre-applied. */
817        speech0 = (float)pow(1-frame_probs[0], beta);
818        music0  = (float)pow(frame_probs[0], beta);
819        if (tonal->count==1)
820        {
821           tonal->pspeech[0]=.5;
822           tonal->pmusic [0]=.5;
823        }
824        /* Updated probability of having only speech (s0) or only music (m0),
825           before considering the new observation. */
826        s0 = tonal->pspeech[0] + tonal->pspeech[1];
827        m0 = tonal->pmusic [0] + tonal->pmusic [1];
828        /* Updates s0 and m0 with instantaneous probability. */
829        tonal->pspeech[0] = s0*(1-tau)*speech0;
830        tonal->pmusic [0] = m0*(1-tau)*music0;
831        /* Propagate the transition probabilities */
832        for (i=1;i<DETECT_SIZE-1;i++)
833        {
834           tonal->pspeech[i] = tonal->pspeech[i+1]*speech0;
835           tonal->pmusic [i] = tonal->pmusic [i+1]*music0;
836        }
837        /* Probability that the latest frame is speech, when all the previous ones were music. */
838        tonal->pspeech[DETECT_SIZE-1] = m0*tau*speech0;
839        /* Probability that the latest frame is music, when all the previous ones were speech. */
840        tonal->pmusic [DETECT_SIZE-1] = s0*tau*music0;
841
842        /* Renormalise probabilities to 1 */
843        for (i=0;i<DETECT_SIZE;i++)
844           psum += tonal->pspeech[i] + tonal->pmusic[i];
845        psum = 1.f/psum;
846        for (i=0;i<DETECT_SIZE;i++)
847        {
848           tonal->pspeech[i] *= psum;
849           tonal->pmusic [i] *= psum;
850        }
851        psum = tonal->pmusic[0];
852        for (i=1;i<DETECT_SIZE;i++)
853           psum += tonal->pspeech[i];
854
855        /* Estimate our confidence in the speech/music decisions */
856        if (frame_probs[1]>.75)
857        {
858           if (tonal->music_prob>.9)
859           {
860              float adapt;
861              adapt = 1.f/(++tonal->music_confidence_count);
862              tonal->music_confidence_count = IMIN(tonal->music_confidence_count, 500);
863              tonal->music_confidence += adapt*MAX16(-.2f,frame_probs[0]-tonal->music_confidence);
864           }
865           if (tonal->music_prob<.1)
866           {
867              float adapt;
868              adapt = 1.f/(++tonal->speech_confidence_count);
869              tonal->speech_confidence_count = IMIN(tonal->speech_confidence_count, 500);
870              tonal->speech_confidence += adapt*MIN16(.2f,frame_probs[0]-tonal->speech_confidence);
871           }
872        } else {
873           if (tonal->music_confidence_count==0)
874              tonal->music_confidence = .9f;
875           if (tonal->speech_confidence_count==0)
876              tonal->speech_confidence = .1f;
877        }
878     }
879     tonal->last_music = tonal->music_prob>.5f;
880 #ifdef MLP_TRAINING
881     for (i=0;i<25;i++)
882        printf("%f ", features[i]);
883     printf("\n");
884 #endif
885
886     info->bandwidth = bandwidth;
887     /*printf("%d %d\n", info->bandwidth, info->opus_bandwidth);*/
888     info->noisiness = frame_noisiness;
889     info->valid = 1;
890     RESTORE_STACK;
891 }
892
893 void run_analysis(TonalityAnalysisState *analysis, const CELTMode *celt_mode, const void *analysis_pcm,
894                  int analysis_frame_size, int frame_size, int c1, int c2, int C, opus_int32 Fs,
895                  int lsb_depth, downmix_func downmix, AnalysisInfo *analysis_info)
896 {
897    int offset;
898    int pcm_len;
899
900    analysis_frame_size -= analysis_frame_size&1;
901    if (analysis_pcm != NULL)
902    {
903       /* Avoid overflow/wrap-around of the analysis buffer */
904       analysis_frame_size = IMIN((DETECT_SIZE-5)*Fs/50, analysis_frame_size);
905
906       pcm_len = analysis_frame_size - analysis->analysis_offset;
907       offset = analysis->analysis_offset;
908       while (pcm_len>0) {
909          tonality_analysis(analysis, celt_mode, analysis_pcm, IMIN(Fs/50, pcm_len), offset, c1, c2, C, lsb_depth, downmix);
910          offset += Fs/50;
911          pcm_len -= Fs/50;
912       }
913       analysis->analysis_offset = analysis_frame_size;
914
915       analysis->analysis_offset -= frame_size;
916    }
917
918    analysis_info->valid = 0;
919    tonality_get_info(analysis, analysis_info, frame_size);
920 }
921
922 #endif /* DISABLE_FLOAT_API */