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