1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/media/libopus/src/analysis.c Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,645 @@
1.4 +/* Copyright (c) 2011 Xiph.Org Foundation
1.5 + Written by Jean-Marc Valin */
1.6 +/*
1.7 + Redistribution and use in source and binary forms, with or without
1.8 + modification, are permitted provided that the following conditions
1.9 + are met:
1.10 +
1.11 + - Redistributions of source code must retain the above copyright
1.12 + notice, this list of conditions and the following disclaimer.
1.13 +
1.14 + - Redistributions in binary form must reproduce the above copyright
1.15 + notice, this list of conditions and the following disclaimer in the
1.16 + documentation and/or other materials provided with the distribution.
1.17 +
1.18 + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
1.19 + ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
1.20 + LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
1.21 + A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR
1.22 + CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
1.23 + EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
1.24 + PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
1.25 + PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
1.26 + LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
1.27 + NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
1.28 + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
1.29 +*/
1.30 +
1.31 +#ifdef HAVE_CONFIG_H
1.32 +#include "config.h"
1.33 +#endif
1.34 +
1.35 +#include "kiss_fft.h"
1.36 +#include "celt.h"
1.37 +#include "modes.h"
1.38 +#include "arch.h"
1.39 +#include "quant_bands.h"
1.40 +#include <stdio.h>
1.41 +#include "analysis.h"
1.42 +#include "mlp.h"
1.43 +#include "stack_alloc.h"
1.44 +
1.45 +extern const MLP net;
1.46 +
1.47 +#ifndef M_PI
1.48 +#define M_PI 3.141592653
1.49 +#endif
1.50 +
1.51 +static const float dct_table[128] = {
1.52 + 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f,
1.53 + 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f,
1.54 + 0.351851f, 0.338330f, 0.311806f, 0.273300f, 0.224292f, 0.166664f, 0.102631f, 0.034654f,
1.55 + -0.034654f,-0.102631f,-0.166664f,-0.224292f,-0.273300f,-0.311806f,-0.338330f,-0.351851f,
1.56 + 0.346760f, 0.293969f, 0.196424f, 0.068975f,-0.068975f,-0.196424f,-0.293969f,-0.346760f,
1.57 + -0.346760f,-0.293969f,-0.196424f,-0.068975f, 0.068975f, 0.196424f, 0.293969f, 0.346760f,
1.58 + 0.338330f, 0.224292f, 0.034654f,-0.166664f,-0.311806f,-0.351851f,-0.273300f,-0.102631f,
1.59 + 0.102631f, 0.273300f, 0.351851f, 0.311806f, 0.166664f,-0.034654f,-0.224292f,-0.338330f,
1.60 + 0.326641f, 0.135299f,-0.135299f,-0.326641f,-0.326641f,-0.135299f, 0.135299f, 0.326641f,
1.61 + 0.326641f, 0.135299f,-0.135299f,-0.326641f,-0.326641f,-0.135299f, 0.135299f, 0.326641f,
1.62 + 0.311806f, 0.034654f,-0.273300f,-0.338330f,-0.102631f, 0.224292f, 0.351851f, 0.166664f,
1.63 + -0.166664f,-0.351851f,-0.224292f, 0.102631f, 0.338330f, 0.273300f,-0.034654f,-0.311806f,
1.64 + 0.293969f,-0.068975f,-0.346760f,-0.196424f, 0.196424f, 0.346760f, 0.068975f,-0.293969f,
1.65 + -0.293969f, 0.068975f, 0.346760f, 0.196424f,-0.196424f,-0.346760f,-0.068975f, 0.293969f,
1.66 + 0.273300f,-0.166664f,-0.338330f, 0.034654f, 0.351851f, 0.102631f,-0.311806f,-0.224292f,
1.67 + 0.224292f, 0.311806f,-0.102631f,-0.351851f,-0.034654f, 0.338330f, 0.166664f,-0.273300f,
1.68 +};
1.69 +
1.70 +static const float analysis_window[240] = {
1.71 + 0.000043f, 0.000171f, 0.000385f, 0.000685f, 0.001071f, 0.001541f, 0.002098f, 0.002739f,
1.72 + 0.003466f, 0.004278f, 0.005174f, 0.006156f, 0.007222f, 0.008373f, 0.009607f, 0.010926f,
1.73 + 0.012329f, 0.013815f, 0.015385f, 0.017037f, 0.018772f, 0.020590f, 0.022490f, 0.024472f,
1.74 + 0.026535f, 0.028679f, 0.030904f, 0.033210f, 0.035595f, 0.038060f, 0.040604f, 0.043227f,
1.75 + 0.045928f, 0.048707f, 0.051564f, 0.054497f, 0.057506f, 0.060591f, 0.063752f, 0.066987f,
1.76 + 0.070297f, 0.073680f, 0.077136f, 0.080665f, 0.084265f, 0.087937f, 0.091679f, 0.095492f,
1.77 + 0.099373f, 0.103323f, 0.107342f, 0.111427f, 0.115579f, 0.119797f, 0.124080f, 0.128428f,
1.78 + 0.132839f, 0.137313f, 0.141849f, 0.146447f, 0.151105f, 0.155823f, 0.160600f, 0.165435f,
1.79 + 0.170327f, 0.175276f, 0.180280f, 0.185340f, 0.190453f, 0.195619f, 0.200838f, 0.206107f,
1.80 + 0.211427f, 0.216797f, 0.222215f, 0.227680f, 0.233193f, 0.238751f, 0.244353f, 0.250000f,
1.81 + 0.255689f, 0.261421f, 0.267193f, 0.273005f, 0.278856f, 0.284744f, 0.290670f, 0.296632f,
1.82 + 0.302628f, 0.308658f, 0.314721f, 0.320816f, 0.326941f, 0.333097f, 0.339280f, 0.345492f,
1.83 + 0.351729f, 0.357992f, 0.364280f, 0.370590f, 0.376923f, 0.383277f, 0.389651f, 0.396044f,
1.84 + 0.402455f, 0.408882f, 0.415325f, 0.421783f, 0.428254f, 0.434737f, 0.441231f, 0.447736f,
1.85 + 0.454249f, 0.460770f, 0.467298f, 0.473832f, 0.480370f, 0.486912f, 0.493455f, 0.500000f,
1.86 + 0.506545f, 0.513088f, 0.519630f, 0.526168f, 0.532702f, 0.539230f, 0.545751f, 0.552264f,
1.87 + 0.558769f, 0.565263f, 0.571746f, 0.578217f, 0.584675f, 0.591118f, 0.597545f, 0.603956f,
1.88 + 0.610349f, 0.616723f, 0.623077f, 0.629410f, 0.635720f, 0.642008f, 0.648271f, 0.654508f,
1.89 + 0.660720f, 0.666903f, 0.673059f, 0.679184f, 0.685279f, 0.691342f, 0.697372f, 0.703368f,
1.90 + 0.709330f, 0.715256f, 0.721144f, 0.726995f, 0.732807f, 0.738579f, 0.744311f, 0.750000f,
1.91 + 0.755647f, 0.761249f, 0.766807f, 0.772320f, 0.777785f, 0.783203f, 0.788573f, 0.793893f,
1.92 + 0.799162f, 0.804381f, 0.809547f, 0.814660f, 0.819720f, 0.824724f, 0.829673f, 0.834565f,
1.93 + 0.839400f, 0.844177f, 0.848895f, 0.853553f, 0.858151f, 0.862687f, 0.867161f, 0.871572f,
1.94 + 0.875920f, 0.880203f, 0.884421f, 0.888573f, 0.892658f, 0.896677f, 0.900627f, 0.904508f,
1.95 + 0.908321f, 0.912063f, 0.915735f, 0.919335f, 0.922864f, 0.926320f, 0.929703f, 0.933013f,
1.96 + 0.936248f, 0.939409f, 0.942494f, 0.945503f, 0.948436f, 0.951293f, 0.954072f, 0.956773f,
1.97 + 0.959396f, 0.961940f, 0.964405f, 0.966790f, 0.969096f, 0.971321f, 0.973465f, 0.975528f,
1.98 + 0.977510f, 0.979410f, 0.981228f, 0.982963f, 0.984615f, 0.986185f, 0.987671f, 0.989074f,
1.99 + 0.990393f, 0.991627f, 0.992778f, 0.993844f, 0.994826f, 0.995722f, 0.996534f, 0.997261f,
1.100 + 0.997902f, 0.998459f, 0.998929f, 0.999315f, 0.999615f, 0.999829f, 0.999957f, 1.000000f,
1.101 +};
1.102 +
1.103 +static const int tbands[NB_TBANDS+1] = {
1.104 + 2, 4, 6, 8, 10, 12, 14, 16, 20, 24, 28, 32, 40, 48, 56, 68, 80, 96, 120
1.105 +};
1.106 +
1.107 +static const int extra_bands[NB_TOT_BANDS+1] = {
1.108 + 1, 2, 4, 6, 8, 10, 12, 14, 16, 20, 24, 28, 32, 40, 48, 56, 68, 80, 96, 120, 160, 200
1.109 +};
1.110 +
1.111 +/*static const float tweight[NB_TBANDS+1] = {
1.112 + .3, .4, .5, .6, .7, .8, .9, 1., 1., 1., 1., 1., 1., 1., .8, .7, .6, .5
1.113 +};*/
1.114 +
1.115 +#define NB_TONAL_SKIP_BANDS 9
1.116 +
1.117 +#define cA 0.43157974f
1.118 +#define cB 0.67848403f
1.119 +#define cC 0.08595542f
1.120 +#define cE ((float)M_PI/2)
1.121 +static OPUS_INLINE float fast_atan2f(float y, float x) {
1.122 + float x2, y2;
1.123 + /* Should avoid underflow on the values we'll get */
1.124 + if (ABS16(x)+ABS16(y)<1e-9f)
1.125 + {
1.126 + x*=1e12f;
1.127 + y*=1e12f;
1.128 + }
1.129 + x2 = x*x;
1.130 + y2 = y*y;
1.131 + if(x2<y2){
1.132 + float den = (y2 + cB*x2) * (y2 + cC*x2);
1.133 + if (den!=0)
1.134 + return -x*y*(y2 + cA*x2) / den + (y<0 ? -cE : cE);
1.135 + else
1.136 + return (y<0 ? -cE : cE);
1.137 + }else{
1.138 + float den = (x2 + cB*y2) * (x2 + cC*y2);
1.139 + if (den!=0)
1.140 + return x*y*(x2 + cA*y2) / den + (y<0 ? -cE : cE) - (x*y<0 ? -cE : cE);
1.141 + else
1.142 + return (y<0 ? -cE : cE) - (x*y<0 ? -cE : cE);
1.143 + }
1.144 +}
1.145 +
1.146 +void tonality_get_info(TonalityAnalysisState *tonal, AnalysisInfo *info_out, int len)
1.147 +{
1.148 + int pos;
1.149 + int curr_lookahead;
1.150 + float psum;
1.151 + int i;
1.152 +
1.153 + pos = tonal->read_pos;
1.154 + curr_lookahead = tonal->write_pos-tonal->read_pos;
1.155 + if (curr_lookahead<0)
1.156 + curr_lookahead += DETECT_SIZE;
1.157 +
1.158 + if (len > 480 && pos != tonal->write_pos)
1.159 + {
1.160 + pos++;
1.161 + if (pos==DETECT_SIZE)
1.162 + pos=0;
1.163 + }
1.164 + if (pos == tonal->write_pos)
1.165 + pos--;
1.166 + if (pos<0)
1.167 + pos = DETECT_SIZE-1;
1.168 + OPUS_COPY(info_out, &tonal->info[pos], 1);
1.169 + tonal->read_subframe += len/120;
1.170 + while (tonal->read_subframe>=4)
1.171 + {
1.172 + tonal->read_subframe -= 4;
1.173 + tonal->read_pos++;
1.174 + }
1.175 + if (tonal->read_pos>=DETECT_SIZE)
1.176 + tonal->read_pos-=DETECT_SIZE;
1.177 +
1.178 + /* Compensate for the delay in the features themselves.
1.179 + FIXME: Need a better estimate the 10 I just made up */
1.180 + curr_lookahead = IMAX(curr_lookahead-10, 0);
1.181 +
1.182 + psum=0;
1.183 + /* Summing the probability of transition patterns that involve music at
1.184 + time (DETECT_SIZE-curr_lookahead-1) */
1.185 + for (i=0;i<DETECT_SIZE-curr_lookahead;i++)
1.186 + psum += tonal->pmusic[i];
1.187 + for (;i<DETECT_SIZE;i++)
1.188 + psum += tonal->pspeech[i];
1.189 + psum = psum*tonal->music_confidence + (1-psum)*tonal->speech_confidence;
1.190 + /*printf("%f %f %f\n", psum, info_out->music_prob, info_out->tonality);*/
1.191 +
1.192 + info_out->music_prob = psum;
1.193 +}
1.194 +
1.195 +void tonality_analysis(TonalityAnalysisState *tonal, AnalysisInfo *info_out, const CELTMode *celt_mode, const void *x, int len, int offset, int c1, int c2, int C, int lsb_depth, downmix_func downmix)
1.196 +{
1.197 + int i, b;
1.198 + const kiss_fft_state *kfft;
1.199 + VARDECL(kiss_fft_cpx, in);
1.200 + VARDECL(kiss_fft_cpx, out);
1.201 + int N = 480, N2=240;
1.202 + float * OPUS_RESTRICT A = tonal->angle;
1.203 + float * OPUS_RESTRICT dA = tonal->d_angle;
1.204 + float * OPUS_RESTRICT d2A = tonal->d2_angle;
1.205 + VARDECL(float, tonality);
1.206 + VARDECL(float, noisiness);
1.207 + float band_tonality[NB_TBANDS];
1.208 + float logE[NB_TBANDS];
1.209 + float BFCC[8];
1.210 + float features[25];
1.211 + float frame_tonality;
1.212 + float max_frame_tonality;
1.213 + /*float tw_sum=0;*/
1.214 + float frame_noisiness;
1.215 + const float pi4 = (float)(M_PI*M_PI*M_PI*M_PI);
1.216 + float slope=0;
1.217 + float frame_stationarity;
1.218 + float relativeE;
1.219 + float frame_probs[2];
1.220 + float alpha, alphaE, alphaE2;
1.221 + float frame_loudness;
1.222 + float bandwidth_mask;
1.223 + int bandwidth=0;
1.224 + float maxE = 0;
1.225 + float noise_floor;
1.226 + int remaining;
1.227 + AnalysisInfo *info;
1.228 + SAVE_STACK;
1.229 +
1.230 + tonal->last_transition++;
1.231 + alpha = 1.f/IMIN(20, 1+tonal->count);
1.232 + alphaE = 1.f/IMIN(50, 1+tonal->count);
1.233 + alphaE2 = 1.f/IMIN(1000, 1+tonal->count);
1.234 +
1.235 + if (tonal->count<4)
1.236 + tonal->music_prob = .5;
1.237 + kfft = celt_mode->mdct.kfft[0];
1.238 + if (tonal->count==0)
1.239 + tonal->mem_fill = 240;
1.240 + downmix(x, &tonal->inmem[tonal->mem_fill], IMIN(len, ANALYSIS_BUF_SIZE-tonal->mem_fill), offset, c1, c2, C);
1.241 + if (tonal->mem_fill+len < ANALYSIS_BUF_SIZE)
1.242 + {
1.243 + tonal->mem_fill += len;
1.244 + /* Don't have enough to update the analysis */
1.245 + RESTORE_STACK;
1.246 + return;
1.247 + }
1.248 + info = &tonal->info[tonal->write_pos++];
1.249 + if (tonal->write_pos>=DETECT_SIZE)
1.250 + tonal->write_pos-=DETECT_SIZE;
1.251 +
1.252 + ALLOC(in, 480, kiss_fft_cpx);
1.253 + ALLOC(out, 480, kiss_fft_cpx);
1.254 + ALLOC(tonality, 240, float);
1.255 + ALLOC(noisiness, 240, float);
1.256 + for (i=0;i<N2;i++)
1.257 + {
1.258 + float w = analysis_window[i];
1.259 + in[i].r = (kiss_fft_scalar)(w*tonal->inmem[i]);
1.260 + in[i].i = (kiss_fft_scalar)(w*tonal->inmem[N2+i]);
1.261 + in[N-i-1].r = (kiss_fft_scalar)(w*tonal->inmem[N-i-1]);
1.262 + in[N-i-1].i = (kiss_fft_scalar)(w*tonal->inmem[N+N2-i-1]);
1.263 + }
1.264 + OPUS_MOVE(tonal->inmem, tonal->inmem+ANALYSIS_BUF_SIZE-240, 240);
1.265 + remaining = len - (ANALYSIS_BUF_SIZE-tonal->mem_fill);
1.266 + downmix(x, &tonal->inmem[240], remaining, offset+ANALYSIS_BUF_SIZE-tonal->mem_fill, c1, c2, C);
1.267 + tonal->mem_fill = 240 + remaining;
1.268 + opus_fft(kfft, in, out);
1.269 +
1.270 + for (i=1;i<N2;i++)
1.271 + {
1.272 + float X1r, X2r, X1i, X2i;
1.273 + float angle, d_angle, d2_angle;
1.274 + float angle2, d_angle2, d2_angle2;
1.275 + float mod1, mod2, avg_mod;
1.276 + X1r = (float)out[i].r+out[N-i].r;
1.277 + X1i = (float)out[i].i-out[N-i].i;
1.278 + X2r = (float)out[i].i+out[N-i].i;
1.279 + X2i = (float)out[N-i].r-out[i].r;
1.280 +
1.281 + angle = (float)(.5f/M_PI)*fast_atan2f(X1i, X1r);
1.282 + d_angle = angle - A[i];
1.283 + d2_angle = d_angle - dA[i];
1.284 +
1.285 + angle2 = (float)(.5f/M_PI)*fast_atan2f(X2i, X2r);
1.286 + d_angle2 = angle2 - angle;
1.287 + d2_angle2 = d_angle2 - d_angle;
1.288 +
1.289 + mod1 = d2_angle - (float)floor(.5+d2_angle);
1.290 + noisiness[i] = ABS16(mod1);
1.291 + mod1 *= mod1;
1.292 + mod1 *= mod1;
1.293 +
1.294 + mod2 = d2_angle2 - (float)floor(.5+d2_angle2);
1.295 + noisiness[i] += ABS16(mod2);
1.296 + mod2 *= mod2;
1.297 + mod2 *= mod2;
1.298 +
1.299 + avg_mod = .25f*(d2A[i]+2.f*mod1+mod2);
1.300 + tonality[i] = 1.f/(1.f+40.f*16.f*pi4*avg_mod)-.015f;
1.301 +
1.302 + A[i] = angle2;
1.303 + dA[i] = d_angle2;
1.304 + d2A[i] = mod2;
1.305 + }
1.306 +
1.307 + frame_tonality = 0;
1.308 + max_frame_tonality = 0;
1.309 + /*tw_sum = 0;*/
1.310 + info->activity = 0;
1.311 + frame_noisiness = 0;
1.312 + frame_stationarity = 0;
1.313 + if (!tonal->count)
1.314 + {
1.315 + for (b=0;b<NB_TBANDS;b++)
1.316 + {
1.317 + tonal->lowE[b] = 1e10;
1.318 + tonal->highE[b] = -1e10;
1.319 + }
1.320 + }
1.321 + relativeE = 0;
1.322 + frame_loudness = 0;
1.323 + for (b=0;b<NB_TBANDS;b++)
1.324 + {
1.325 + float E=0, tE=0, nE=0;
1.326 + float L1, L2;
1.327 + float stationarity;
1.328 + for (i=tbands[b];i<tbands[b+1];i++)
1.329 + {
1.330 + float binE = out[i].r*(float)out[i].r + out[N-i].r*(float)out[N-i].r
1.331 + + out[i].i*(float)out[i].i + out[N-i].i*(float)out[N-i].i;
1.332 +#ifdef FIXED_POINT
1.333 + /* FIXME: It's probably best to change the BFCC filter initial state instead */
1.334 + binE *= 5.55e-17f;
1.335 +#endif
1.336 + E += binE;
1.337 + tE += binE*tonality[i];
1.338 + nE += binE*2.f*(.5f-noisiness[i]);
1.339 + }
1.340 + tonal->E[tonal->E_count][b] = E;
1.341 + frame_noisiness += nE/(1e-15f+E);
1.342 +
1.343 + frame_loudness += (float)sqrt(E+1e-10f);
1.344 + logE[b] = (float)log(E+1e-10f);
1.345 + tonal->lowE[b] = MIN32(logE[b], tonal->lowE[b]+.01f);
1.346 + tonal->highE[b] = MAX32(logE[b], tonal->highE[b]-.1f);
1.347 + if (tonal->highE[b] < tonal->lowE[b]+1.f)
1.348 + {
1.349 + tonal->highE[b]+=.5f;
1.350 + tonal->lowE[b]-=.5f;
1.351 + }
1.352 + relativeE += (logE[b]-tonal->lowE[b])/(1e-15f+tonal->highE[b]-tonal->lowE[b]);
1.353 +
1.354 + L1=L2=0;
1.355 + for (i=0;i<NB_FRAMES;i++)
1.356 + {
1.357 + L1 += (float)sqrt(tonal->E[i][b]);
1.358 + L2 += tonal->E[i][b];
1.359 + }
1.360 +
1.361 + stationarity = MIN16(0.99f,L1/(float)sqrt(1e-15+NB_FRAMES*L2));
1.362 + stationarity *= stationarity;
1.363 + stationarity *= stationarity;
1.364 + frame_stationarity += stationarity;
1.365 + /*band_tonality[b] = tE/(1e-15+E)*/;
1.366 + band_tonality[b] = MAX16(tE/(1e-15f+E), stationarity*tonal->prev_band_tonality[b]);
1.367 +#if 0
1.368 + if (b>=NB_TONAL_SKIP_BANDS)
1.369 + {
1.370 + frame_tonality += tweight[b]*band_tonality[b];
1.371 + tw_sum += tweight[b];
1.372 + }
1.373 +#else
1.374 + frame_tonality += band_tonality[b];
1.375 + if (b>=NB_TBANDS-NB_TONAL_SKIP_BANDS)
1.376 + frame_tonality -= band_tonality[b-NB_TBANDS+NB_TONAL_SKIP_BANDS];
1.377 +#endif
1.378 + max_frame_tonality = MAX16(max_frame_tonality, (1.f+.03f*(b-NB_TBANDS))*frame_tonality);
1.379 + slope += band_tonality[b]*(b-8);
1.380 + /*printf("%f %f ", band_tonality[b], stationarity);*/
1.381 + tonal->prev_band_tonality[b] = band_tonality[b];
1.382 + }
1.383 +
1.384 + bandwidth_mask = 0;
1.385 + bandwidth = 0;
1.386 + maxE = 0;
1.387 + noise_floor = 5.7e-4f/(1<<(IMAX(0,lsb_depth-8)));
1.388 +#ifdef FIXED_POINT
1.389 + noise_floor *= 1<<(15+SIG_SHIFT);
1.390 +#endif
1.391 + noise_floor *= noise_floor;
1.392 + for (b=0;b<NB_TOT_BANDS;b++)
1.393 + {
1.394 + float E=0;
1.395 + int band_start, band_end;
1.396 + /* Keep a margin of 300 Hz for aliasing */
1.397 + band_start = extra_bands[b];
1.398 + band_end = extra_bands[b+1];
1.399 + for (i=band_start;i<band_end;i++)
1.400 + {
1.401 + float binE = out[i].r*(float)out[i].r + out[N-i].r*(float)out[N-i].r
1.402 + + out[i].i*(float)out[i].i + out[N-i].i*(float)out[N-i].i;
1.403 + E += binE;
1.404 + }
1.405 + maxE = MAX32(maxE, E);
1.406 + tonal->meanE[b] = MAX32((1-alphaE2)*tonal->meanE[b], E);
1.407 + E = MAX32(E, tonal->meanE[b]);
1.408 + /* Use a simple follower with 13 dB/Bark slope for spreading function */
1.409 + bandwidth_mask = MAX32(.05f*bandwidth_mask, E);
1.410 + /* Consider the band "active" only if all these conditions are met:
1.411 + 1) less than 10 dB below the simple follower
1.412 + 2) less than 90 dB below the peak band (maximal masking possible considering
1.413 + both the ATH and the loudness-dependent slope of the spreading function)
1.414 + 3) above the PCM quantization noise floor
1.415 + */
1.416 + if (E>.1*bandwidth_mask && E*1e9f > maxE && E > noise_floor*(band_end-band_start))
1.417 + bandwidth = b;
1.418 + }
1.419 + if (tonal->count<=2)
1.420 + bandwidth = 20;
1.421 + frame_loudness = 20*(float)log10(frame_loudness);
1.422 + tonal->Etracker = MAX32(tonal->Etracker-.03f, frame_loudness);
1.423 + tonal->lowECount *= (1-alphaE);
1.424 + if (frame_loudness < tonal->Etracker-30)
1.425 + tonal->lowECount += alphaE;
1.426 +
1.427 + for (i=0;i<8;i++)
1.428 + {
1.429 + float sum=0;
1.430 + for (b=0;b<16;b++)
1.431 + sum += dct_table[i*16+b]*logE[b];
1.432 + BFCC[i] = sum;
1.433 + }
1.434 +
1.435 + frame_stationarity /= NB_TBANDS;
1.436 + relativeE /= NB_TBANDS;
1.437 + if (tonal->count<10)
1.438 + relativeE = .5;
1.439 + frame_noisiness /= NB_TBANDS;
1.440 +#if 1
1.441 + info->activity = frame_noisiness + (1-frame_noisiness)*relativeE;
1.442 +#else
1.443 + info->activity = .5*(1+frame_noisiness-frame_stationarity);
1.444 +#endif
1.445 + frame_tonality = (max_frame_tonality/(NB_TBANDS-NB_TONAL_SKIP_BANDS));
1.446 + frame_tonality = MAX16(frame_tonality, tonal->prev_tonality*.8f);
1.447 + tonal->prev_tonality = frame_tonality;
1.448 +
1.449 + slope /= 8*8;
1.450 + info->tonality_slope = slope;
1.451 +
1.452 + tonal->E_count = (tonal->E_count+1)%NB_FRAMES;
1.453 + tonal->count++;
1.454 + info->tonality = frame_tonality;
1.455 +
1.456 + for (i=0;i<4;i++)
1.457 + 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];
1.458 +
1.459 + for (i=0;i<4;i++)
1.460 + tonal->cmean[i] = (1-alpha)*tonal->cmean[i] + alpha*BFCC[i];
1.461 +
1.462 + for (i=0;i<4;i++)
1.463 + features[4+i] = 0.63246f*(BFCC[i]-tonal->mem[i+24]) + 0.31623f*(tonal->mem[i]-tonal->mem[i+16]);
1.464 + for (i=0;i<3;i++)
1.465 + 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];
1.466 +
1.467 + if (tonal->count > 5)
1.468 + {
1.469 + for (i=0;i<9;i++)
1.470 + tonal->std[i] = (1-alpha)*tonal->std[i] + alpha*features[i]*features[i];
1.471 + }
1.472 +
1.473 + for (i=0;i<8;i++)
1.474 + {
1.475 + tonal->mem[i+24] = tonal->mem[i+16];
1.476 + tonal->mem[i+16] = tonal->mem[i+8];
1.477 + tonal->mem[i+8] = tonal->mem[i];
1.478 + tonal->mem[i] = BFCC[i];
1.479 + }
1.480 + for (i=0;i<9;i++)
1.481 + features[11+i] = (float)sqrt(tonal->std[i]);
1.482 + features[20] = info->tonality;
1.483 + features[21] = info->activity;
1.484 + features[22] = frame_stationarity;
1.485 + features[23] = info->tonality_slope;
1.486 + features[24] = tonal->lowECount;
1.487 +
1.488 +#ifndef DISABLE_FLOAT_API
1.489 + mlp_process(&net, features, frame_probs);
1.490 + frame_probs[0] = .5f*(frame_probs[0]+1);
1.491 + /* Curve fitting between the MLP probability and the actual probability */
1.492 + frame_probs[0] = .01f + 1.21f*frame_probs[0]*frame_probs[0] - .23f*(float)pow(frame_probs[0], 10);
1.493 + /* Probability of active audio (as opposed to silence) */
1.494 + frame_probs[1] = .5f*frame_probs[1]+.5f;
1.495 + /* Consider that silence has a 50-50 probability. */
1.496 + frame_probs[0] = frame_probs[1]*frame_probs[0] + (1-frame_probs[1])*.5f;
1.497 +
1.498 + /*printf("%f %f ", frame_probs[0], frame_probs[1]);*/
1.499 + {
1.500 + /* Probability of state transition */
1.501 + float tau;
1.502 + /* Represents independence of the MLP probabilities, where
1.503 + beta=1 means fully independent. */
1.504 + float beta;
1.505 + /* Denormalized probability of speech (p0) and music (p1) after update */
1.506 + float p0, p1;
1.507 + /* Probabilities for "all speech" and "all music" */
1.508 + float s0, m0;
1.509 + /* Probability sum for renormalisation */
1.510 + float psum;
1.511 + /* Instantaneous probability of speech and music, with beta pre-applied. */
1.512 + float speech0;
1.513 + float music0;
1.514 +
1.515 + /* One transition every 3 minutes of active audio */
1.516 + tau = .00005f*frame_probs[1];
1.517 + beta = .05f;
1.518 + if (1) {
1.519 + /* Adapt beta based on how "unexpected" the new prob is */
1.520 + float p, q;
1.521 + p = MAX16(.05f,MIN16(.95f,frame_probs[0]));
1.522 + q = MAX16(.05f,MIN16(.95f,tonal->music_prob));
1.523 + beta = .01f+.05f*ABS16(p-q)/(p*(1-q)+q*(1-p));
1.524 + }
1.525 + /* p0 and p1 are the probabilities of speech and music at this frame
1.526 + using only information from previous frame and applying the
1.527 + state transition model */
1.528 + p0 = (1-tonal->music_prob)*(1-tau) + tonal->music_prob *tau;
1.529 + p1 = tonal->music_prob *(1-tau) + (1-tonal->music_prob)*tau;
1.530 + /* We apply the current probability with exponent beta to work around
1.531 + the fact that the probability estimates aren't independent. */
1.532 + p0 *= (float)pow(1-frame_probs[0], beta);
1.533 + p1 *= (float)pow(frame_probs[0], beta);
1.534 + /* Normalise the probabilities to get the Marokv probability of music. */
1.535 + tonal->music_prob = p1/(p0+p1);
1.536 + info->music_prob = tonal->music_prob;
1.537 +
1.538 + /* This chunk of code deals with delayed decision. */
1.539 + psum=1e-20f;
1.540 + /* Instantaneous probability of speech and music, with beta pre-applied. */
1.541 + speech0 = (float)pow(1-frame_probs[0], beta);
1.542 + music0 = (float)pow(frame_probs[0], beta);
1.543 + if (tonal->count==1)
1.544 + {
1.545 + tonal->pspeech[0]=.5;
1.546 + tonal->pmusic [0]=.5;
1.547 + }
1.548 + /* Updated probability of having only speech (s0) or only music (m0),
1.549 + before considering the new observation. */
1.550 + s0 = tonal->pspeech[0] + tonal->pspeech[1];
1.551 + m0 = tonal->pmusic [0] + tonal->pmusic [1];
1.552 + /* Updates s0 and m0 with instantaneous probability. */
1.553 + tonal->pspeech[0] = s0*(1-tau)*speech0;
1.554 + tonal->pmusic [0] = m0*(1-tau)*music0;
1.555 + /* Propagate the transition probabilities */
1.556 + for (i=1;i<DETECT_SIZE-1;i++)
1.557 + {
1.558 + tonal->pspeech[i] = tonal->pspeech[i+1]*speech0;
1.559 + tonal->pmusic [i] = tonal->pmusic [i+1]*music0;
1.560 + }
1.561 + /* Probability that the latest frame is speech, when all the previous ones were music. */
1.562 + tonal->pspeech[DETECT_SIZE-1] = m0*tau*speech0;
1.563 + /* Probability that the latest frame is music, when all the previous ones were speech. */
1.564 + tonal->pmusic [DETECT_SIZE-1] = s0*tau*music0;
1.565 +
1.566 + /* Renormalise probabilities to 1 */
1.567 + for (i=0;i<DETECT_SIZE;i++)
1.568 + psum += tonal->pspeech[i] + tonal->pmusic[i];
1.569 + psum = 1.f/psum;
1.570 + for (i=0;i<DETECT_SIZE;i++)
1.571 + {
1.572 + tonal->pspeech[i] *= psum;
1.573 + tonal->pmusic [i] *= psum;
1.574 + }
1.575 + psum = tonal->pmusic[0];
1.576 + for (i=1;i<DETECT_SIZE;i++)
1.577 + psum += tonal->pspeech[i];
1.578 +
1.579 + /* Estimate our confidence in the speech/music decisions */
1.580 + if (frame_probs[1]>.75)
1.581 + {
1.582 + if (tonal->music_prob>.9)
1.583 + {
1.584 + float adapt;
1.585 + adapt = 1.f/(++tonal->music_confidence_count);
1.586 + tonal->music_confidence_count = IMIN(tonal->music_confidence_count, 500);
1.587 + tonal->music_confidence += adapt*MAX16(-.2f,frame_probs[0]-tonal->music_confidence);
1.588 + }
1.589 + if (tonal->music_prob<.1)
1.590 + {
1.591 + float adapt;
1.592 + adapt = 1.f/(++tonal->speech_confidence_count);
1.593 + tonal->speech_confidence_count = IMIN(tonal->speech_confidence_count, 500);
1.594 + tonal->speech_confidence += adapt*MIN16(.2f,frame_probs[0]-tonal->speech_confidence);
1.595 + }
1.596 + } else {
1.597 + if (tonal->music_confidence_count==0)
1.598 + tonal->music_confidence = .9f;
1.599 + if (tonal->speech_confidence_count==0)
1.600 + tonal->speech_confidence = .1f;
1.601 + }
1.602 + }
1.603 + if (tonal->last_music != (tonal->music_prob>.5f))
1.604 + tonal->last_transition=0;
1.605 + tonal->last_music = tonal->music_prob>.5f;
1.606 +#else
1.607 + info->music_prob = 0;
1.608 +#endif
1.609 + /*for (i=0;i<25;i++)
1.610 + printf("%f ", features[i]);
1.611 + printf("\n");*/
1.612 +
1.613 + info->bandwidth = bandwidth;
1.614 + /*printf("%d %d\n", info->bandwidth, info->opus_bandwidth);*/
1.615 + info->noisiness = frame_noisiness;
1.616 + info->valid = 1;
1.617 + if (info_out!=NULL)
1.618 + OPUS_COPY(info_out, info, 1);
1.619 + RESTORE_STACK;
1.620 +}
1.621 +
1.622 +void run_analysis(TonalityAnalysisState *analysis, const CELTMode *celt_mode, const void *analysis_pcm,
1.623 + int analysis_frame_size, int frame_size, int c1, int c2, int C, opus_int32 Fs,
1.624 + int lsb_depth, downmix_func downmix, AnalysisInfo *analysis_info)
1.625 +{
1.626 + int offset;
1.627 + int pcm_len;
1.628 +
1.629 + if (analysis_pcm != NULL)
1.630 + {
1.631 + /* Avoid overflow/wrap-around of the analysis buffer */
1.632 + analysis_frame_size = IMIN((DETECT_SIZE-5)*Fs/100, analysis_frame_size);
1.633 +
1.634 + pcm_len = analysis_frame_size - analysis->analysis_offset;
1.635 + offset = analysis->analysis_offset;
1.636 + do {
1.637 + tonality_analysis(analysis, NULL, celt_mode, analysis_pcm, IMIN(480, pcm_len), offset, c1, c2, C, lsb_depth, downmix);
1.638 + offset += 480;
1.639 + pcm_len -= 480;
1.640 + } while (pcm_len>0);
1.641 + analysis->analysis_offset = analysis_frame_size;
1.642 +
1.643 + analysis->analysis_offset -= frame_size;
1.644 + }
1.645 +
1.646 + analysis_info->valid = 0;
1.647 + tonality_get_info(analysis, analysis_info, frame_size);
1.648 +}
