Adds automatic bandwidth detection
[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 #include "kiss_fft.h"
33 #include "celt.h"
34 #include "modes.h"
35 #include "arch.h"
36 #include "quant_bands.h"
37 #include <stdio.h>
38 #ifndef FIXED_POINT
39 #include "mlp.c"
40 #include "mlp_data.c"
41 #endif
42
43 #ifndef M_PI
44 #define M_PI 3.141592653
45 #endif
46
47 float dct_table[128] = {
48         0.250000, 0.250000, 0.250000, 0.250000, 0.250000, 0.250000, 0.250000, 0.250000,
49         0.250000, 0.250000, 0.250000, 0.250000, 0.250000, 0.250000, 0.250000, 0.250000,
50         0.351851, 0.338330, 0.311806, 0.273300, 0.224292, 0.166664, 0.102631, 0.034654,
51         -0.034654, -0.102631, -0.166664, -0.224292, -0.273300, -0.311806, -0.338330, -0.351851,
52         0.346760, 0.293969, 0.196424, 0.068975, -0.068975, -0.196424, -0.293969, -0.346760,
53         -0.346760, -0.293969, -0.196424, -0.068975, 0.068975, 0.196424, 0.293969, 0.346760,
54         0.338330, 0.224292, 0.034654, -0.166664, -0.311806, -0.351851, -0.273300, -0.102631,
55         0.102631, 0.273300, 0.351851, 0.311806, 0.166664, -0.034654, -0.224292, -0.338330,
56         0.326641, 0.135299, -0.135299, -0.326641, -0.326641, -0.135299, 0.135299, 0.326641,
57         0.326641, 0.135299, -0.135299, -0.326641, -0.326641, -0.135299, 0.135299, 0.326641,
58         0.311806, 0.034654, -0.273300, -0.338330, -0.102631, 0.224292, 0.351851, 0.166664,
59         -0.166664, -0.351851, -0.224292, 0.102631, 0.338330, 0.273300, -0.034654, -0.311806,
60         0.293969, -0.068975, -0.346760, -0.196424, 0.196424, 0.346760, 0.068975, -0.293969,
61         -0.293969, 0.068975, 0.346760, 0.196424, -0.196424, -0.346760, -0.068975, 0.293969,
62         0.273300, -0.166664, -0.338330, 0.034654, 0.351851, 0.102631, -0.311806, -0.224292,
63         0.224292, 0.311806, -0.102631, -0.351851, -0.034654, 0.338330, 0.166664, -0.273300,
64 };
65
66 #define NB_FRAMES 8
67
68 #define NB_TBANDS 18
69 static const int tbands[NB_TBANDS+1] = {
70       2, 4, 6, 8, 10, 12, 14, 16, 20, 24, 28, 32, 40, 48, 56, 68, 80, 96, 120
71 };
72
73 #define NB_TONAL_SKIP_BANDS 8
74
75 typedef struct {
76    float angle[240];
77    float d_angle[240];
78    float d2_angle[240];
79    float prev_band_tonality[NB_TBANDS];
80    float prev_tonality;
81    float E[NB_FRAMES][NB_TBANDS];
82    float lowE[NB_TBANDS], highE[NB_TBANDS];
83    float meanE[NB_TBANDS], meanRE[NB_TBANDS];
84    float mem[32];
85    float cmean[8];
86    float std[9];
87    float music_prob;
88    float Etracker;
89    float lowECount;
90    int E_count;
91    int last_music;
92    int last_transition;
93    int count;
94    int opus_bandwidth;
95 } TonalityAnalysisState;
96
97 void tonality_analysis(TonalityAnalysisState *tonal, AnalysisInfo *info, CELTEncoder *celt_enc, const opus_val16 *x, int C)
98 {
99     int i, b;
100     const CELTMode *mode;
101     const kiss_fft_state *kfft;
102     kiss_fft_cpx in[480], out[480];
103     int N = 480, N2=240;
104     float * restrict A = tonal->angle;
105     float * restrict dA = tonal->d_angle;
106     float * restrict d2A = tonal->d2_angle;
107     float tonality[240];
108     float noisiness[240];
109     float band_tonality[NB_TBANDS];
110     float logE[NB_TBANDS];
111     float BFCC[8];
112     float features[100];
113     float frame_tonality;
114     float frame_noisiness;
115     const float pi4 = M_PI*M_PI*M_PI*M_PI;
116     float slope=0;
117     float frame_stationarity;
118     float relativeE;
119     float frame_prob;
120     float alpha, alphaE, alphaE2;
121     float frame_loudness;
122     float bandwidth_mask;
123     int bandwidth=0;
124     float bandE[NB_TBANDS];
125     celt_encoder_ctl(celt_enc, CELT_GET_MODE(&mode));
126
127     tonal->last_transition++;
128     alpha = 1.f/IMIN(20, 1+tonal->count);
129     alphaE = 1.f/IMIN(50, 1+tonal->count);
130     alphaE2 = 1.f/IMIN(6000, 1+tonal->count);
131
132     if (tonal->count<4)
133        tonal->music_prob = .5;
134     kfft = mode->mdct.kfft[0];
135     if (C==1)
136     {
137        for (i=0;i<N2;i++)
138        {
139           float w = .5-.5*cos(M_PI*(i+1)/N2);
140           in[i].r = MULT16_16(w, x[i]);
141           in[i].i = MULT16_16(w, x[N-N2+i]);
142           in[N-i-1].r = MULT16_16(w, x[N-i-1]);
143           in[N-i-1].i = MULT16_16(w, x[2*N-N2-i-1]);
144        }
145     } else {
146        for (i=0;i<N2;i++)
147        {
148           float w = .5-.5*cos(M_PI*(i+1)/N2);
149           in[i].r = MULT16_16(w, x[2*i]+x[2*i+1]);
150           in[i].i = MULT16_16(w, x[2*(N-N2+i)]+x[2*(N-N2+i)+1]);
151           in[N-i-1].r = MULT16_16(w, x[2*(N-i-1)]+x[2*(N-i-1)+1]);
152           in[N-i-1].i = MULT16_16(w, x[2*(2*N-N2-i-1)]+x[2*(2*N-N2-i-1)+1]);
153        }
154     }
155     opus_fft(kfft, in, out);
156
157     for (i=1;i<N2;i++)
158     {
159        float X1r, X2r, X1i, X2i;
160        float angle, d_angle, d2_angle;
161        float angle2, d_angle2, d2_angle2;
162        float mod1, mod2, avg_mod;
163        X1r = out[i].r+out[N-i].r;
164        X1i = out[i].i-out[N-i].i;
165        X2r = out[i].i+out[N-i].i;
166        X2i = out[N-i].r-out[i].r;
167
168        angle = (.5/M_PI)*atan2(X1i, X1r);
169        d_angle = angle - A[i];
170        d2_angle = d_angle - dA[i];
171
172        angle2 = (.5/M_PI)*atan2(X2i, X2r);
173        d_angle2 = angle2 - angle;
174        d2_angle2 = d_angle2 - d_angle;
175
176        mod1 = d2_angle - floor(.5+d2_angle);
177        noisiness[i] = fabs(mod1);
178        mod1 *= mod1;
179        mod1 *= mod1;
180
181        mod2 = d2_angle2 - floor(.5+d2_angle2);
182        noisiness[i] += fabs(mod2);
183        mod2 *= mod2;
184        mod2 *= mod2;
185
186        avg_mod = .25*(d2A[i]+2*mod1+mod2);
187        tonality[i] = 1./(1+40*16*pi4*avg_mod)-.015;
188
189        A[i] = angle2;
190        dA[i] = d_angle2;
191        d2A[i] = mod2;
192     }
193
194     frame_tonality = 0;
195     info->activity = 0;
196     frame_noisiness = 0;
197     frame_stationarity = 0;
198     if (!tonal->count)
199     {
200        for (b=0;b<NB_TBANDS;b++)
201        {
202           tonal->lowE[b] = 1e10;
203           tonal->highE[b] = -1e10;
204        }
205     }
206     relativeE = 0;
207     info->boost_amount[0]=info->boost_amount[1]=0;
208     info->boost_band[0]=info->boost_band[1]=0;
209     frame_loudness = 0;
210     bandwidth_mask = 0;
211     for (b=0;b<NB_TBANDS;b++)
212     {
213        float E=0, tE=0, nE=0;
214        float L1, L2;
215        float stationarity;
216        for (i=tbands[b];i<tbands[b+1];i++)
217        {
218           float binE = out[i].r*out[i].r + out[N-i].r*out[N-i].r
219                      + out[i].i*out[i].i + out[N-i].i*out[N-i].i;
220           E += binE;
221           tE += binE*tonality[i];
222           nE += binE*2*(.5-noisiness[i]);
223        }
224        bandE[b] = E;
225        tonal->E[tonal->E_count][b] = E;
226        frame_noisiness += nE/(1e-15+E);
227
228        frame_loudness += sqrt(E+1e-10);
229        /* Add a reasonable noise floor */
230        tonal->meanE[b] = (1-alphaE2)*tonal->meanE[b] + alphaE2*E;
231        tonal->meanRE[b] = (1-alphaE2)*tonal->meanRE[b] + alphaE2*sqrt(E);
232        /* 13 dB slope for spreading function */
233        bandwidth_mask = MAX32(.05*bandwidth_mask, E);
234        /* Checks if band looks like stationary noise or if it's below a (trivial) masking curve */
235        if (tonal->meanRE[b]*tonal->meanRE[b] < tonal->meanE[b]*.95 && E>.1*bandwidth_mask)
236           bandwidth = b;
237        logE[b] = log(E+1e-10);
238        tonal->lowE[b] = MIN32(logE[b], tonal->lowE[b]+.01);
239        tonal->highE[b] = MAX32(logE[b], tonal->highE[b]-.1);
240        if (tonal->highE[b] < tonal->lowE[b]+1)
241        {
242           tonal->highE[b]+=.5;
243           tonal->lowE[b]-=.5;
244        }
245        relativeE += (logE[b]-tonal->lowE[b])/(EPSILON+tonal->highE[b]-tonal->lowE[b]);
246
247        L1=L2=0;
248        for (i=0;i<NB_FRAMES;i++)
249        {
250           L1 += sqrt(tonal->E[i][b]);
251           L2 += tonal->E[i][b];
252        }
253
254        stationarity = MIN16(0.99,L1/sqrt(EPSILON+NB_FRAMES*L2));
255        stationarity *= stationarity;
256        stationarity *= stationarity;
257        frame_stationarity += stationarity;
258        /*band_tonality[b] = tE/(1e-15+E)*/;
259        band_tonality[b] = MAX16(tE/(EPSILON+E), stationarity*tonal->prev_band_tonality[b]);
260        if (b>=NB_TONAL_SKIP_BANDS)
261           frame_tonality += band_tonality[b];
262        slope += band_tonality[b]*(b-8);
263        if (band_tonality[b] > info->boost_amount[1] && b>=7 && b < NB_TBANDS-1)
264        {
265           if (band_tonality[b] > info->boost_amount[0])
266           {
267              info->boost_amount[1] = info->boost_amount[0];
268              info->boost_band[1] = info->boost_band[0];
269              info->boost_amount[0] = band_tonality[b];
270              info->boost_band[0] = b;
271           } else {
272              info->boost_amount[1] = band_tonality[b];
273              info->boost_band[1] = b;
274           }
275        }
276        tonal->prev_band_tonality[b] = band_tonality[b];
277     }
278
279     frame_loudness = 20*log10(frame_loudness);
280     tonal->Etracker = MAX32(tonal->Etracker-.03, frame_loudness);
281     tonal->lowECount *= (1-alphaE);
282     if (frame_loudness < tonal->Etracker-30)
283        tonal->lowECount += alphaE;
284
285     for (i=0;i<8;i++)
286     {
287        float sum=0;
288        for (b=0;b<16;b++)
289           sum += dct_table[i*16+b]*logE[b];
290        BFCC[i] = sum;
291     }
292
293     frame_stationarity /= NB_TBANDS;
294     relativeE /= NB_TBANDS;
295     if (tonal->count<10)
296        relativeE = .5;
297     frame_noisiness /= NB_TBANDS;
298 #if 1
299     info->activity = frame_noisiness + (1-frame_noisiness)*relativeE;
300 #else
301     info->activity = .5*(1+frame_noisiness-frame_stationarity);
302 #endif
303     frame_tonality /= NB_TBANDS-NB_TONAL_SKIP_BANDS;
304     frame_tonality = MAX16(frame_tonality, tonal->prev_tonality*.8);
305     tonal->prev_tonality = frame_tonality;
306     info->boost_amount[0] -= frame_tonality+.2;
307     info->boost_amount[1] -= frame_tonality+.2;
308     if (band_tonality[info->boost_band[0]] < band_tonality[info->boost_band[0]+1]+.15
309         || band_tonality[info->boost_band[0]] < band_tonality[info->boost_band[0]-1]+.15)
310        info->boost_amount[0]=0;
311     if (band_tonality[info->boost_band[1]] < band_tonality[info->boost_band[1]+1]+.15
312         || band_tonality[info->boost_band[1]] < band_tonality[info->boost_band[1]-1]+.15)
313        info->boost_amount[1]=0;
314
315     slope /= 8*8;
316     info->tonality_slope = slope;
317
318     tonal->E_count = (tonal->E_count+1)%NB_FRAMES;
319     tonal->count++;
320     info->tonality = frame_tonality;
321
322     for (i=0;i<4;i++)
323        features[i] = -0.12299*(BFCC[i]+tonal->mem[i+24]) + 0.49195*(tonal->mem[i]+tonal->mem[i+16]) + 0.69693*tonal->mem[i+8] - 1.4349*tonal->cmean[i];
324
325     for (i=0;i<4;i++)
326        tonal->cmean[i] = (1-alpha)*tonal->cmean[i] + alpha*BFCC[i];
327
328     for (i=0;i<4;i++)
329         features[4+i] = 0.63246*(BFCC[i]-tonal->mem[i+24]) + 0.31623*(tonal->mem[i]-tonal->mem[i+16]);
330     for (i=0;i<3;i++)
331         features[8+i] = 0.53452*(BFCC[i]+tonal->mem[i+24]) - 0.26726*(tonal->mem[i]+tonal->mem[i+16]) -0.53452*tonal->mem[i+8];
332
333     if (tonal->count > 5)
334     {
335        for (i=0;i<9;i++)
336           tonal->std[i] = (1-alpha)*tonal->std[i] + alpha*features[i]*features[i];
337     }
338
339     for (i=0;i<8;i++)
340     {
341        tonal->mem[i+24] = tonal->mem[i+16];
342        tonal->mem[i+16] = tonal->mem[i+8];
343        tonal->mem[i+8] = tonal->mem[i];
344        tonal->mem[i] = BFCC[i];
345     }
346     for (i=0;i<9;i++)
347        features[11+i] = sqrt(tonal->std[i]);
348     features[20] = info->tonality;
349     features[21] = info->activity;
350     features[22] = frame_stationarity;
351     features[23] = info->tonality_slope;
352     features[24] = tonal->lowECount;
353
354 #ifndef FIXED_POINT
355     mlp_process(&net, features, &frame_prob);
356     /* Adds a "probability dead zone", with a cap on certainty */
357     frame_prob = .90*frame_prob*frame_prob*frame_prob;
358
359     frame_prob = .5*(frame_prob+1);
360
361     /*printf("%f\n", frame_prob);*/
362     {
363        float tau, beta;
364        float p0, p1;
365        float max_certainty;
366        /* One transition every 3 minutes */
367        tau = .00005;
368        beta = .1;
369        max_certainty = 1.f/(10+1*tonal->last_transition);
370        p0 = (1-tonal->music_prob)*(1-tau) +    tonal->music_prob *tau;
371        p1 =    tonal->music_prob *(1-tau) + (1-tonal->music_prob)*tau;
372        p0 *= pow(1-frame_prob, beta);
373        p1 *= pow(frame_prob, beta);
374        tonal->music_prob = MAX16(max_certainty, MIN16(1-max_certainty, p1/(p0+p1)));
375        info->music_prob = tonal->music_prob;
376        /*printf("%f %f\n", frame_prob, info->music_prob);*/
377     }
378     if (tonal->last_music != (tonal->music_prob>.5))
379        tonal->last_transition=0;
380     tonal->last_music = tonal->music_prob>.5;
381 #else
382     info->music_prob = 0;
383 #endif
384     /*for (i=0;i<25;i++)
385        printf("%f ", features[i]);
386     printf("\n");*/
387
388     /* FIXME: Can't detect SWB for now because the last band ends at 12 kHz */
389     if (bandwidth == NB_TBANDS-1 || tonal->count<100)
390     {
391        tonal->opus_bandwidth = OPUS_BANDWIDTH_FULLBAND;
392     } else {
393        int close_enough = 0;
394        if (bandE[bandwidth-1] < 3000*bandE[NB_TBANDS-1] && bandwidth < NB_TBANDS-1)
395           close_enough=1;
396        if (bandwidth<=11 || (bandwidth==12 && close_enough))
397           tonal->opus_bandwidth = OPUS_BANDWIDTH_NARROWBAND;
398        else if (bandwidth<=13)
399           tonal->opus_bandwidth = OPUS_BANDWIDTH_MEDIUMBAND;
400        else if (bandwidth<=15 || (bandwidth==16 && close_enough))
401           tonal->opus_bandwidth = OPUS_BANDWIDTH_WIDEBAND;
402     }
403     info->valid = 1;
404 }