8f2971e49222a39eabc7cb85ee1328a78ca24081
[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 mem[32];
84    float cmean[8];
85    float std[9];
86    float music_prob;
87    float Etracker;
88    float lowECount;
89    int E_count;
90    int last_music;
91    int last_transition;
92    int count;
93 } TonalityAnalysisState;
94
95 void tonality_analysis(TonalityAnalysisState *tonal, AnalysisInfo *info, CELTEncoder *celt_enc, const opus_val16 *x, int C)
96 {
97     int i, b;
98     const CELTMode *mode;
99     const kiss_fft_state *kfft;
100     kiss_fft_cpx in[480], out[480];
101     int N = 480, N2=240;
102     float * restrict A = tonal->angle;
103     float * restrict dA = tonal->d_angle;
104     float * restrict d2A = tonal->d2_angle;
105     float tonality[240];
106     float noisiness[240];
107     float band_tonality[NB_TBANDS];
108     float logE[NB_TBANDS];
109     float BFCC[8];
110     float features[100];
111     float frame_tonality;
112     float frame_noisiness;
113     const float pi4 = M_PI*M_PI*M_PI*M_PI;
114     float slope=0;
115     float frame_stationarity;
116     float relativeE;
117     float frame_prob;
118     float alpha, alphaE;
119     float frame_loudness;
120     celt_encoder_ctl(celt_enc, CELT_GET_MODE(&mode));
121
122     tonal->last_transition++;
123     alpha = 1.f/IMIN(20, 1+tonal->count);
124     alphaE = 1.f/IMIN(50, 1+tonal->count);
125
126     if (tonal->count<4)
127        tonal->music_prob = .5;
128     kfft = mode->mdct.kfft[0];
129     if (C==1)
130     {
131        for (i=0;i<N2;i++)
132        {
133           float w = .5-.5*cos(M_PI*(i+1)/N2);
134           in[i].r = MULT16_16(w, x[i]);
135           in[i].i = MULT16_16(w, x[N-N2+i]);
136           in[N-i-1].r = MULT16_16(w, x[N-i-1]);
137           in[N-i-1].i = MULT16_16(w, x[2*N-N2-i-1]);
138        }
139     } else {
140        for (i=0;i<N2;i++)
141        {
142           float w = .5-.5*cos(M_PI*(i+1)/N2);
143           in[i].r = MULT16_16(w, x[2*i]+x[2*i+1]);
144           in[i].i = MULT16_16(w, x[2*(N-N2+i)]+x[2*(N-N2+i)+1]);
145           in[N-i-1].r = MULT16_16(w, x[2*(N-i-1)]+x[2*(N-i-1)+1]);
146           in[N-i-1].i = MULT16_16(w, x[2*(2*N-N2-i-1)]+x[2*(2*N-N2-i-1)+1]);
147        }
148     }
149     opus_fft(kfft, in, out);
150
151     for (i=1;i<N2;i++)
152     {
153        float X1r, X2r, X1i, X2i;
154        float angle, d_angle, d2_angle;
155        float angle2, d_angle2, d2_angle2;
156        float mod1, mod2, avg_mod;
157        X1r = out[i].r+out[N-i].r;
158        X1i = out[i].i-out[N-i].i;
159        X2r = out[i].i+out[N-i].i;
160        X2i = out[N-i].r-out[i].r;
161
162        angle = (.5/M_PI)*atan2(X1i, X1r);
163        d_angle = angle - A[i];
164        d2_angle = d_angle - dA[i];
165
166        angle2 = (.5/M_PI)*atan2(X2i, X2r);
167        d_angle2 = angle2 - angle;
168        d2_angle2 = d_angle2 - d_angle;
169
170        mod1 = d2_angle - floor(.5+d2_angle);
171        noisiness[i] = fabs(mod1);
172        mod1 *= mod1;
173        mod1 *= mod1;
174
175        mod2 = d2_angle2 - floor(.5+d2_angle2);
176        noisiness[i] += fabs(mod2);
177        mod2 *= mod2;
178        mod2 *= mod2;
179
180        avg_mod = .25*(d2A[i]+2*mod1+mod2);
181        tonality[i] = 1./(1+40*16*pi4*avg_mod)-.015;
182
183        A[i] = angle2;
184        dA[i] = d_angle2;
185        d2A[i] = mod2;
186     }
187
188     frame_tonality = 0;
189     info->activity = 0;
190     frame_noisiness = 0;
191     frame_stationarity = 0;
192     if (!tonal->count)
193     {
194        for (b=0;b<NB_TBANDS;b++)
195        {
196           tonal->lowE[b] = 1e10;
197           tonal->highE[b] = -1e10;
198        }
199     }
200     relativeE = 0;
201     info->boost_amount[0]=info->boost_amount[1]=0;
202     info->boost_band[0]=info->boost_band[1]=0;
203     frame_loudness = 0;
204     for (b=0;b<NB_TBANDS;b++)
205     {
206        float E=0, tE=0, nE=0;
207        float L1, L2;
208        float stationarity;
209        for (i=tbands[b];i<tbands[b+1];i++)
210        {
211           float binE = out[i].r*out[i].r + out[N-i].r*out[N-i].r
212                      + out[i].i*out[i].i + out[N-i].i*out[N-i].i;
213           E += binE;
214           tE += binE*tonality[i];
215           nE += binE*2*(.5-noisiness[i]);
216        }
217        tonal->E[tonal->E_count][b] = E;
218        frame_noisiness += nE/(1e-15+E);
219
220        frame_loudness += sqrt(E+1e-10);
221        /* Add a reasonable noise floor */
222        logE[b] = log(E+1e-10);
223        tonal->lowE[b] = MIN32(logE[b], tonal->lowE[b]+.01);
224        tonal->highE[b] = MAX32(logE[b], tonal->highE[b]-.1);
225        if (tonal->highE[b] < tonal->lowE[b]+1)
226        {
227           tonal->highE[b]+=.5;
228           tonal->lowE[b]-=.5;
229        }
230        relativeE += (logE[b]-tonal->lowE[b])/(EPSILON+tonal->highE[b]-tonal->lowE[b]);
231
232        L1=L2=0;
233        for (i=0;i<NB_FRAMES;i++)
234        {
235           L1 += sqrt(tonal->E[i][b]);
236           L2 += tonal->E[i][b];
237        }
238
239        stationarity = MIN16(0.99,L1/sqrt(EPSILON+NB_FRAMES*L2));
240        stationarity *= stationarity;
241        stationarity *= stationarity;
242        frame_stationarity += stationarity;
243        /*band_tonality[b] = tE/(1e-15+E)*/;
244        band_tonality[b] = MAX16(tE/(EPSILON+E), stationarity*tonal->prev_band_tonality[b]);
245        if (b>=NB_TONAL_SKIP_BANDS)
246           frame_tonality += band_tonality[b];
247        slope += band_tonality[b]*(b-8);
248        if (band_tonality[b] > info->boost_amount[1] && b>=7 && b < NB_TBANDS-1)
249        {
250           if (band_tonality[b] > info->boost_amount[0])
251           {
252              info->boost_amount[1] = info->boost_amount[0];
253              info->boost_band[1] = info->boost_band[0];
254              info->boost_amount[0] = band_tonality[b];
255              info->boost_band[0] = b;
256           } else {
257              info->boost_amount[1] = band_tonality[b];
258              info->boost_band[1] = b;
259           }
260        }
261        tonal->prev_band_tonality[b] = band_tonality[b];
262     }
263     frame_loudness = 20*log10(frame_loudness);
264     tonal->Etracker = MAX32(tonal->Etracker-.03, frame_loudness);
265     tonal->lowECount *= (1-alphaE);
266     if (frame_loudness < tonal->Etracker-30)
267        tonal->lowECount += alphaE;
268
269     for (i=0;i<8;i++)
270     {
271        float sum=0;
272        for (b=0;b<16;b++)
273           sum += dct_table[i*16+b]*logE[b];
274        BFCC[i] = sum;
275     }
276
277     frame_stationarity /= NB_TBANDS;
278     relativeE /= NB_TBANDS;
279     if (tonal->count<10)
280        relativeE = .5;
281     frame_noisiness /= NB_TBANDS;
282 #if 1
283     info->activity = frame_noisiness + (1-frame_noisiness)*relativeE;
284 #else
285     info->activity = .5*(1+frame_noisiness-frame_stationarity);
286 #endif
287     frame_tonality /= NB_TBANDS-NB_TONAL_SKIP_BANDS;
288     frame_tonality = MAX16(frame_tonality, tonal->prev_tonality*.8);
289     tonal->prev_tonality = frame_tonality;
290     info->boost_amount[0] -= frame_tonality+.2;
291     info->boost_amount[1] -= frame_tonality+.2;
292     if (band_tonality[info->boost_band[0]] < band_tonality[info->boost_band[0]+1]+.15
293         || band_tonality[info->boost_band[0]] < band_tonality[info->boost_band[0]-1]+.15)
294        info->boost_amount[0]=0;
295     if (band_tonality[info->boost_band[1]] < band_tonality[info->boost_band[1]+1]+.15
296         || band_tonality[info->boost_band[1]] < band_tonality[info->boost_band[1]-1]+.15)
297        info->boost_amount[1]=0;
298
299     slope /= 8*8;
300     info->tonality_slope = slope;
301
302     tonal->E_count = (tonal->E_count+1)%NB_FRAMES;
303     tonal->count++;
304     info->tonality = frame_tonality;
305
306     for (i=0;i<4;i++)
307        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];
308
309     for (i=0;i<4;i++)
310        tonal->cmean[i] = (1-alpha)*tonal->cmean[i] + alpha*BFCC[i];
311
312     for (i=0;i<4;i++)
313         features[4+i] = 0.63246*(BFCC[i]-tonal->mem[i+24]) + 0.31623*(tonal->mem[i]-tonal->mem[i+16]);
314     for (i=0;i<3;i++)
315         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];
316
317     if (tonal->count > 5)
318     {
319        for (i=0;i<9;i++)
320           tonal->std[i] = (1-alpha)*tonal->std[i] + alpha*features[i]*features[i];
321     }
322
323     for (i=0;i<8;i++)
324     {
325        tonal->mem[i+24] = tonal->mem[i+16];
326        tonal->mem[i+16] = tonal->mem[i+8];
327        tonal->mem[i+8] = tonal->mem[i];
328        tonal->mem[i] = BFCC[i];
329     }
330     for (i=0;i<9;i++)
331        features[11+i] = sqrt(tonal->std[i]);
332     features[20] = info->tonality;
333     features[21] = info->activity;
334     features[22] = frame_stationarity;
335     features[23] = info->tonality_slope;
336     features[24] = tonal->lowECount;
337
338 #ifndef FIXED_POINT
339     mlp_process(&net, features, &frame_prob);
340     /* Adds a "probability dead zone", with a cap on certainty */
341     frame_prob = .90*frame_prob*frame_prob*frame_prob;
342
343     frame_prob = .5*(frame_prob+1);
344
345     /*printf("%f\n", frame_prob);*/
346     {
347        float tau, beta;
348        float p0, p1;
349        float max_certainty;
350        /* One transition every 3 minutes */
351        tau = .00005;
352        beta = .1;
353        max_certainty = 1.f/(10+1*tonal->last_transition);
354        p0 = (1-tonal->music_prob)*(1-tau) +    tonal->music_prob *tau;
355        p1 =    tonal->music_prob *(1-tau) + (1-tonal->music_prob)*tau;
356        p0 *= pow(1-frame_prob, beta);
357        p1 *= pow(frame_prob, beta);
358        tonal->music_prob = MAX16(max_certainty, MIN16(1-max_certainty, p1/(p0+p1)));
359        info->music_prob = tonal->music_prob;
360        /*printf("%f %f\n", frame_prob, info->music_prob);*/
361     }
362     if (tonal->last_music != (tonal->music_prob>.5))
363        tonal->last_transition=0;
364     tonal->last_music = tonal->music_prob>.5;
365 #else
366     info->music_prob = 0;
367 #endif
368     /*for (i=0;i<25;i++)
369        printf("%f ", features[i]);
370     printf("\n");*/
371
372     info->valid = 1;
373 }