media/libopus/src/analysis.c

Wed, 31 Dec 2014 06:09:35 +0100

author
Michael Schloh von Bennewitz <michael@schloh.com>
date
Wed, 31 Dec 2014 06:09:35 +0100
changeset 0
6474c204b198
permissions
-rw-r--r--

Cloned upstream origin tor-browser at tor-browser-31.3.0esr-4.5-1-build1
revision ID fc1c9ff7c1b2defdbc039f12214767608f46423f for hacking purpose.

michael@0 1 /* Copyright (c) 2011 Xiph.Org Foundation
michael@0 2 Written by Jean-Marc Valin */
michael@0 3 /*
michael@0 4 Redistribution and use in source and binary forms, with or without
michael@0 5 modification, are permitted provided that the following conditions
michael@0 6 are met:
michael@0 7
michael@0 8 - Redistributions of source code must retain the above copyright
michael@0 9 notice, this list of conditions and the following disclaimer.
michael@0 10
michael@0 11 - Redistributions in binary form must reproduce the above copyright
michael@0 12 notice, this list of conditions and the following disclaimer in the
michael@0 13 documentation and/or other materials provided with the distribution.
michael@0 14
michael@0 15 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
michael@0 16 ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
michael@0 17 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
michael@0 18 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR
michael@0 19 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
michael@0 20 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
michael@0 21 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
michael@0 22 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
michael@0 23 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
michael@0 24 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
michael@0 25 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
michael@0 26 */
michael@0 27
michael@0 28 #ifdef HAVE_CONFIG_H
michael@0 29 #include "config.h"
michael@0 30 #endif
michael@0 31
michael@0 32 #include "kiss_fft.h"
michael@0 33 #include "celt.h"
michael@0 34 #include "modes.h"
michael@0 35 #include "arch.h"
michael@0 36 #include "quant_bands.h"
michael@0 37 #include <stdio.h>
michael@0 38 #include "analysis.h"
michael@0 39 #include "mlp.h"
michael@0 40 #include "stack_alloc.h"
michael@0 41
michael@0 42 extern const MLP net;
michael@0 43
michael@0 44 #ifndef M_PI
michael@0 45 #define M_PI 3.141592653
michael@0 46 #endif
michael@0 47
michael@0 48 static const float dct_table[128] = {
michael@0 49 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f,
michael@0 50 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f, 0.250000f,
michael@0 51 0.351851f, 0.338330f, 0.311806f, 0.273300f, 0.224292f, 0.166664f, 0.102631f, 0.034654f,
michael@0 52 -0.034654f,-0.102631f,-0.166664f,-0.224292f,-0.273300f,-0.311806f,-0.338330f,-0.351851f,
michael@0 53 0.346760f, 0.293969f, 0.196424f, 0.068975f,-0.068975f,-0.196424f,-0.293969f,-0.346760f,
michael@0 54 -0.346760f,-0.293969f,-0.196424f,-0.068975f, 0.068975f, 0.196424f, 0.293969f, 0.346760f,
michael@0 55 0.338330f, 0.224292f, 0.034654f,-0.166664f,-0.311806f,-0.351851f,-0.273300f,-0.102631f,
michael@0 56 0.102631f, 0.273300f, 0.351851f, 0.311806f, 0.166664f,-0.034654f,-0.224292f,-0.338330f,
michael@0 57 0.326641f, 0.135299f,-0.135299f,-0.326641f,-0.326641f,-0.135299f, 0.135299f, 0.326641f,
michael@0 58 0.326641f, 0.135299f,-0.135299f,-0.326641f,-0.326641f,-0.135299f, 0.135299f, 0.326641f,
michael@0 59 0.311806f, 0.034654f,-0.273300f,-0.338330f,-0.102631f, 0.224292f, 0.351851f, 0.166664f,
michael@0 60 -0.166664f,-0.351851f,-0.224292f, 0.102631f, 0.338330f, 0.273300f,-0.034654f,-0.311806f,
michael@0 61 0.293969f,-0.068975f,-0.346760f,-0.196424f, 0.196424f, 0.346760f, 0.068975f,-0.293969f,
michael@0 62 -0.293969f, 0.068975f, 0.346760f, 0.196424f,-0.196424f,-0.346760f,-0.068975f, 0.293969f,
michael@0 63 0.273300f,-0.166664f,-0.338330f, 0.034654f, 0.351851f, 0.102631f,-0.311806f,-0.224292f,
michael@0 64 0.224292f, 0.311806f,-0.102631f,-0.351851f,-0.034654f, 0.338330f, 0.166664f,-0.273300f,
michael@0 65 };
michael@0 66
michael@0 67 static const float analysis_window[240] = {
michael@0 68 0.000043f, 0.000171f, 0.000385f, 0.000685f, 0.001071f, 0.001541f, 0.002098f, 0.002739f,
michael@0 69 0.003466f, 0.004278f, 0.005174f, 0.006156f, 0.007222f, 0.008373f, 0.009607f, 0.010926f,
michael@0 70 0.012329f, 0.013815f, 0.015385f, 0.017037f, 0.018772f, 0.020590f, 0.022490f, 0.024472f,
michael@0 71 0.026535f, 0.028679f, 0.030904f, 0.033210f, 0.035595f, 0.038060f, 0.040604f, 0.043227f,
michael@0 72 0.045928f, 0.048707f, 0.051564f, 0.054497f, 0.057506f, 0.060591f, 0.063752f, 0.066987f,
michael@0 73 0.070297f, 0.073680f, 0.077136f, 0.080665f, 0.084265f, 0.087937f, 0.091679f, 0.095492f,
michael@0 74 0.099373f, 0.103323f, 0.107342f, 0.111427f, 0.115579f, 0.119797f, 0.124080f, 0.128428f,
michael@0 75 0.132839f, 0.137313f, 0.141849f, 0.146447f, 0.151105f, 0.155823f, 0.160600f, 0.165435f,
michael@0 76 0.170327f, 0.175276f, 0.180280f, 0.185340f, 0.190453f, 0.195619f, 0.200838f, 0.206107f,
michael@0 77 0.211427f, 0.216797f, 0.222215f, 0.227680f, 0.233193f, 0.238751f, 0.244353f, 0.250000f,
michael@0 78 0.255689f, 0.261421f, 0.267193f, 0.273005f, 0.278856f, 0.284744f, 0.290670f, 0.296632f,
michael@0 79 0.302628f, 0.308658f, 0.314721f, 0.320816f, 0.326941f, 0.333097f, 0.339280f, 0.345492f,
michael@0 80 0.351729f, 0.357992f, 0.364280f, 0.370590f, 0.376923f, 0.383277f, 0.389651f, 0.396044f,
michael@0 81 0.402455f, 0.408882f, 0.415325f, 0.421783f, 0.428254f, 0.434737f, 0.441231f, 0.447736f,
michael@0 82 0.454249f, 0.460770f, 0.467298f, 0.473832f, 0.480370f, 0.486912f, 0.493455f, 0.500000f,
michael@0 83 0.506545f, 0.513088f, 0.519630f, 0.526168f, 0.532702f, 0.539230f, 0.545751f, 0.552264f,
michael@0 84 0.558769f, 0.565263f, 0.571746f, 0.578217f, 0.584675f, 0.591118f, 0.597545f, 0.603956f,
michael@0 85 0.610349f, 0.616723f, 0.623077f, 0.629410f, 0.635720f, 0.642008f, 0.648271f, 0.654508f,
michael@0 86 0.660720f, 0.666903f, 0.673059f, 0.679184f, 0.685279f, 0.691342f, 0.697372f, 0.703368f,
michael@0 87 0.709330f, 0.715256f, 0.721144f, 0.726995f, 0.732807f, 0.738579f, 0.744311f, 0.750000f,
michael@0 88 0.755647f, 0.761249f, 0.766807f, 0.772320f, 0.777785f, 0.783203f, 0.788573f, 0.793893f,
michael@0 89 0.799162f, 0.804381f, 0.809547f, 0.814660f, 0.819720f, 0.824724f, 0.829673f, 0.834565f,
michael@0 90 0.839400f, 0.844177f, 0.848895f, 0.853553f, 0.858151f, 0.862687f, 0.867161f, 0.871572f,
michael@0 91 0.875920f, 0.880203f, 0.884421f, 0.888573f, 0.892658f, 0.896677f, 0.900627f, 0.904508f,
michael@0 92 0.908321f, 0.912063f, 0.915735f, 0.919335f, 0.922864f, 0.926320f, 0.929703f, 0.933013f,
michael@0 93 0.936248f, 0.939409f, 0.942494f, 0.945503f, 0.948436f, 0.951293f, 0.954072f, 0.956773f,
michael@0 94 0.959396f, 0.961940f, 0.964405f, 0.966790f, 0.969096f, 0.971321f, 0.973465f, 0.975528f,
michael@0 95 0.977510f, 0.979410f, 0.981228f, 0.982963f, 0.984615f, 0.986185f, 0.987671f, 0.989074f,
michael@0 96 0.990393f, 0.991627f, 0.992778f, 0.993844f, 0.994826f, 0.995722f, 0.996534f, 0.997261f,
michael@0 97 0.997902f, 0.998459f, 0.998929f, 0.999315f, 0.999615f, 0.999829f, 0.999957f, 1.000000f,
michael@0 98 };
michael@0 99
michael@0 100 static const int tbands[NB_TBANDS+1] = {
michael@0 101 2, 4, 6, 8, 10, 12, 14, 16, 20, 24, 28, 32, 40, 48, 56, 68, 80, 96, 120
michael@0 102 };
michael@0 103
michael@0 104 static const int extra_bands[NB_TOT_BANDS+1] = {
michael@0 105 1, 2, 4, 6, 8, 10, 12, 14, 16, 20, 24, 28, 32, 40, 48, 56, 68, 80, 96, 120, 160, 200
michael@0 106 };
michael@0 107
michael@0 108 /*static const float tweight[NB_TBANDS+1] = {
michael@0 109 .3, .4, .5, .6, .7, .8, .9, 1., 1., 1., 1., 1., 1., 1., .8, .7, .6, .5
michael@0 110 };*/
michael@0 111
michael@0 112 #define NB_TONAL_SKIP_BANDS 9
michael@0 113
michael@0 114 #define cA 0.43157974f
michael@0 115 #define cB 0.67848403f
michael@0 116 #define cC 0.08595542f
michael@0 117 #define cE ((float)M_PI/2)
michael@0 118 static OPUS_INLINE float fast_atan2f(float y, float x) {
michael@0 119 float x2, y2;
michael@0 120 /* Should avoid underflow on the values we'll get */
michael@0 121 if (ABS16(x)+ABS16(y)<1e-9f)
michael@0 122 {
michael@0 123 x*=1e12f;
michael@0 124 y*=1e12f;
michael@0 125 }
michael@0 126 x2 = x*x;
michael@0 127 y2 = y*y;
michael@0 128 if(x2<y2){
michael@0 129 float den = (y2 + cB*x2) * (y2 + cC*x2);
michael@0 130 if (den!=0)
michael@0 131 return -x*y*(y2 + cA*x2) / den + (y<0 ? -cE : cE);
michael@0 132 else
michael@0 133 return (y<0 ? -cE : cE);
michael@0 134 }else{
michael@0 135 float den = (x2 + cB*y2) * (x2 + cC*y2);
michael@0 136 if (den!=0)
michael@0 137 return x*y*(x2 + cA*y2) / den + (y<0 ? -cE : cE) - (x*y<0 ? -cE : cE);
michael@0 138 else
michael@0 139 return (y<0 ? -cE : cE) - (x*y<0 ? -cE : cE);
michael@0 140 }
michael@0 141 }
michael@0 142
michael@0 143 void tonality_get_info(TonalityAnalysisState *tonal, AnalysisInfo *info_out, int len)
michael@0 144 {
michael@0 145 int pos;
michael@0 146 int curr_lookahead;
michael@0 147 float psum;
michael@0 148 int i;
michael@0 149
michael@0 150 pos = tonal->read_pos;
michael@0 151 curr_lookahead = tonal->write_pos-tonal->read_pos;
michael@0 152 if (curr_lookahead<0)
michael@0 153 curr_lookahead += DETECT_SIZE;
michael@0 154
michael@0 155 if (len > 480 && pos != tonal->write_pos)
michael@0 156 {
michael@0 157 pos++;
michael@0 158 if (pos==DETECT_SIZE)
michael@0 159 pos=0;
michael@0 160 }
michael@0 161 if (pos == tonal->write_pos)
michael@0 162 pos--;
michael@0 163 if (pos<0)
michael@0 164 pos = DETECT_SIZE-1;
michael@0 165 OPUS_COPY(info_out, &tonal->info[pos], 1);
michael@0 166 tonal->read_subframe += len/120;
michael@0 167 while (tonal->read_subframe>=4)
michael@0 168 {
michael@0 169 tonal->read_subframe -= 4;
michael@0 170 tonal->read_pos++;
michael@0 171 }
michael@0 172 if (tonal->read_pos>=DETECT_SIZE)
michael@0 173 tonal->read_pos-=DETECT_SIZE;
michael@0 174
michael@0 175 /* Compensate for the delay in the features themselves.
michael@0 176 FIXME: Need a better estimate the 10 I just made up */
michael@0 177 curr_lookahead = IMAX(curr_lookahead-10, 0);
michael@0 178
michael@0 179 psum=0;
michael@0 180 /* Summing the probability of transition patterns that involve music at
michael@0 181 time (DETECT_SIZE-curr_lookahead-1) */
michael@0 182 for (i=0;i<DETECT_SIZE-curr_lookahead;i++)
michael@0 183 psum += tonal->pmusic[i];
michael@0 184 for (;i<DETECT_SIZE;i++)
michael@0 185 psum += tonal->pspeech[i];
michael@0 186 psum = psum*tonal->music_confidence + (1-psum)*tonal->speech_confidence;
michael@0 187 /*printf("%f %f %f\n", psum, info_out->music_prob, info_out->tonality);*/
michael@0 188
michael@0 189 info_out->music_prob = psum;
michael@0 190 }
michael@0 191
michael@0 192 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)
michael@0 193 {
michael@0 194 int i, b;
michael@0 195 const kiss_fft_state *kfft;
michael@0 196 VARDECL(kiss_fft_cpx, in);
michael@0 197 VARDECL(kiss_fft_cpx, out);
michael@0 198 int N = 480, N2=240;
michael@0 199 float * OPUS_RESTRICT A = tonal->angle;
michael@0 200 float * OPUS_RESTRICT dA = tonal->d_angle;
michael@0 201 float * OPUS_RESTRICT d2A = tonal->d2_angle;
michael@0 202 VARDECL(float, tonality);
michael@0 203 VARDECL(float, noisiness);
michael@0 204 float band_tonality[NB_TBANDS];
michael@0 205 float logE[NB_TBANDS];
michael@0 206 float BFCC[8];
michael@0 207 float features[25];
michael@0 208 float frame_tonality;
michael@0 209 float max_frame_tonality;
michael@0 210 /*float tw_sum=0;*/
michael@0 211 float frame_noisiness;
michael@0 212 const float pi4 = (float)(M_PI*M_PI*M_PI*M_PI);
michael@0 213 float slope=0;
michael@0 214 float frame_stationarity;
michael@0 215 float relativeE;
michael@0 216 float frame_probs[2];
michael@0 217 float alpha, alphaE, alphaE2;
michael@0 218 float frame_loudness;
michael@0 219 float bandwidth_mask;
michael@0 220 int bandwidth=0;
michael@0 221 float maxE = 0;
michael@0 222 float noise_floor;
michael@0 223 int remaining;
michael@0 224 AnalysisInfo *info;
michael@0 225 SAVE_STACK;
michael@0 226
michael@0 227 tonal->last_transition++;
michael@0 228 alpha = 1.f/IMIN(20, 1+tonal->count);
michael@0 229 alphaE = 1.f/IMIN(50, 1+tonal->count);
michael@0 230 alphaE2 = 1.f/IMIN(1000, 1+tonal->count);
michael@0 231
michael@0 232 if (tonal->count<4)
michael@0 233 tonal->music_prob = .5;
michael@0 234 kfft = celt_mode->mdct.kfft[0];
michael@0 235 if (tonal->count==0)
michael@0 236 tonal->mem_fill = 240;
michael@0 237 downmix(x, &tonal->inmem[tonal->mem_fill], IMIN(len, ANALYSIS_BUF_SIZE-tonal->mem_fill), offset, c1, c2, C);
michael@0 238 if (tonal->mem_fill+len < ANALYSIS_BUF_SIZE)
michael@0 239 {
michael@0 240 tonal->mem_fill += len;
michael@0 241 /* Don't have enough to update the analysis */
michael@0 242 RESTORE_STACK;
michael@0 243 return;
michael@0 244 }
michael@0 245 info = &tonal->info[tonal->write_pos++];
michael@0 246 if (tonal->write_pos>=DETECT_SIZE)
michael@0 247 tonal->write_pos-=DETECT_SIZE;
michael@0 248
michael@0 249 ALLOC(in, 480, kiss_fft_cpx);
michael@0 250 ALLOC(out, 480, kiss_fft_cpx);
michael@0 251 ALLOC(tonality, 240, float);
michael@0 252 ALLOC(noisiness, 240, float);
michael@0 253 for (i=0;i<N2;i++)
michael@0 254 {
michael@0 255 float w = analysis_window[i];
michael@0 256 in[i].r = (kiss_fft_scalar)(w*tonal->inmem[i]);
michael@0 257 in[i].i = (kiss_fft_scalar)(w*tonal->inmem[N2+i]);
michael@0 258 in[N-i-1].r = (kiss_fft_scalar)(w*tonal->inmem[N-i-1]);
michael@0 259 in[N-i-1].i = (kiss_fft_scalar)(w*tonal->inmem[N+N2-i-1]);
michael@0 260 }
michael@0 261 OPUS_MOVE(tonal->inmem, tonal->inmem+ANALYSIS_BUF_SIZE-240, 240);
michael@0 262 remaining = len - (ANALYSIS_BUF_SIZE-tonal->mem_fill);
michael@0 263 downmix(x, &tonal->inmem[240], remaining, offset+ANALYSIS_BUF_SIZE-tonal->mem_fill, c1, c2, C);
michael@0 264 tonal->mem_fill = 240 + remaining;
michael@0 265 opus_fft(kfft, in, out);
michael@0 266
michael@0 267 for (i=1;i<N2;i++)
michael@0 268 {
michael@0 269 float X1r, X2r, X1i, X2i;
michael@0 270 float angle, d_angle, d2_angle;
michael@0 271 float angle2, d_angle2, d2_angle2;
michael@0 272 float mod1, mod2, avg_mod;
michael@0 273 X1r = (float)out[i].r+out[N-i].r;
michael@0 274 X1i = (float)out[i].i-out[N-i].i;
michael@0 275 X2r = (float)out[i].i+out[N-i].i;
michael@0 276 X2i = (float)out[N-i].r-out[i].r;
michael@0 277
michael@0 278 angle = (float)(.5f/M_PI)*fast_atan2f(X1i, X1r);
michael@0 279 d_angle = angle - A[i];
michael@0 280 d2_angle = d_angle - dA[i];
michael@0 281
michael@0 282 angle2 = (float)(.5f/M_PI)*fast_atan2f(X2i, X2r);
michael@0 283 d_angle2 = angle2 - angle;
michael@0 284 d2_angle2 = d_angle2 - d_angle;
michael@0 285
michael@0 286 mod1 = d2_angle - (float)floor(.5+d2_angle);
michael@0 287 noisiness[i] = ABS16(mod1);
michael@0 288 mod1 *= mod1;
michael@0 289 mod1 *= mod1;
michael@0 290
michael@0 291 mod2 = d2_angle2 - (float)floor(.5+d2_angle2);
michael@0 292 noisiness[i] += ABS16(mod2);
michael@0 293 mod2 *= mod2;
michael@0 294 mod2 *= mod2;
michael@0 295
michael@0 296 avg_mod = .25f*(d2A[i]+2.f*mod1+mod2);
michael@0 297 tonality[i] = 1.f/(1.f+40.f*16.f*pi4*avg_mod)-.015f;
michael@0 298
michael@0 299 A[i] = angle2;
michael@0 300 dA[i] = d_angle2;
michael@0 301 d2A[i] = mod2;
michael@0 302 }
michael@0 303
michael@0 304 frame_tonality = 0;
michael@0 305 max_frame_tonality = 0;
michael@0 306 /*tw_sum = 0;*/
michael@0 307 info->activity = 0;
michael@0 308 frame_noisiness = 0;
michael@0 309 frame_stationarity = 0;
michael@0 310 if (!tonal->count)
michael@0 311 {
michael@0 312 for (b=0;b<NB_TBANDS;b++)
michael@0 313 {
michael@0 314 tonal->lowE[b] = 1e10;
michael@0 315 tonal->highE[b] = -1e10;
michael@0 316 }
michael@0 317 }
michael@0 318 relativeE = 0;
michael@0 319 frame_loudness = 0;
michael@0 320 for (b=0;b<NB_TBANDS;b++)
michael@0 321 {
michael@0 322 float E=0, tE=0, nE=0;
michael@0 323 float L1, L2;
michael@0 324 float stationarity;
michael@0 325 for (i=tbands[b];i<tbands[b+1];i++)
michael@0 326 {
michael@0 327 float binE = out[i].r*(float)out[i].r + out[N-i].r*(float)out[N-i].r
michael@0 328 + out[i].i*(float)out[i].i + out[N-i].i*(float)out[N-i].i;
michael@0 329 #ifdef FIXED_POINT
michael@0 330 /* FIXME: It's probably best to change the BFCC filter initial state instead */
michael@0 331 binE *= 5.55e-17f;
michael@0 332 #endif
michael@0 333 E += binE;
michael@0 334 tE += binE*tonality[i];
michael@0 335 nE += binE*2.f*(.5f-noisiness[i]);
michael@0 336 }
michael@0 337 tonal->E[tonal->E_count][b] = E;
michael@0 338 frame_noisiness += nE/(1e-15f+E);
michael@0 339
michael@0 340 frame_loudness += (float)sqrt(E+1e-10f);
michael@0 341 logE[b] = (float)log(E+1e-10f);
michael@0 342 tonal->lowE[b] = MIN32(logE[b], tonal->lowE[b]+.01f);
michael@0 343 tonal->highE[b] = MAX32(logE[b], tonal->highE[b]-.1f);
michael@0 344 if (tonal->highE[b] < tonal->lowE[b]+1.f)
michael@0 345 {
michael@0 346 tonal->highE[b]+=.5f;
michael@0 347 tonal->lowE[b]-=.5f;
michael@0 348 }
michael@0 349 relativeE += (logE[b]-tonal->lowE[b])/(1e-15f+tonal->highE[b]-tonal->lowE[b]);
michael@0 350
michael@0 351 L1=L2=0;
michael@0 352 for (i=0;i<NB_FRAMES;i++)
michael@0 353 {
michael@0 354 L1 += (float)sqrt(tonal->E[i][b]);
michael@0 355 L2 += tonal->E[i][b];
michael@0 356 }
michael@0 357
michael@0 358 stationarity = MIN16(0.99f,L1/(float)sqrt(1e-15+NB_FRAMES*L2));
michael@0 359 stationarity *= stationarity;
michael@0 360 stationarity *= stationarity;
michael@0 361 frame_stationarity += stationarity;
michael@0 362 /*band_tonality[b] = tE/(1e-15+E)*/;
michael@0 363 band_tonality[b] = MAX16(tE/(1e-15f+E), stationarity*tonal->prev_band_tonality[b]);
michael@0 364 #if 0
michael@0 365 if (b>=NB_TONAL_SKIP_BANDS)
michael@0 366 {
michael@0 367 frame_tonality += tweight[b]*band_tonality[b];
michael@0 368 tw_sum += tweight[b];
michael@0 369 }
michael@0 370 #else
michael@0 371 frame_tonality += band_tonality[b];
michael@0 372 if (b>=NB_TBANDS-NB_TONAL_SKIP_BANDS)
michael@0 373 frame_tonality -= band_tonality[b-NB_TBANDS+NB_TONAL_SKIP_BANDS];
michael@0 374 #endif
michael@0 375 max_frame_tonality = MAX16(max_frame_tonality, (1.f+.03f*(b-NB_TBANDS))*frame_tonality);
michael@0 376 slope += band_tonality[b]*(b-8);
michael@0 377 /*printf("%f %f ", band_tonality[b], stationarity);*/
michael@0 378 tonal->prev_band_tonality[b] = band_tonality[b];
michael@0 379 }
michael@0 380
michael@0 381 bandwidth_mask = 0;
michael@0 382 bandwidth = 0;
michael@0 383 maxE = 0;
michael@0 384 noise_floor = 5.7e-4f/(1<<(IMAX(0,lsb_depth-8)));
michael@0 385 #ifdef FIXED_POINT
michael@0 386 noise_floor *= 1<<(15+SIG_SHIFT);
michael@0 387 #endif
michael@0 388 noise_floor *= noise_floor;
michael@0 389 for (b=0;b<NB_TOT_BANDS;b++)
michael@0 390 {
michael@0 391 float E=0;
michael@0 392 int band_start, band_end;
michael@0 393 /* Keep a margin of 300 Hz for aliasing */
michael@0 394 band_start = extra_bands[b];
michael@0 395 band_end = extra_bands[b+1];
michael@0 396 for (i=band_start;i<band_end;i++)
michael@0 397 {
michael@0 398 float binE = out[i].r*(float)out[i].r + out[N-i].r*(float)out[N-i].r
michael@0 399 + out[i].i*(float)out[i].i + out[N-i].i*(float)out[N-i].i;
michael@0 400 E += binE;
michael@0 401 }
michael@0 402 maxE = MAX32(maxE, E);
michael@0 403 tonal->meanE[b] = MAX32((1-alphaE2)*tonal->meanE[b], E);
michael@0 404 E = MAX32(E, tonal->meanE[b]);
michael@0 405 /* Use a simple follower with 13 dB/Bark slope for spreading function */
michael@0 406 bandwidth_mask = MAX32(.05f*bandwidth_mask, E);
michael@0 407 /* Consider the band "active" only if all these conditions are met:
michael@0 408 1) less than 10 dB below the simple follower
michael@0 409 2) less than 90 dB below the peak band (maximal masking possible considering
michael@0 410 both the ATH and the loudness-dependent slope of the spreading function)
michael@0 411 3) above the PCM quantization noise floor
michael@0 412 */
michael@0 413 if (E>.1*bandwidth_mask && E*1e9f > maxE && E > noise_floor*(band_end-band_start))
michael@0 414 bandwidth = b;
michael@0 415 }
michael@0 416 if (tonal->count<=2)
michael@0 417 bandwidth = 20;
michael@0 418 frame_loudness = 20*(float)log10(frame_loudness);
michael@0 419 tonal->Etracker = MAX32(tonal->Etracker-.03f, frame_loudness);
michael@0 420 tonal->lowECount *= (1-alphaE);
michael@0 421 if (frame_loudness < tonal->Etracker-30)
michael@0 422 tonal->lowECount += alphaE;
michael@0 423
michael@0 424 for (i=0;i<8;i++)
michael@0 425 {
michael@0 426 float sum=0;
michael@0 427 for (b=0;b<16;b++)
michael@0 428 sum += dct_table[i*16+b]*logE[b];
michael@0 429 BFCC[i] = sum;
michael@0 430 }
michael@0 431
michael@0 432 frame_stationarity /= NB_TBANDS;
michael@0 433 relativeE /= NB_TBANDS;
michael@0 434 if (tonal->count<10)
michael@0 435 relativeE = .5;
michael@0 436 frame_noisiness /= NB_TBANDS;
michael@0 437 #if 1
michael@0 438 info->activity = frame_noisiness + (1-frame_noisiness)*relativeE;
michael@0 439 #else
michael@0 440 info->activity = .5*(1+frame_noisiness-frame_stationarity);
michael@0 441 #endif
michael@0 442 frame_tonality = (max_frame_tonality/(NB_TBANDS-NB_TONAL_SKIP_BANDS));
michael@0 443 frame_tonality = MAX16(frame_tonality, tonal->prev_tonality*.8f);
michael@0 444 tonal->prev_tonality = frame_tonality;
michael@0 445
michael@0 446 slope /= 8*8;
michael@0 447 info->tonality_slope = slope;
michael@0 448
michael@0 449 tonal->E_count = (tonal->E_count+1)%NB_FRAMES;
michael@0 450 tonal->count++;
michael@0 451 info->tonality = frame_tonality;
michael@0 452
michael@0 453 for (i=0;i<4;i++)
michael@0 454 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];
michael@0 455
michael@0 456 for (i=0;i<4;i++)
michael@0 457 tonal->cmean[i] = (1-alpha)*tonal->cmean[i] + alpha*BFCC[i];
michael@0 458
michael@0 459 for (i=0;i<4;i++)
michael@0 460 features[4+i] = 0.63246f*(BFCC[i]-tonal->mem[i+24]) + 0.31623f*(tonal->mem[i]-tonal->mem[i+16]);
michael@0 461 for (i=0;i<3;i++)
michael@0 462 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];
michael@0 463
michael@0 464 if (tonal->count > 5)
michael@0 465 {
michael@0 466 for (i=0;i<9;i++)
michael@0 467 tonal->std[i] = (1-alpha)*tonal->std[i] + alpha*features[i]*features[i];
michael@0 468 }
michael@0 469
michael@0 470 for (i=0;i<8;i++)
michael@0 471 {
michael@0 472 tonal->mem[i+24] = tonal->mem[i+16];
michael@0 473 tonal->mem[i+16] = tonal->mem[i+8];
michael@0 474 tonal->mem[i+8] = tonal->mem[i];
michael@0 475 tonal->mem[i] = BFCC[i];
michael@0 476 }
michael@0 477 for (i=0;i<9;i++)
michael@0 478 features[11+i] = (float)sqrt(tonal->std[i]);
michael@0 479 features[20] = info->tonality;
michael@0 480 features[21] = info->activity;
michael@0 481 features[22] = frame_stationarity;
michael@0 482 features[23] = info->tonality_slope;
michael@0 483 features[24] = tonal->lowECount;
michael@0 484
michael@0 485 #ifndef DISABLE_FLOAT_API
michael@0 486 mlp_process(&net, features, frame_probs);
michael@0 487 frame_probs[0] = .5f*(frame_probs[0]+1);
michael@0 488 /* Curve fitting between the MLP probability and the actual probability */
michael@0 489 frame_probs[0] = .01f + 1.21f*frame_probs[0]*frame_probs[0] - .23f*(float)pow(frame_probs[0], 10);
michael@0 490 /* Probability of active audio (as opposed to silence) */
michael@0 491 frame_probs[1] = .5f*frame_probs[1]+.5f;
michael@0 492 /* Consider that silence has a 50-50 probability. */
michael@0 493 frame_probs[0] = frame_probs[1]*frame_probs[0] + (1-frame_probs[1])*.5f;
michael@0 494
michael@0 495 /*printf("%f %f ", frame_probs[0], frame_probs[1]);*/
michael@0 496 {
michael@0 497 /* Probability of state transition */
michael@0 498 float tau;
michael@0 499 /* Represents independence of the MLP probabilities, where
michael@0 500 beta=1 means fully independent. */
michael@0 501 float beta;
michael@0 502 /* Denormalized probability of speech (p0) and music (p1) after update */
michael@0 503 float p0, p1;
michael@0 504 /* Probabilities for "all speech" and "all music" */
michael@0 505 float s0, m0;
michael@0 506 /* Probability sum for renormalisation */
michael@0 507 float psum;
michael@0 508 /* Instantaneous probability of speech and music, with beta pre-applied. */
michael@0 509 float speech0;
michael@0 510 float music0;
michael@0 511
michael@0 512 /* One transition every 3 minutes of active audio */
michael@0 513 tau = .00005f*frame_probs[1];
michael@0 514 beta = .05f;
michael@0 515 if (1) {
michael@0 516 /* Adapt beta based on how "unexpected" the new prob is */
michael@0 517 float p, q;
michael@0 518 p = MAX16(.05f,MIN16(.95f,frame_probs[0]));
michael@0 519 q = MAX16(.05f,MIN16(.95f,tonal->music_prob));
michael@0 520 beta = .01f+.05f*ABS16(p-q)/(p*(1-q)+q*(1-p));
michael@0 521 }
michael@0 522 /* p0 and p1 are the probabilities of speech and music at this frame
michael@0 523 using only information from previous frame and applying the
michael@0 524 state transition model */
michael@0 525 p0 = (1-tonal->music_prob)*(1-tau) + tonal->music_prob *tau;
michael@0 526 p1 = tonal->music_prob *(1-tau) + (1-tonal->music_prob)*tau;
michael@0 527 /* We apply the current probability with exponent beta to work around
michael@0 528 the fact that the probability estimates aren't independent. */
michael@0 529 p0 *= (float)pow(1-frame_probs[0], beta);
michael@0 530 p1 *= (float)pow(frame_probs[0], beta);
michael@0 531 /* Normalise the probabilities to get the Marokv probability of music. */
michael@0 532 tonal->music_prob = p1/(p0+p1);
michael@0 533 info->music_prob = tonal->music_prob;
michael@0 534
michael@0 535 /* This chunk of code deals with delayed decision. */
michael@0 536 psum=1e-20f;
michael@0 537 /* Instantaneous probability of speech and music, with beta pre-applied. */
michael@0 538 speech0 = (float)pow(1-frame_probs[0], beta);
michael@0 539 music0 = (float)pow(frame_probs[0], beta);
michael@0 540 if (tonal->count==1)
michael@0 541 {
michael@0 542 tonal->pspeech[0]=.5;
michael@0 543 tonal->pmusic [0]=.5;
michael@0 544 }
michael@0 545 /* Updated probability of having only speech (s0) or only music (m0),
michael@0 546 before considering the new observation. */
michael@0 547 s0 = tonal->pspeech[0] + tonal->pspeech[1];
michael@0 548 m0 = tonal->pmusic [0] + tonal->pmusic [1];
michael@0 549 /* Updates s0 and m0 with instantaneous probability. */
michael@0 550 tonal->pspeech[0] = s0*(1-tau)*speech0;
michael@0 551 tonal->pmusic [0] = m0*(1-tau)*music0;
michael@0 552 /* Propagate the transition probabilities */
michael@0 553 for (i=1;i<DETECT_SIZE-1;i++)
michael@0 554 {
michael@0 555 tonal->pspeech[i] = tonal->pspeech[i+1]*speech0;
michael@0 556 tonal->pmusic [i] = tonal->pmusic [i+1]*music0;
michael@0 557 }
michael@0 558 /* Probability that the latest frame is speech, when all the previous ones were music. */
michael@0 559 tonal->pspeech[DETECT_SIZE-1] = m0*tau*speech0;
michael@0 560 /* Probability that the latest frame is music, when all the previous ones were speech. */
michael@0 561 tonal->pmusic [DETECT_SIZE-1] = s0*tau*music0;
michael@0 562
michael@0 563 /* Renormalise probabilities to 1 */
michael@0 564 for (i=0;i<DETECT_SIZE;i++)
michael@0 565 psum += tonal->pspeech[i] + tonal->pmusic[i];
michael@0 566 psum = 1.f/psum;
michael@0 567 for (i=0;i<DETECT_SIZE;i++)
michael@0 568 {
michael@0 569 tonal->pspeech[i] *= psum;
michael@0 570 tonal->pmusic [i] *= psum;
michael@0 571 }
michael@0 572 psum = tonal->pmusic[0];
michael@0 573 for (i=1;i<DETECT_SIZE;i++)
michael@0 574 psum += tonal->pspeech[i];
michael@0 575
michael@0 576 /* Estimate our confidence in the speech/music decisions */
michael@0 577 if (frame_probs[1]>.75)
michael@0 578 {
michael@0 579 if (tonal->music_prob>.9)
michael@0 580 {
michael@0 581 float adapt;
michael@0 582 adapt = 1.f/(++tonal->music_confidence_count);
michael@0 583 tonal->music_confidence_count = IMIN(tonal->music_confidence_count, 500);
michael@0 584 tonal->music_confidence += adapt*MAX16(-.2f,frame_probs[0]-tonal->music_confidence);
michael@0 585 }
michael@0 586 if (tonal->music_prob<.1)
michael@0 587 {
michael@0 588 float adapt;
michael@0 589 adapt = 1.f/(++tonal->speech_confidence_count);
michael@0 590 tonal->speech_confidence_count = IMIN(tonal->speech_confidence_count, 500);
michael@0 591 tonal->speech_confidence += adapt*MIN16(.2f,frame_probs[0]-tonal->speech_confidence);
michael@0 592 }
michael@0 593 } else {
michael@0 594 if (tonal->music_confidence_count==0)
michael@0 595 tonal->music_confidence = .9f;
michael@0 596 if (tonal->speech_confidence_count==0)
michael@0 597 tonal->speech_confidence = .1f;
michael@0 598 }
michael@0 599 }
michael@0 600 if (tonal->last_music != (tonal->music_prob>.5f))
michael@0 601 tonal->last_transition=0;
michael@0 602 tonal->last_music = tonal->music_prob>.5f;
michael@0 603 #else
michael@0 604 info->music_prob = 0;
michael@0 605 #endif
michael@0 606 /*for (i=0;i<25;i++)
michael@0 607 printf("%f ", features[i]);
michael@0 608 printf("\n");*/
michael@0 609
michael@0 610 info->bandwidth = bandwidth;
michael@0 611 /*printf("%d %d\n", info->bandwidth, info->opus_bandwidth);*/
michael@0 612 info->noisiness = frame_noisiness;
michael@0 613 info->valid = 1;
michael@0 614 if (info_out!=NULL)
michael@0 615 OPUS_COPY(info_out, info, 1);
michael@0 616 RESTORE_STACK;
michael@0 617 }
michael@0 618
michael@0 619 void run_analysis(TonalityAnalysisState *analysis, const CELTMode *celt_mode, const void *analysis_pcm,
michael@0 620 int analysis_frame_size, int frame_size, int c1, int c2, int C, opus_int32 Fs,
michael@0 621 int lsb_depth, downmix_func downmix, AnalysisInfo *analysis_info)
michael@0 622 {
michael@0 623 int offset;
michael@0 624 int pcm_len;
michael@0 625
michael@0 626 if (analysis_pcm != NULL)
michael@0 627 {
michael@0 628 /* Avoid overflow/wrap-around of the analysis buffer */
michael@0 629 analysis_frame_size = IMIN((DETECT_SIZE-5)*Fs/100, analysis_frame_size);
michael@0 630
michael@0 631 pcm_len = analysis_frame_size - analysis->analysis_offset;
michael@0 632 offset = analysis->analysis_offset;
michael@0 633 do {
michael@0 634 tonality_analysis(analysis, NULL, celt_mode, analysis_pcm, IMIN(480, pcm_len), offset, c1, c2, C, lsb_depth, downmix);
michael@0 635 offset += 480;
michael@0 636 pcm_len -= 480;
michael@0 637 } while (pcm_len>0);
michael@0 638 analysis->analysis_offset = analysis_frame_size;
michael@0 639
michael@0 640 analysis->analysis_offset -= frame_size;
michael@0 641 }
michael@0 642
michael@0 643 analysis_info->valid = 0;
michael@0 644 tonality_get_info(analysis, analysis_info, frame_size);
michael@0 645 }

mercurial