Cleans up the most ugly parts of the analysis code
[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 #define NB_FRAMES 8
40
41 #define NB_TBANDS 17
42 static const int tbands[NB_TBANDS+1] = {
43       4, 6, 8, 10, 12, 14, 16, 20, 24, 32, 40, 48, 56, 68, 80, 96, 120, 156
44 };
45
46 typedef struct {
47    float angle[240];
48    float d_angle[240];
49    float d2_angle[240];
50    float prev_band_tonality[NB_TBANDS];
51    float prev_tonality;
52    float E[NB_FRAMES][NB_TBANDS];
53    float lowE[NB_TBANDS], highE[NB_TBANDS];
54    int E_count;
55    int count;
56 } TonalityAnalysisState;
57
58 void tonality_analysis(TonalityAnalysisState *tonal, AnalysisInfo *info, CELTEncoder *celt_enc, const opus_val16 *x, int C)
59 {
60     int i, b;
61     const CELTMode *mode;
62     const kiss_fft_state *kfft;
63     kiss_fft_cpx in[480], out[480];
64     int N = 480, N2=240;
65     float * restrict A = tonal->angle;
66     float * restrict dA = tonal->d_angle;
67     float * restrict d2A = tonal->d2_angle;
68     float tonality[240];
69     float noisiness[240];
70     float band_tonality[NB_TBANDS];
71     float frame_tonality;
72     float frame_noisiness;
73     const float pi4 = M_PI*M_PI*M_PI*M_PI;
74     float slope=0;
75     float frame_stationarity;
76     float relativeE;
77     celt_encoder_ctl(celt_enc, CELT_GET_MODE(&mode));
78
79     kfft = mode->mdct.kfft[0];
80     if (C==1)
81     {
82        for (i=0;i<N2;i++)
83        {
84           float w = .5-.5*cos(M_PI*(i+1)/N2);
85           in[i].r = MULT16_16(w, x[i]);
86           in[i].i = MULT16_16(w, x[N-N2+i]);
87           in[N-i-1].r = MULT16_16(w, x[N-i-1]);
88           in[N-i-1].i = MULT16_16(w, x[2*N-N2-i-1]);
89        }
90     } else {
91        for (i=0;i<N2;i++)
92        {
93           float w = .5-.5*cos(M_PI*(i+1)/N2);
94           in[i].r = MULT16_16(w, x[2*i]+x[2*i+1]);
95           in[i].i = MULT16_16(w, x[2*(N-N2+i)]+x[2*(N-N2+i)+1]);
96           in[N-i-1].r = MULT16_16(w, x[2*(N-i-1)]+x[2*(N-i-1)+1]);
97           in[N-i-1].i = MULT16_16(w, x[2*(2*N-N2-i-1)]+x[2*(2*N-N2-i-1)+1]);
98        }
99     }
100     opus_fft(kfft, in, out);
101
102     for (i=1;i<N2;i++)
103     {
104        float X1r, X2r, X1i, X2i;
105        float angle, d_angle, d2_angle;
106        float angle2, d_angle2, d2_angle2;
107        float mod1, mod2, avg_mod;
108        X1r = out[i].r+out[N-i].r;
109        X1i = out[i].i-out[N-i].i;
110        X2r = out[i].i+out[N-i].i;
111        X2i = out[N-i].r-out[i].r;
112
113        angle = (.5/M_PI)*atan2(X1i, X1r);
114        d_angle = angle - A[i];
115        d2_angle = d_angle - dA[i];
116
117        angle2 = (.5/M_PI)*atan2(X2i, X2r);
118        d_angle2 = angle2 - angle;
119        d2_angle2 = d_angle2 - d_angle;
120
121        mod1 = d2_angle - floor(.5+d2_angle);
122        noisiness[i] = fabs(mod1);
123        mod1 *= mod1;
124        mod1 *= mod1;
125
126        mod2 = d2_angle2 - floor(.5+d2_angle2);
127        noisiness[i] += fabs(mod2);
128        mod2 *= mod2;
129        mod2 *= mod2;
130
131        avg_mod = .25*(d2A[i]+2*mod1+mod2);
132        tonality[i] = 1./(1+40*16*pi4*avg_mod)-.015;
133
134        A[i] = angle2;
135        dA[i] = d_angle2;
136        d2A[i] = mod2;
137     }
138
139     frame_tonality = 0;
140     info->activity = 0;
141     frame_noisiness = 0;
142     frame_stationarity = 0;
143     if (!tonal->count)
144     {
145        for (b=0;b<NB_TBANDS;b++)
146        {
147           tonal->lowE[b] = 1e10;
148           tonal->highE[b] = -1e10;
149        }
150     }
151     relativeE = 0;
152     info->boost_amount[0]=info->boost_amount[1]=0;
153     info->boost_band[0]=info->boost_band[1]=0;
154     for (b=0;b<NB_TBANDS;b++)
155     {
156        float E=0, tE=0, nE=0, logE;
157        float L1, L2;
158        float stationarity;
159        for (i=tbands[b];i<tbands[b+1];i++)
160        {
161           float binE = out[i].r*out[i].r + out[N-i].r*out[N-i].r
162                      + out[i].i*out[i].i + out[N-i].i*out[N-i].i;
163           E += binE;
164           tE += binE*tonality[i];
165           nE += binE*2*(.5-noisiness[i]);
166        }
167        tonal->E[tonal->E_count][b] = E;
168        frame_noisiness += nE/(1e-15+E);
169
170        logE = log(E+EPSILON);
171        tonal->lowE[b] = MIN32(logE, tonal->lowE[b]+.01);
172        tonal->highE[b] = MAX32(logE, tonal->highE[b]-.1);
173        if (tonal->highE[b] < tonal->lowE[b]+1)
174        {
175           tonal->highE[b]+=.5;
176           tonal->lowE[b]-=.5;
177        }
178        relativeE += (logE-tonal->lowE[b])/(EPSILON+tonal->highE[b]-tonal->lowE[b]);
179
180        L1=L2=0;
181        for (i=0;i<NB_FRAMES;i++)
182        {
183           L1 += sqrt(tonal->E[i][b]);
184           L2 += tonal->E[i][b];
185        }
186
187        stationarity = MIN16(0.99,L1/sqrt(EPSILON+NB_FRAMES*L2));
188        stationarity *= stationarity;
189        stationarity *= stationarity;
190        frame_stationarity += stationarity;
191        /*band_tonality[b] = tE/(1e-15+E)*/;
192        band_tonality[b] = MAX16(tE/(EPSILON+E), stationarity*tonal->prev_band_tonality[b]);
193        if (b>=7)
194           frame_tonality += band_tonality[b];
195        slope += band_tonality[b]*(b-8);
196        if (band_tonality[b] > info->boost_amount[1] && b>=7 && b < NB_TBANDS-1)
197        {
198           if (band_tonality[b] > info->boost_amount[0])
199           {
200              info->boost_amount[1] = info->boost_amount[0];
201              info->boost_band[1] = info->boost_band[0];
202              info->boost_amount[0] = band_tonality[b];
203              info->boost_band[0] = b;
204           } else {
205              info->boost_amount[1] = band_tonality[b];
206              info->boost_band[1] = b;
207           }
208        }
209        tonal->prev_band_tonality[b] = band_tonality[b];
210     }
211     frame_stationarity /= NB_TBANDS;
212     relativeE /= NB_TBANDS;
213     if (tonal->count<10)
214        relativeE = .5;
215     frame_noisiness /= NB_TBANDS;
216 #if 1
217     info->activity = frame_noisiness + (1-frame_noisiness)*relativeE;
218 #else
219     info->activity = .5*(1+frame_noisiness-frame_stationarity);
220 #endif
221     frame_tonality /= NB_TBANDS-7;
222     frame_tonality = MAX16(frame_tonality, tonal->prev_tonality*.8);
223     tonal->prev_tonality = frame_tonality;
224     info->boost_amount[0] -= frame_tonality+.2;
225     info->boost_amount[1] -= frame_tonality+.2;
226     if (band_tonality[info->boost_band[0]] < band_tonality[info->boost_band[0]+1]+.15
227         || band_tonality[info->boost_band[0]] < band_tonality[info->boost_band[0]-1]+.15)
228        info->boost_amount[0]=0;
229     if (band_tonality[info->boost_band[1]] < band_tonality[info->boost_band[1]+1]+.15
230         || band_tonality[info->boost_band[1]] < band_tonality[info->boost_band[1]-1]+.15)
231        info->boost_amount[1]=0;
232
233     slope /= 8*8;
234     info->tonality_slope = slope;
235
236     tonal->E_count = (tonal->E_count+1)%NB_FRAMES;
237     tonal->count++;
238     info->tonality = frame_tonality;
239     info->valid = 1;
240 }