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