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