1 /* Copyright (c) 2017-2018 Mozilla */
2 /*
3 Redistribution and use in source and binary forms, with or without
4 modification, are permitted provided that the following conditions
5 are met:
6
7 - Redistributions of source code must retain the above copyright
8 notice, this list of conditions and the following disclaimer.
9
10 - Redistributions in binary form must reproduce the above copyright
11 notice, this list of conditions and the following disclaimer in the
12 documentation and/or other materials provided with the distribution.
13
14 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
15 ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
16 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
17 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR
18 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
19 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
20 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
21 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
22 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
23 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 */
26
27 #ifdef HAVE_CONFIG_H
28 #include "config.h"
29 #endif
30
31 #include <stdlib.h>
32 #include <string.h>
33 #include <stdio.h>
34 #include "kiss_fft.h"
35 #include <math.h>
36 #include "freq.h"
37 #include "pitch.h"
38 #include "arch.h"
39 #include "burg.h"
40 #include <assert.h>
41 #include "os_support.h"
42
43 #define SQUARE(x) ((x)*(x))
44
45 static const opus_int16 eband5ms[] = {
46 /*0 200 400 600 800 1k 1.2 1.4 1.6 2k 2.4 2.8 3.2 4k 4.8 5.6 6.8 8k*/
47 0, 1, 2, 3, 4, 5, 6, 7, 8, 10, 12, 14, 16, 20, 24, 28, 34, 40
48 };
49
50 static const float compensation[] = {
51 0.8f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 0.666667f, 0.5f, 0.5f, 0.5f, 0.333333f, 0.25f, 0.25f, 0.2f, 0.166667f, 0.173913f
52 };
53
54
55 extern const kiss_fft_state kfft;
56 extern const float half_window[OVERLAP_SIZE];
57 extern const float dct_table[NB_BANDS*NB_BANDS];
58
59
compute_band_energy_inverse(float * bandE,const kiss_fft_cpx * X)60 static void compute_band_energy_inverse(float *bandE, const kiss_fft_cpx *X) {
61 int i;
62 float sum[NB_BANDS] = {0};
63 for (i=0;i<NB_BANDS-1;i++)
64 {
65 int j;
66 int band_size;
67 band_size = (eband5ms[i+1]-eband5ms[i])*WINDOW_SIZE_5MS;
68 for (j=0;j<band_size;j++) {
69 float tmp;
70 float frac = (float)j/band_size;
71 tmp = SQUARE(X[(eband5ms[i]*WINDOW_SIZE_5MS) + j].r);
72 tmp += SQUARE(X[(eband5ms[i]*WINDOW_SIZE_5MS) + j].i);
73 tmp = 1.f/(tmp + 1e-9);
74 sum[i] += (1-frac)*tmp;
75 sum[i+1] += frac*tmp;
76 }
77 }
78 sum[0] *= 2;
79 sum[NB_BANDS-1] *= 2;
80 for (i=0;i<NB_BANDS;i++)
81 {
82 bandE[i] = sum[i];
83 }
84 }
85
lpcn_lpc(opus_val16 * lpc,opus_val16 * rc,const opus_val32 * ac,int p)86 static float lpcn_lpc(
87 opus_val16 *lpc, /* out: [0...p-1] LPC coefficients */
88 opus_val16 *rc,
89 const opus_val32 *ac, /* in: [0...p] autocorrelation values */
90 int p
91 )
92 {
93 int i, j;
94 opus_val32 r;
95 opus_val32 error = ac[0];
96
97 OPUS_CLEAR(lpc, p);
98 OPUS_CLEAR(rc, p);
99 if (ac[0] != 0)
100 {
101 for (i = 0; i < p; i++) {
102 /* Sum up this iteration's reflection coefficient */
103 opus_val32 rr = 0;
104 for (j = 0; j < i; j++)
105 rr += MULT32_32_Q31(lpc[j],ac[i - j]);
106 rr += SHR32(ac[i + 1],3);
107 r = -SHL32(rr,3)/error;
108 rc[i] = r;
109 /* Update LPC coefficients and total error */
110 lpc[i] = SHR32(r,3);
111 for (j = 0; j < (i+1)>>1; j++)
112 {
113 opus_val32 tmp1, tmp2;
114 tmp1 = lpc[j];
115 tmp2 = lpc[i-1-j];
116 lpc[j] = tmp1 + MULT32_32_Q31(r,tmp2);
117 lpc[i-1-j] = tmp2 + MULT32_32_Q31(r,tmp1);
118 }
119
120 error = error - MULT32_32_Q31(MULT32_32_Q31(r,r),error);
121 /* Bail out once we get 30 dB gain */
122 if (error<.001f*ac[0])
123 break;
124 }
125 }
126 return error;
127 }
128
129
130
lpcn_compute_band_energy(float * bandE,const kiss_fft_cpx * X)131 void lpcn_compute_band_energy(float *bandE, const kiss_fft_cpx *X) {
132 int i;
133 float sum[NB_BANDS] = {0};
134 for (i=0;i<NB_BANDS-1;i++)
135 {
136 int j;
137 int band_size;
138 band_size = (eband5ms[i+1]-eband5ms[i])*WINDOW_SIZE_5MS;
139 for (j=0;j<band_size;j++) {
140 float tmp;
141 float frac = (float)j/band_size;
142 tmp = SQUARE(X[(eband5ms[i]*WINDOW_SIZE_5MS) + j].r);
143 tmp += SQUARE(X[(eband5ms[i]*WINDOW_SIZE_5MS) + j].i);
144 sum[i] += (1-frac)*tmp;
145 sum[i+1] += frac*tmp;
146 }
147 }
148 sum[0] *= 2;
149 sum[NB_BANDS-1] *= 2;
150 for (i=0;i<NB_BANDS;i++)
151 {
152 bandE[i] = sum[i];
153 }
154 }
155
compute_burg_cepstrum(const float * pcm,float * burg_cepstrum,int len,int order)156 static void compute_burg_cepstrum(const float *pcm, float *burg_cepstrum, int len, int order) {
157 int i;
158 float burg_in[FRAME_SIZE];
159 float burg_lpc[LPC_ORDER];
160 float x[WINDOW_SIZE];
161 float Eburg[NB_BANDS];
162 float g;
163 kiss_fft_cpx LPC[FREQ_SIZE];
164 float Ly[NB_BANDS];
165 float logMax = -2;
166 float follow = -2;
167 assert(order <= LPC_ORDER);
168 assert(len <= FRAME_SIZE);
169 for (i=0;i<len-1;i++) burg_in[i] = pcm[i+1] - PREEMPHASIS*pcm[i];
170 g = silk_burg_analysis(burg_lpc, burg_in, 1e-3, len-1, 1, order);
171 g /= len - 2*(order-1);
172 OPUS_CLEAR(x, WINDOW_SIZE);
173 x[0] = 1;
174 for (i=0;i<order;i++) x[i+1] = -burg_lpc[i]*pow(.995, i+1);
175 forward_transform(LPC, x);
176 compute_band_energy_inverse(Eburg, LPC);
177 for (i=0;i<NB_BANDS;i++) Eburg[i] *= .45*g*(1.f/((float)WINDOW_SIZE*WINDOW_SIZE*WINDOW_SIZE));
178 for (i=0;i<NB_BANDS;i++) {
179 Ly[i] = log10(1e-2+Eburg[i]);
180 Ly[i] = MAX16(logMax-8, MAX16(follow-2.5, Ly[i]));
181 logMax = MAX16(logMax, Ly[i]);
182 follow = MAX16(follow-2.5, Ly[i]);
183 }
184 dct(burg_cepstrum, Ly);
185 burg_cepstrum[0] += - 4;
186 }
187
burg_cepstral_analysis(float * ceps,const float * x)188 void burg_cepstral_analysis(float *ceps, const float *x) {
189 int i;
190 compute_burg_cepstrum(x, &ceps[0 ], FRAME_SIZE/2, LPC_ORDER);
191 compute_burg_cepstrum(&x[FRAME_SIZE/2], &ceps[NB_BANDS], FRAME_SIZE/2, LPC_ORDER);
192 for (i=0;i<NB_BANDS;i++) {
193 float c0, c1;
194 c0 = ceps[i];
195 c1 = ceps[NB_BANDS+i];
196 ceps[i ] = .5*(c0+c1);
197 ceps[NB_BANDS+i] = (c0-c1);
198 }
199 }
200
201
interp_band_gain(float * g,const float * bandE)202 static void interp_band_gain(float *g, const float *bandE) {
203 int i;
204 memset(g, 0, FREQ_SIZE);
205 for (i=0;i<NB_BANDS-1;i++)
206 {
207 int j;
208 int band_size;
209 band_size = (eband5ms[i+1]-eband5ms[i])*WINDOW_SIZE_5MS;
210 for (j=0;j<band_size;j++) {
211 float frac = (float)j/band_size;
212 g[(eband5ms[i]*WINDOW_SIZE_5MS) + j] = (1-frac)*bandE[i] + frac*bandE[i+1];
213 }
214 }
215 }
216
217
dct(float * out,const float * in)218 void dct(float *out, const float *in) {
219 int i;
220 for (i=0;i<NB_BANDS;i++) {
221 int j;
222 float sum = 0;
223 for (j=0;j<NB_BANDS;j++) {
224 sum += in[j] * dct_table[j*NB_BANDS + i];
225 }
226 out[i] = sum*sqrt(2./NB_BANDS);
227 }
228 }
229
idct(float * out,const float * in)230 static void idct(float *out, const float *in) {
231 int i;
232 for (i=0;i<NB_BANDS;i++) {
233 int j;
234 float sum = 0;
235 for (j=0;j<NB_BANDS;j++) {
236 sum += in[j] * dct_table[i*NB_BANDS + j];
237 }
238 out[i] = sum*sqrt(2./NB_BANDS);
239 }
240 }
241
forward_transform(kiss_fft_cpx * out,const float * in)242 void forward_transform(kiss_fft_cpx *out, const float *in) {
243 int i;
244 kiss_fft_cpx x[WINDOW_SIZE];
245 kiss_fft_cpx y[WINDOW_SIZE];
246 for (i=0;i<WINDOW_SIZE;i++) {
247 x[i].r = in[i];
248 x[i].i = 0;
249 }
250 opus_fft(&kfft, x, y, 0);
251 for (i=0;i<FREQ_SIZE;i++) {
252 out[i] = y[i];
253 }
254 }
255
inverse_transform(float * out,const kiss_fft_cpx * in)256 static void inverse_transform(float *out, const kiss_fft_cpx *in) {
257 int i;
258 kiss_fft_cpx x[WINDOW_SIZE];
259 kiss_fft_cpx y[WINDOW_SIZE];
260 for (i=0;i<FREQ_SIZE;i++) {
261 x[i] = in[i];
262 }
263 for (;i<WINDOW_SIZE;i++) {
264 x[i].r = x[WINDOW_SIZE - i].r;
265 x[i].i = -x[WINDOW_SIZE - i].i;
266 }
267 opus_fft(&kfft, x, y, 0);
268 /* output in reverse order for IFFT. */
269 out[0] = WINDOW_SIZE*y[0].r;
270 for (i=1;i<WINDOW_SIZE;i++) {
271 out[i] = WINDOW_SIZE*y[WINDOW_SIZE - i].r;
272 }
273 }
274
lpc_from_bands(float * lpc,const float * Ex)275 static float lpc_from_bands(float *lpc, const float *Ex)
276 {
277 int i;
278 float e;
279 float ac[LPC_ORDER+1];
280 float rc[LPC_ORDER];
281 float Xr[FREQ_SIZE];
282 kiss_fft_cpx X_auto[FREQ_SIZE];
283 float x_auto[WINDOW_SIZE];
284 interp_band_gain(Xr, Ex);
285 Xr[FREQ_SIZE-1] = 0;
286 OPUS_CLEAR(X_auto, FREQ_SIZE);
287 for (i=0;i<FREQ_SIZE;i++) X_auto[i].r = Xr[i];
288 inverse_transform(x_auto, X_auto);
289 for (i=0;i<LPC_ORDER+1;i++) ac[i] = x_auto[i];
290
291 /* -40 dB noise floor. */
292 ac[0] += ac[0]*1e-4 + 320/12/38.;
293 /* Lag windowing. */
294 for (i=1;i<LPC_ORDER+1;i++) ac[i] *= (1 - 6e-5*i*i);
295 e = lpcn_lpc(lpc, rc, ac, LPC_ORDER);
296 return e;
297 }
298
lpc_weighting(float * lpc,float gamma)299 void lpc_weighting(float *lpc, float gamma)
300 {
301 int i;
302 float gamma_i = gamma;
303 for (i = 0; i < LPC_ORDER; i++)
304 {
305 lpc[i] *= gamma_i;
306 gamma_i *= gamma;
307 }
308 }
309
lpc_from_cepstrum(float * lpc,const float * cepstrum)310 float lpc_from_cepstrum(float *lpc, const float *cepstrum)
311 {
312 int i;
313 float Ex[NB_BANDS];
314 float tmp[NB_BANDS];
315 OPUS_COPY(tmp, cepstrum, NB_BANDS);
316 tmp[0] += 4;
317 idct(Ex, tmp);
318 for (i=0;i<NB_BANDS;i++) Ex[i] = pow(10.f, Ex[i])*compensation[i];
319 return lpc_from_bands(lpc, Ex);
320 }
321
apply_window(float * x)322 void apply_window(float *x) {
323 int i;
324 for (i=0;i<OVERLAP_SIZE;i++) {
325 x[i] *= half_window[i];
326 x[WINDOW_SIZE - 1 - i] *= half_window[i];
327 }
328 }
329