Adds an analysis function to control VBR
[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    int E_count;
54 } TonalityAnalysisState;
55
56 int boost_band[2];
57 float boost_amount[2];
58
59 float tonality_analysis(TonalityAnalysisState *tonal, CELTEncoder *celt_enc, const opus_val16 *x, int C, float *tslope)
60 {
61     int i, b;
62     const CELTMode *mode;
63     const kiss_fft_state *kfft;
64     kiss_fft_cpx in[480], out[480];
65     const opus_val16 *window;
66     int overlap = 240;
67     int N = 480, N2=240;
68     float * restrict A = tonal->angle;
69     float * restrict dA = tonal->d_angle;
70     float * restrict d2A = tonal->d2_angle;
71     float tonality[240];
72     float band_tonality[NB_TBANDS];
73     float frame_tonality;
74     const float pi4 = M_PI*M_PI*M_PI*M_PI;
75     float slope=0;
76     float max_tonality=-1;
77     int max_band=0;
78     celt_encoder_ctl(celt_enc, CELT_GET_MODE(&mode));
79
80     kfft = mode->mdct.kfft[0];
81     window = mode->window;
82     if (C==1)
83     {
84        for (i=0;i<N2;i++)
85        {
86           float w = .5-.5*cos(M_PI*(i+1)/N2);
87           in[i].r = MULT16_16(w, x[i]);
88           in[i].i = MULT16_16(w, x[N-N2+i]);
89           in[N-i-1].r = MULT16_16(w, x[N-i-1]);
90           in[N-i-1].i = MULT16_16(w, x[2*N-N2-i-1]);
91        }
92     } else {
93        for (i=0;i<N2;i++)
94        {
95           float w = .5-.5*cos(M_PI*(i+1)/N2);
96           in[i].r = MULT16_16(w, x[2*i]+x[2*i+1]);
97           in[i].i = MULT16_16(w, x[2*(N-N2+i)]+x[2*(N-N2+i)+1]);
98           in[N-i-1].r = MULT16_16(w, x[2*(N-i-1)]+x[2*(N-i-1)+1]);
99           in[N-i-1].i = MULT16_16(w, x[2*(2*N-N2-i-1)]+x[2*(2*N-N2-i-1)+1]);
100        }
101     }
102     opus_fft(kfft, in, out);
103
104     for (i=1;i<N2;i++)
105     {
106        float X1r, X2r, X1i, X2i;
107        float angle, d_angle, d2_angle;
108        float angle2, d_angle2, d2_angle2;
109        float mod1, mod2, avg_mod;
110        X1r = out[i].r+out[N-i].r;
111        X1i = out[i].i-out[N-i].i;
112        X2r = out[i].i+out[N-i].i;
113        X2i = out[N-i].r-out[i].r;
114        //printf("%f\n", X1r);
115        angle = (.5/M_PI)*atan2(X1i, X1r);
116        d_angle = angle - A[i];
117        d2_angle = d_angle - dA[i];
118
119        angle2 = (.5/M_PI)*atan2(X2i, X2r);
120        d_angle2 = angle2 - angle;
121        d2_angle2 = d_angle2 - d_angle;
122        //printf("%f ", angle2);
123
124        //printf("%f ", d2_angle);
125        mod1 = d2_angle - floor(.5+d2_angle);
126        //printf("%f ", mod1);
127        mod1 *= mod1;
128        mod1 *= mod1;
129        mod2 = d2_angle2 - floor(.5+d2_angle2);
130        mod2 *= mod2;
131        mod2 *= mod2;
132
133        avg_mod = .25*(d2A[i]+2*mod1+mod2);
134        tonality[i] = 1./(1+40*16*pi4*avg_mod)-.015;
135
136        A[i] = angle2;
137        dA[i] = d_angle2;
138        d2A[i] = mod2;
139     }
140
141     frame_tonality = 0;
142     for (b=0;b<NB_TBANDS;b++)
143     {
144        float E=0, tE=0;
145        float L1, L2;
146        float stationarity;
147        for (i=tbands[b];i<tbands[b+1];i++)
148        {
149           float binE = out[i].r*out[i].r + out[N-i].r*out[N-i].r
150                      + out[i].i*out[i].i + out[N-i].i*out[N-i].i;
151           E += binE;
152           tE += binE*tonality[i];
153        }
154        tonal->E[tonal->E_count][b] = E;
155        L1=L2=0;
156        for (i=0;i<NB_FRAMES;i++)
157        {
158           L1 += sqrt(tonal->E[i][b]);
159           L2 += tonal->E[i][b];
160        }
161
162        stationarity = MIN16(0.99,L1/sqrt(EPSILON+NB_FRAMES*L2));
163        stationarity *= stationarity;
164        stationarity *= stationarity;
165        //fprintf(stderr, "%f %f %f\n", L1, L2, stationarity);
166        //fprintf(stderr, "%f %f\n", tE, E);
167        //fprintf(stderr, "%f %f\n", stationarity, );
168        //band_tonality[b] = tE/(1e-15+E);
169        band_tonality[b] = MAX16(tE/(1e-15+E), stationarity*tonal->prev_band_tonality[b]);
170        //if (band_tonality[b]>1)
171        //   printf("%f %f %f\n", L1, L2, stationarity);
172        //fprintf(stdout, "%f ", band_tonality[b]);
173        if (b>=7)
174           frame_tonality += band_tonality[b];
175        slope += band_tonality[b]*(b-8);
176        if (band_tonality[b] > boost_amount[1] && b>=7 && b < NB_TBANDS-1)
177        {
178           if (band_tonality[b] > boost_amount[0])
179           {
180              boost_amount[1] = boost_amount[0];
181              boost_band[1] = boost_band[0];
182              boost_amount[0] = band_tonality[b];
183              boost_band[0] = b;
184           } else {
185              boost_amount[1] = band_tonality[b];
186              boost_band[1] = b;
187           }
188        }
189        tonal->prev_band_tonality[b] = band_tonality[b];
190     }
191     frame_tonality /= NB_TBANDS-7;
192     frame_tonality = MAX16(frame_tonality, tonal->prev_tonality*.8);
193     //fprintf(stdout, "%f\n", frame_tonality);
194     tonal->prev_tonality = frame_tonality;
195     boost_amount[0] -= frame_tonality+.2;
196     boost_amount[1] -= frame_tonality+.2;
197     if (band_tonality[boost_band[0]] < band_tonality[boost_band[0]+1]+.15
198         || band_tonality[boost_band[0]] < band_tonality[boost_band[0]-1]+.15)
199        boost_amount[0]=0;
200     if (band_tonality[boost_band[1]] < band_tonality[boost_band[1]+1]+.15
201         || band_tonality[boost_band[1]] < band_tonality[boost_band[1]-1]+.15)
202        boost_amount[1]=0;
203
204     //boost_band = 16;
205     //boost_amount = .6;
206     //printf("%d %f %f\n", max_band, max_tonality, frame_tonality);
207     slope /= 8*8;
208     *tslope = slope;
209     //fprintf(stdout, "%f %f\n", frame_tonality, slope);
210
211     tonal->E_count = (tonal->E_count+1)%NB_FRAMES;
212     return frame_tonality;
213 }