a55e8a944c9bfbaba50e18e87481acae7894267c
[ffmpeg.git] / libavcodec / ra144enc.c
1 /*
2  * Real Audio 1.0 (14.4K) encoder
3  * Copyright (c) 2010 Francesco Lavra <francescolavra@interfree.it>
4  *
5  * This file is part of FFmpeg.
6  *
7  * FFmpeg is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public
9  * License as published by the Free Software Foundation; either
10  * version 2.1 of the License, or (at your option) any later version.
11  *
12  * FFmpeg is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
15  * Lesser General Public License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with FFmpeg; if not, write to the Free Software
19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
20  */
21
22 /**
23  * @file libavcodec/ra144enc.c
24  * Real Audio 1.0 (14.4K) encoder
25  * @author Francesco Lavra <francescolavra@interfree.it>
26  */
27
28 #include <values.h>
29 #include <float.h>
30
31 #include "avcodec.h"
32 #include "put_bits.h"
33 #include "lpc.h"
34 #include "celp_filters.h"
35 #include "ra144.h"
36
37
38 static av_cold int ra144_encode_init(AVCodecContext * avctx)
39 {
40     RA144Context *ractx;
41
42     if (avctx->sample_fmt != SAMPLE_FMT_S16) {
43         av_log(avctx, AV_LOG_ERROR, "invalid sample format\n");
44         return -1;
45     }
46     if (avctx->channels != 1) {
47         av_log(avctx, AV_LOG_ERROR, "invalid number of channels: %d\n",
48                avctx->channels);
49         return -1;
50     }
51     avctx->frame_size = NBLOCKS * BLOCKSIZE;
52     avctx->bit_rate = 8000;
53     ractx = avctx->priv_data;
54     ractx->lpc_coef[0] = ractx->lpc_tables[0];
55     ractx->lpc_coef[1] = ractx->lpc_tables[1];
56     ractx->avctx = avctx;
57     dsputil_init(&ractx->dsp, avctx);
58     return 0;
59 }
60
61
62 /**
63  * Quantizes a value by searching a sorted table for the element with the
64  * nearest value
65  *
66  * @param value value to quantize
67  * @param table array containing the quantization table
68  * @param size size of the quantization table
69  * @return index of the quantization table corresponding to the element with the
70  *         nearest value
71  */
72 static int quantize(int value, const int16_t *table, unsigned int size)
73 {
74     unsigned int low = 0, high = size - 1;
75
76     while (1) {
77         int index = (low + high) >> 1;
78         int error = table[index] - value;
79
80         if (index == low)
81             return table[high] + error > value ? low : high;
82         if (error > 0) {
83             high = index;
84         } else {
85             low = index;
86         }
87     }
88 }
89
90
91 /**
92  * Orthogonalizes a vector to another vector
93  *
94  * @param v vector to orthogonalize
95  * @param u vector against which orthogonalization is performed
96  */
97 static void orthogonalize(float *v, const float *u)
98 {
99     int i;
100     float num = 0, den = 0;
101
102     for (i = 0; i < BLOCKSIZE; i++) {
103         num += v[i] * u[i];
104         den += u[i] * u[i];
105     }
106     num /= den;
107     for (i = 0; i < BLOCKSIZE; i++)
108         v[i] -= num * u[i];
109 }
110
111
112 /**
113  * Calculates match score and gain of an LPC-filtered vector with respect to
114  * input data, possibly othogonalizing it to up to 2 other vectors
115  *
116  * @param work array used to calculate the filtered vector
117  * @param coefs coefficients of the LPC filter
118  * @param vect original vector
119  * @param ortho1 first vector against which orthogonalization is performed
120  * @param ortho2 second vector against which orthogonalization is performed
121  * @param data input data
122  * @param score pointer to variable where match score is returned
123  * @param gain pointer to variable where gain is returned
124  */
125 static void get_match_score(float *work, const float *coefs, float *vect,
126                             const float *ortho1, const float *ortho2,
127                             const float *data, float *score, float *gain)
128 {
129     float c, g;
130     int i;
131
132     ff_celp_lp_synthesis_filterf(work, coefs, vect, BLOCKSIZE, LPC_ORDER);
133     if (ortho1)
134         orthogonalize(work, ortho1);
135     if (ortho2)
136         orthogonalize(work, ortho2);
137     c = g = 0;
138     for (i = 0; i < BLOCKSIZE; i++) {
139         g += work[i] * work[i];
140         c += data[i] * work[i];
141     }
142     if (c <= 0) {
143         *score = 0;
144         return;
145     }
146     *gain = c / g;
147     *score = *gain * c;
148 }
149
150
151 /**
152  * Creates a vector from the adaptive codebook at a given lag value
153  *
154  * @param vect array where vector is stored
155  * @param cb adaptive codebook
156  * @param lag lag value
157  */
158 static void create_adapt_vect(float *vect, const int16_t *cb, int lag)
159 {
160     int i;
161
162     cb += BUFFERSIZE - lag;
163     for (i = 0; i < FFMIN(BLOCKSIZE, lag); i++)
164         vect[i] = cb[i];
165     if (lag < BLOCKSIZE)
166         for (i = 0; i < BLOCKSIZE - lag; i++)
167             vect[lag + i] = cb[i];
168 }
169
170
171 /**
172  * Searches the adaptive codebook for the best entry and gain and removes its
173  * contribution from input data
174  *
175  * @param adapt_cb array from which the adaptive codebook is extracted
176  * @param work array used to calculate LPC-filtered vectors
177  * @param coefs coefficients of the LPC filter
178  * @param data input data
179  * @return index of the best entry of the adaptive codebook
180  */
181 static int adaptive_cb_search(const int16_t *adapt_cb, float *work,
182                               const float *coefs, float *data)
183 {
184     int i, best_vect;
185     float score, gain, best_score, best_gain;
186     float exc[BLOCKSIZE];
187
188     gain = best_score = 0;
189     for (i = BLOCKSIZE / 2; i <= BUFFERSIZE; i++) {
190         create_adapt_vect(exc, adapt_cb, i);
191         get_match_score(work, coefs, exc, NULL, NULL, data, &score, &gain);
192         if (score > best_score) {
193             best_score = score;
194             best_vect = i;
195             best_gain = gain;
196         }
197     }
198     if (!best_score)
199         return 0;
200
201     /**
202      * Re-calculate the filtered vector from the vector with maximum match score
203      * and remove its contribution from input data.
204      */
205     create_adapt_vect(exc, adapt_cb, best_vect);
206     ff_celp_lp_synthesis_filterf(work, coefs, exc, BLOCKSIZE, LPC_ORDER);
207     for (i = 0; i < BLOCKSIZE; i++)
208         data[i] -= best_gain * work[i];
209     return (best_vect - BLOCKSIZE / 2 + 1);
210 }
211
212
213 /**
214  * Finds the best vector of a fixed codebook by applying an LPC filter to
215  * codebook entries, possibly othogonalizing them to up to 2 other vectors and
216  * matching the results with input data
217  *
218  * @param work array used to calculate the filtered vectors
219  * @param coefs coefficients of the LPC filter
220  * @param cb fixed codebook
221  * @param ortho1 first vector against which orthogonalization is performed
222  * @param ortho2 second vector against which orthogonalization is performed
223  * @param data input data
224  * @param idx pointer to variable where the index of the best codebook entry is
225  *        returned
226  * @param gain pointer to variable where the gain of the best codebook entry is
227  *        returned
228  */
229 static void find_best_vect(float *work, const float *coefs,
230                            const int8_t cb[][BLOCKSIZE], const float *ortho1,
231                            const float *ortho2, float *data, int *idx,
232                            float *gain)
233 {
234     int i, j;
235     float g, score, best_score;
236     float vect[BLOCKSIZE];
237
238     *idx = *gain = best_score = 0;
239     for (i = 0; i < FIXED_CB_SIZE; i++) {
240         for (j = 0; j < BLOCKSIZE; j++)
241             vect[j] = cb[i][j];
242         get_match_score(work, coefs, vect, ortho1, ortho2, data, &score, &g);
243         if (score > best_score) {
244             best_score = score;
245             *idx = i;
246             *gain = g;
247         }
248     }
249 }
250
251
252 /**
253  * Searches the two fixed codebooks for the best entry and gain
254  *
255  * @param work array used to calculate LPC-filtered vectors
256  * @param coefs coefficients of the LPC filter
257  * @param data input data
258  * @param cba_idx index of the best entry of the adaptive codebook
259  * @param cb1_idx pointer to variable where the index of the best entry of the
260  *        first fixed codebook is returned
261  * @param cb2_idx pointer to variable where the index of the best entry of the
262  *        second fixed codebook is returned
263  */
264 static void fixed_cb_search(float *work, const float *coefs, float *data,
265                             int cba_idx, int *cb1_idx, int *cb2_idx)
266 {
267     int i, ortho_cb1;
268     float gain;
269     float cba_vect[BLOCKSIZE], cb1_vect[BLOCKSIZE];
270     float vect[BLOCKSIZE];
271
272     /**
273      * The filtered vector from the adaptive codebook can be retrieved from
274      * work, because this function is called just after adaptive_cb_search().
275      */
276     if (cba_idx)
277         memcpy(cba_vect, work, sizeof(cba_vect));
278
279     find_best_vect(work, coefs, ff_cb1_vects, cba_idx ? cba_vect : NULL, NULL,
280                    data, cb1_idx, &gain);
281
282     /**
283      * Re-calculate the filtered vector from the vector with maximum match score
284      * and remove its contribution from input data.
285      */
286     if (gain) {
287         for (i = 0; i < BLOCKSIZE; i++)
288             vect[i] = ff_cb1_vects[*cb1_idx][i];
289         ff_celp_lp_synthesis_filterf(work, coefs, vect, BLOCKSIZE, LPC_ORDER);
290         if (cba_idx)
291             orthogonalize(work, cba_vect);
292         for (i = 0; i < BLOCKSIZE; i++)
293             data[i] -= gain * work[i];
294         memcpy(cb1_vect, work, sizeof(cb1_vect));
295         ortho_cb1 = 1;
296     } else
297         ortho_cb1 = 0;
298
299     find_best_vect(work, coefs, ff_cb2_vects, cba_idx ? cba_vect : NULL,
300                    ortho_cb1 ? cb1_vect : NULL, data, cb2_idx, &gain);
301 }
302
303
304 /**
305  * Encodes a subblock of the current frame
306  *
307  * @param ractx encoder context
308  * @param sblock_data input data of the subblock
309  * @param lpc_coefs coefficients of the LPC filter
310  * @param rms RMS of the reflection coefficients
311  * @param pb pointer to PutBitContext of the current frame
312  */
313 static void ra144_encode_subblock(RA144Context *ractx,
314                                   const int16_t *sblock_data,
315                                   const int16_t *lpc_coefs, unsigned int rms,
316                                   PutBitContext *pb)
317 {
318     float data[BLOCKSIZE], work[LPC_ORDER + BLOCKSIZE];
319     float coefs[LPC_ORDER];
320     float zero[BLOCKSIZE], cba[BLOCKSIZE], cb1[BLOCKSIZE], cb2[BLOCKSIZE];
321     int16_t cba_vect[BLOCKSIZE];
322     int cba_idx, cb1_idx, cb2_idx, gain;
323     int i, n, m[3];
324     float g[3];
325     float error, best_error;
326
327     for (i = 0; i < LPC_ORDER; i++) {
328         work[i] = ractx->curr_sblock[BLOCKSIZE + i];
329         coefs[i] = lpc_coefs[i] * (1/4096.0);
330     }
331
332     /**
333      * Calculate the zero-input response of the LPC filter and subtract it from
334      * input data.
335      */
336     memset(data, 0, sizeof(data));
337     ff_celp_lp_synthesis_filterf(work + LPC_ORDER, coefs, data, BLOCKSIZE,
338                                  LPC_ORDER);
339     for (i = 0; i < BLOCKSIZE; i++) {
340         zero[i] = work[LPC_ORDER + i];
341         data[i] = sblock_data[i] - zero[i];
342     }
343
344     /**
345      * Codebook search is performed without taking into account the contribution
346      * of the previous subblock, since it has been just subtracted from input
347      * data.
348      */
349     memset(work, 0, LPC_ORDER * sizeof(*work));
350
351     cba_idx = adaptive_cb_search(ractx->adapt_cb, work + LPC_ORDER, coefs,
352                                  data);
353     if (cba_idx) {
354         /**
355          * The filtered vector from the adaptive codebook can be retrieved from
356          * work, see implementation of adaptive_cb_search().
357          */
358         memcpy(cba, work + LPC_ORDER, sizeof(cba));
359
360         ff_copy_and_dup(cba_vect, ractx->adapt_cb, cba_idx + BLOCKSIZE / 2 - 1);
361         m[0] = (ff_irms(cba_vect) * rms) >> 12;
362     }
363     fixed_cb_search(work + LPC_ORDER, coefs, data, cba_idx, &cb1_idx, &cb2_idx);
364     for (i = 0; i < BLOCKSIZE; i++) {
365         cb1[i] = ff_cb1_vects[cb1_idx][i];
366         cb2[i] = ff_cb2_vects[cb2_idx][i];
367     }
368     ff_celp_lp_synthesis_filterf(work + LPC_ORDER, coefs, cb1, BLOCKSIZE,
369                                  LPC_ORDER);
370     memcpy(cb1, work + LPC_ORDER, sizeof(cb1));
371     m[1] = (ff_cb1_base[cb1_idx] * rms) >> 8;
372     ff_celp_lp_synthesis_filterf(work + LPC_ORDER, coefs, cb2, BLOCKSIZE,
373                                  LPC_ORDER);
374     memcpy(cb2, work + LPC_ORDER, sizeof(cb2));
375     m[2] = (ff_cb2_base[cb2_idx] * rms) >> 8;
376     best_error = FLT_MAX;
377     gain = 0;
378     for (n = 0; n < 256; n++) {
379         g[1] = ((ff_gain_val_tab[n][1] * m[1]) >> ff_gain_exp_tab[n]) *
380                (1/4096.0);
381         g[2] = ((ff_gain_val_tab[n][2] * m[2]) >> ff_gain_exp_tab[n]) *
382                (1/4096.0);
383         error = 0;
384         if (cba_idx) {
385             g[0] = ((ff_gain_val_tab[n][0] * m[0]) >> ff_gain_exp_tab[n]) *
386                    (1/4096.0);
387             for (i = 0; i < BLOCKSIZE; i++) {
388                 data[i] = zero[i] + g[0] * cba[i] + g[1] * cb1[i] +
389                           g[2] * cb2[i];
390                 error += (data[i] - sblock_data[i]) *
391                          (data[i] - sblock_data[i]);
392             }
393         } else {
394             for (i = 0; i < BLOCKSIZE; i++) {
395                 data[i] = zero[i] + g[1] * cb1[i] + g[2] * cb2[i];
396                 error += (data[i] - sblock_data[i]) *
397                          (data[i] - sblock_data[i]);
398             }
399         }
400         if (error < best_error) {
401             best_error = error;
402             gain = n;
403         }
404     }
405     put_bits(pb, 7, cba_idx);
406     put_bits(pb, 8, gain);
407     put_bits(pb, 7, cb1_idx);
408     put_bits(pb, 7, cb2_idx);
409     ff_subblock_synthesis(ractx, lpc_coefs, cba_idx, cb1_idx, cb2_idx, rms,
410                           gain);
411 }
412
413
414 static int ra144_encode_frame(AVCodecContext *avctx, uint8_t *frame,
415                               int buf_size, void *data)
416 {
417     static const uint8_t sizes[LPC_ORDER] = {64, 32, 32, 16, 16, 8, 8, 8, 8, 4};
418     static const uint8_t bit_sizes[LPC_ORDER] = {6, 5, 5, 4, 4, 3, 3, 3, 3, 2};
419     RA144Context *ractx;
420     PutBitContext pb;
421     int32_t lpc_data[NBLOCKS * BLOCKSIZE];
422     int32_t lpc_coefs[LPC_ORDER][MAX_LPC_ORDER];
423     int shift[LPC_ORDER];
424     int16_t block_coefs[NBLOCKS][LPC_ORDER];
425     int lpc_refl[LPC_ORDER];    /**< reflection coefficients of the frame */
426     unsigned int refl_rms[NBLOCKS]; /**< RMS of the reflection coefficients */
427     int energy = 0;
428     int i, idx;
429
430     if (buf_size < FRAMESIZE) {
431         av_log(avctx, AV_LOG_ERROR, "output buffer too small\n");
432         return 0;
433     }
434     ractx = avctx->priv_data;
435
436     /**
437      * Since the LPC coefficients are calculated on a frame centered over the
438      * fourth subframe, to encode a given frame, data from the next frame is
439      * needed. In each call to this function, the previous frame (whose data are
440      * saved in the encoder context) is encoded, and data from the current frame
441      * are saved in the encoder context to be used in the next function call.
442      */
443     for (i = 0; i < (2 * BLOCKSIZE + BLOCKSIZE / 2); i++) {
444         lpc_data[i] = ractx->curr_block[BLOCKSIZE + BLOCKSIZE / 2 + i];
445         energy += (lpc_data[i] * lpc_data[i]) >> 4;
446     }
447     for (i = 2 * BLOCKSIZE + BLOCKSIZE / 2; i < NBLOCKS * BLOCKSIZE; i++) {
448         lpc_data[i] = *((int16_t *)data + i - 2 * BLOCKSIZE - BLOCKSIZE / 2) >>
449                       2;
450         energy += (lpc_data[i] * lpc_data[i]) >> 4;
451     }
452     energy = ff_energy_tab[quantize(ff_t_sqrt(energy >> 5) >> 10, ff_energy_tab,
453                                     32)];
454
455     ff_lpc_calc_coefs(&ractx->dsp, lpc_data, NBLOCKS * BLOCKSIZE, LPC_ORDER,
456                       LPC_ORDER, 16, lpc_coefs, shift, 1, ORDER_METHOD_EST, 12,
457                       0);
458     for (i = 0; i < LPC_ORDER; i++)
459         block_coefs[NBLOCKS - 1][i] = -(lpc_coefs[LPC_ORDER - 1][i] <<
460                                         (12 - shift[LPC_ORDER - 1]));
461
462     /**
463      * TODO: apply perceptual weighting of the input speech through bandwidth
464      * expansion of the LPC filter.
465      */
466
467     if (ff_eval_refl(lpc_refl, block_coefs[NBLOCKS - 1], avctx)) {
468         /**
469          * The filter is unstable: use the coefficients of the previous frame.
470          */
471         ff_int_to_int16(block_coefs[NBLOCKS - 1], ractx->lpc_coef[1]);
472         ff_eval_refl(lpc_refl, block_coefs[NBLOCKS - 1], avctx);
473     }
474     init_put_bits(&pb, frame, buf_size);
475     for (i = 0; i < LPC_ORDER; i++) {
476         idx = quantize(lpc_refl[i], ff_lpc_refl_cb[i], sizes[i]);
477         put_bits(&pb, bit_sizes[i], idx);
478         lpc_refl[i] = ff_lpc_refl_cb[i][idx];
479     }
480     ractx->lpc_refl_rms[0] = ff_rms(lpc_refl);
481     ff_eval_coefs(ractx->lpc_coef[0], lpc_refl);
482     refl_rms[0] = ff_interp(ractx, block_coefs[0], 1, 1, ractx->old_energy);
483     refl_rms[1] = ff_interp(ractx, block_coefs[1], 2,
484                             energy <= ractx->old_energy,
485                             ff_t_sqrt(energy * ractx->old_energy) >> 12);
486     refl_rms[2] = ff_interp(ractx, block_coefs[2], 3, 0, energy);
487     refl_rms[3] = ff_rescale_rms(ractx->lpc_refl_rms[0], energy);
488     ff_int_to_int16(block_coefs[NBLOCKS - 1], ractx->lpc_coef[0]);
489     put_bits(&pb, 5, quantize(energy, ff_energy_tab, 32));
490     for (i = 0; i < NBLOCKS; i++)
491         ra144_encode_subblock(ractx, ractx->curr_block + i * BLOCKSIZE,
492                               block_coefs[i], refl_rms[i], &pb);
493     flush_put_bits(&pb);
494     ractx->old_energy = energy;
495     ractx->lpc_refl_rms[1] = ractx->lpc_refl_rms[0];
496     FFSWAP(unsigned int *, ractx->lpc_coef[0], ractx->lpc_coef[1]);
497     for (i = 0; i < NBLOCKS * BLOCKSIZE; i++)
498         ractx->curr_block[i] = *((int16_t *)data + i) >> 2;
499     return FRAMESIZE;
500 }
501
502
503 AVCodec ra_144_encoder =
504 {
505     "real_144",
506     CODEC_TYPE_AUDIO,
507     CODEC_ID_RA_144,
508     sizeof(RA144Context),
509     ra144_encode_init,
510     ra144_encode_frame,
511     .long_name = NULL_IF_CONFIG_SMALL("RealAudio 1.0 (14.4K) encoder"),
512 };