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