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