6c1db2f329652777262dab1be3d5f04dcae7cb28
[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
39 #ifndef M_PI
40 #define M_PI 3.141592653
41 #endif
42
43 float dct_table[128] = {
44         0.250000, 0.250000, 0.250000, 0.250000, 0.250000, 0.250000, 0.250000, 0.250000,
45         0.250000, 0.250000, 0.250000, 0.250000, 0.250000, 0.250000, 0.250000, 0.250000,
46         0.351851, 0.338330, 0.311806, 0.273300, 0.224292, 0.166664, 0.102631, 0.034654,
47         -0.034654, -0.102631, -0.166664, -0.224292, -0.273300, -0.311806, -0.338330, -0.351851,
48         0.346760, 0.293969, 0.196424, 0.068975, -0.068975, -0.196424, -0.293969, -0.346760,
49         -0.346760, -0.293969, -0.196424, -0.068975, 0.068975, 0.196424, 0.293969, 0.346760,
50         0.338330, 0.224292, 0.034654, -0.166664, -0.311806, -0.351851, -0.273300, -0.102631,
51         0.102631, 0.273300, 0.351851, 0.311806, 0.166664, -0.034654, -0.224292, -0.338330,
52         0.326641, 0.135299, -0.135299, -0.326641, -0.326641, -0.135299, 0.135299, 0.326641,
53         0.326641, 0.135299, -0.135299, -0.326641, -0.326641, -0.135299, 0.135299, 0.326641,
54         0.311806, 0.034654, -0.273300, -0.338330, -0.102631, 0.224292, 0.351851, 0.166664,
55         -0.166664, -0.351851, -0.224292, 0.102631, 0.338330, 0.273300, -0.034654, -0.311806,
56         0.293969, -0.068975, -0.346760, -0.196424, 0.196424, 0.346760, 0.068975, -0.293969,
57         -0.293969, 0.068975, 0.346760, 0.196424, -0.196424, -0.346760, -0.068975, 0.293969,
58         0.273300, -0.166664, -0.338330, 0.034654, 0.351851, 0.102631, -0.311806, -0.224292,
59         0.224292, 0.311806, -0.102631, -0.351851, -0.034654, 0.338330, 0.166664, -0.273300,
60 };
61
62 #define NB_FRAMES 8
63
64 #define NB_TBANDS 18
65 static const int tbands[NB_TBANDS+1] = {
66       2, 4, 6, 8, 10, 12, 14, 16, 20, 24, 28, 32, 40, 48, 56, 68, 80, 96, 120
67 };
68
69 #define NB_TONAL_SKIP_BANDS 8
70
71 typedef struct {
72    float angle[240];
73    float d_angle[240];
74    float d2_angle[240];
75    float prev_band_tonality[NB_TBANDS];
76    float prev_tonality;
77    float E[NB_FRAMES][NB_TBANDS];
78    float lowE[NB_TBANDS], highE[NB_TBANDS];
79    float mem[32];
80    int E_count;
81    int count;
82 } TonalityAnalysisState;
83
84 void tonality_analysis(TonalityAnalysisState *tonal, AnalysisInfo *info, CELTEncoder *celt_enc, const opus_val16 *x, int C)
85 {
86     int i, b;
87     const CELTMode *mode;
88     const kiss_fft_state *kfft;
89     kiss_fft_cpx in[480], out[480];
90     int N = 480, N2=240;
91     float * restrict A = tonal->angle;
92     float * restrict dA = tonal->d_angle;
93     float * restrict d2A = tonal->d2_angle;
94     float tonality[240];
95     float noisiness[240];
96     float band_tonality[NB_TBANDS];
97     float logE[NB_TBANDS];
98     float BFCC[8];
99     float features[27];
100     float frame_tonality;
101     float frame_noisiness;
102     const float pi4 = M_PI*M_PI*M_PI*M_PI;
103     float slope=0;
104     float frame_stationarity;
105     float relativeE;
106     celt_encoder_ctl(celt_enc, CELT_GET_MODE(&mode));
107
108     kfft = mode->mdct.kfft[0];
109     if (C==1)
110     {
111        for (i=0;i<N2;i++)
112        {
113           float w = .5-.5*cos(M_PI*(i+1)/N2);
114           in[i].r = MULT16_16(w, x[i]);
115           in[i].i = MULT16_16(w, x[N-N2+i]);
116           in[N-i-1].r = MULT16_16(w, x[N-i-1]);
117           in[N-i-1].i = MULT16_16(w, x[2*N-N2-i-1]);
118        }
119     } else {
120        for (i=0;i<N2;i++)
121        {
122           float w = .5-.5*cos(M_PI*(i+1)/N2);
123           in[i].r = MULT16_16(w, x[2*i]+x[2*i+1]);
124           in[i].i = MULT16_16(w, x[2*(N-N2+i)]+x[2*(N-N2+i)+1]);
125           in[N-i-1].r = MULT16_16(w, x[2*(N-i-1)]+x[2*(N-i-1)+1]);
126           in[N-i-1].i = MULT16_16(w, x[2*(2*N-N2-i-1)]+x[2*(2*N-N2-i-1)+1]);
127        }
128     }
129     opus_fft(kfft, in, out);
130
131     for (i=1;i<N2;i++)
132     {
133        float X1r, X2r, X1i, X2i;
134        float angle, d_angle, d2_angle;
135        float angle2, d_angle2, d2_angle2;
136        float mod1, mod2, avg_mod;
137        X1r = out[i].r+out[N-i].r;
138        X1i = out[i].i-out[N-i].i;
139        X2r = out[i].i+out[N-i].i;
140        X2i = out[N-i].r-out[i].r;
141
142        angle = (.5/M_PI)*atan2(X1i, X1r);
143        d_angle = angle - A[i];
144        d2_angle = d_angle - dA[i];
145
146        angle2 = (.5/M_PI)*atan2(X2i, X2r);
147        d_angle2 = angle2 - angle;
148        d2_angle2 = d_angle2 - d_angle;
149
150        mod1 = d2_angle - floor(.5+d2_angle);
151        noisiness[i] = fabs(mod1);
152        mod1 *= mod1;
153        mod1 *= mod1;
154
155        mod2 = d2_angle2 - floor(.5+d2_angle2);
156        noisiness[i] += fabs(mod2);
157        mod2 *= mod2;
158        mod2 *= mod2;
159
160        avg_mod = .25*(d2A[i]+2*mod1+mod2);
161        tonality[i] = 1./(1+40*16*pi4*avg_mod)-.015;
162
163        A[i] = angle2;
164        dA[i] = d_angle2;
165        d2A[i] = mod2;
166     }
167
168     frame_tonality = 0;
169     info->activity = 0;
170     frame_noisiness = 0;
171     frame_stationarity = 0;
172     if (!tonal->count)
173     {
174        for (b=0;b<NB_TBANDS;b++)
175        {
176           tonal->lowE[b] = 1e10;
177           tonal->highE[b] = -1e10;
178        }
179     }
180     relativeE = 0;
181     info->boost_amount[0]=info->boost_amount[1]=0;
182     info->boost_band[0]=info->boost_band[1]=0;
183     for (b=0;b<NB_TBANDS;b++)
184     {
185        float E=0, tE=0, nE=0;
186        float L1, L2;
187        float stationarity;
188        for (i=tbands[b];i<tbands[b+1];i++)
189        {
190           float binE = out[i].r*out[i].r + out[N-i].r*out[N-i].r
191                      + out[i].i*out[i].i + out[N-i].i*out[N-i].i;
192           E += binE;
193           tE += binE*tonality[i];
194           nE += binE*2*(.5-noisiness[i]);
195        }
196        tonal->E[tonal->E_count][b] = E;
197        frame_noisiness += nE/(1e-15+E);
198
199        logE[b] = log(E+EPSILON);
200        tonal->lowE[b] = MIN32(logE[b], tonal->lowE[b]+.01);
201        tonal->highE[b] = MAX32(logE[b], tonal->highE[b]-.1);
202        if (tonal->highE[b] < tonal->lowE[b]+1)
203        {
204           tonal->highE[b]+=.5;
205           tonal->lowE[b]-=.5;
206        }
207        relativeE += (logE[b]-tonal->lowE[b])/(EPSILON+tonal->highE[b]-tonal->lowE[b]);
208
209        L1=L2=0;
210        for (i=0;i<NB_FRAMES;i++)
211        {
212           L1 += sqrt(tonal->E[i][b]);
213           L2 += tonal->E[i][b];
214        }
215
216        stationarity = MIN16(0.99,L1/sqrt(EPSILON+NB_FRAMES*L2));
217        stationarity *= stationarity;
218        stationarity *= stationarity;
219        frame_stationarity += stationarity;
220        /*band_tonality[b] = tE/(1e-15+E)*/;
221        band_tonality[b] = MAX16(tE/(EPSILON+E), stationarity*tonal->prev_band_tonality[b]);
222        if (b>=NB_TONAL_SKIP_BANDS)
223           frame_tonality += band_tonality[b];
224        slope += band_tonality[b]*(b-8);
225        if (band_tonality[b] > info->boost_amount[1] && b>=7 && b < NB_TBANDS-1)
226        {
227           if (band_tonality[b] > info->boost_amount[0])
228           {
229              info->boost_amount[1] = info->boost_amount[0];
230              info->boost_band[1] = info->boost_band[0];
231              info->boost_amount[0] = band_tonality[b];
232              info->boost_band[0] = b;
233           } else {
234              info->boost_amount[1] = band_tonality[b];
235              info->boost_band[1] = b;
236           }
237        }
238        tonal->prev_band_tonality[b] = band_tonality[b];
239     }
240
241     for (i=0;i<8;i++)
242     {
243        float sum=0;
244        for (b=0;b<16;b++)
245           sum += dct_table[i*16+b]*logE[b];
246        BFCC[i] = sum;
247     }
248
249     frame_stationarity /= NB_TBANDS;
250     relativeE /= NB_TBANDS;
251     if (tonal->count<10)
252        relativeE = .5;
253     frame_noisiness /= NB_TBANDS;
254 #if 1
255     info->activity = frame_noisiness + (1-frame_noisiness)*relativeE;
256 #else
257     info->activity = .5*(1+frame_noisiness-frame_stationarity);
258 #endif
259     frame_tonality /= NB_TBANDS-NB_TONAL_SKIP_BANDS;
260     frame_tonality = MAX16(frame_tonality, tonal->prev_tonality*.8);
261     tonal->prev_tonality = frame_tonality;
262     info->boost_amount[0] -= frame_tonality+.2;
263     info->boost_amount[1] -= frame_tonality+.2;
264     if (band_tonality[info->boost_band[0]] < band_tonality[info->boost_band[0]+1]+.15
265         || band_tonality[info->boost_band[0]] < band_tonality[info->boost_band[0]-1]+.15)
266        info->boost_amount[0]=0;
267     if (band_tonality[info->boost_band[1]] < band_tonality[info->boost_band[1]+1]+.15
268         || band_tonality[info->boost_band[1]] < band_tonality[info->boost_band[1]-1]+.15)
269        info->boost_amount[1]=0;
270
271     slope /= 8*8;
272     info->tonality_slope = slope;
273
274     tonal->E_count = (tonal->E_count+1)%NB_FRAMES;
275     tonal->count++;
276     info->tonality = frame_tonality;
277
278     for (i=1;i<8;i++)
279         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];
280
281     for (i=0;i<8;i++)
282         features[7+i] = 0.63246*(BFCC[i]-tonal->mem[i+24]) + 0.31623*(tonal->mem[i]-tonal->mem[i+16]);
283     for (i=0;i<8;i++)
284         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];
285     for (i=0;i<8;i++)
286     {
287        tonal->mem[i+24] = tonal->mem[i+16];
288        tonal->mem[i+16] = tonal->mem[i+8];
289        tonal->mem[i+8] = tonal->mem[i];
290        tonal->mem[i] = BFCC[i];
291     }
292     features[23] = info->tonality;
293     features[24] = info->tonality_slope;
294     features[25] = info->activity;
295     features[26] = frame_stationarity;
296
297     /*for (i=0;i<27;i++)
298        printf("%f ", features[i]);
299     printf("\n");*/
300
301     info->valid = 1;
302 }