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