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