X-Git-Url: https://git.sesse.net/?a=blobdiff_plain;f=libavcodec%2Fra288.c;h=abc7df51bd65f7c620c910b74141ab3adc8a0ff9;hb=ca6e7708b42e7d33ba3053bcd447d52a077bca25;hp=d68498746275c78625e9a490f8f1b44ea04d6406;hpb=c1fadf50129751812c1011e3e8e6a51044e9ea40;p=ffmpeg

diff --git a/libavcodec/ra288.c b/libavcodec/ra288.c
index d6849874627..abc7df51bd6 100644
--- a/libavcodec/ra288.c
+++ b/libavcodec/ra288.c
@@ -20,263 +20,186 @@
  */
 
 #include "avcodec.h"
+#define ALT_BITSTREAM_READER_LE
+#include "get_bits.h"
 #include "ra288.h"
+#include "lpc.h"
+#include "celp_math.h"
+#include "celp_filters.h"
 
 typedef struct {
-    float history[8];
-    float output[40];
-    float pr1[36];
-    float pr2[10];
-    int   phase, phasep;
-
-    float st1a[111], st1b[37], st1[37];
-    float st2a[38], st2b[11], st2[11];
-    float sb[41];
-    float lhist[10];
-} Real288_internal;
-
-static void prodsum(float *tgt, float *src, int len, int n);
-static void co(int n, int i, int j, float *in, float *out, float *st1, float *st2, const float *table);
-static int pred(float *in, float *tgt, int n);
-static void colmult(float *tgt, float *m1, const float *m2, int n);
-
-
-/* initial decode */
-static void unpack(unsigned short *tgt, const unsigned char *src,
-                   unsigned int len)
-{
-    int x, y, z;
-    int n, temp;
-    int buffer[len];
+    float sp_lpc[36];      ///< LPC coefficients for speech data (spec: A)
+    float gain_lpc[10];    ///< LPC coefficients for gain        (spec: GB)
 
-    for (x=0; x < len; tgt[x++] = 0)
-        buffer[x] = 9 + (x & 1);
+    /** speech data history                                      (spec: SB).
+     *  Its first 70 coefficients are updated only at backward filtering.
+     */
+    float sp_hist[111];
 
-    for (x=y=z=0; x < len/*was 38*/; x++) {
-        n = buffer[y] - z;
-        temp = src[x];
+    /// speech part of the gain autocorrelation                  (spec: REXP)
+    float sp_rec[37];
 
-        if (n < 8)
-            temp &= 255 >> (8 - n);
+    /** log-gain history                                         (spec: SBLG).
+     *  Its first 28 coefficients are updated only at backward filtering.
+     */
+    float gain_hist[38];
 
-        tgt[y] += temp << z;
+    /// recursive part of the gain autocorrelation               (spec: REXPLG)
+    float gain_rec[11];
+} RA288Context;
 
-        if (n <= 8) {
-            tgt[++y] += src[x] >> n;
-            z = 8 - n;
-        } else
-            z += 8;
-  }
+static av_cold int ra288_decode_init(AVCodecContext *avctx)
+{
+    avctx->sample_fmt = SAMPLE_FMT_FLT;
+    return 0;
 }
 
-static void update(Real288_internal *glob)
+static void apply_window(float *tgt, const float *m1, const float *m2, int n)
 {
-    int x,y;
-    float buffer1[40], temp1[37];
-    float buffer2[8], temp2[11];
-
-    for (x=0, y=glob->phasep+5; x < 40; buffer1[x++] = glob->output[(y++)%40]);
-
-    co(36, 40, 35, buffer1, temp1, glob->st1a, glob->st1b, table1);
-
-    if (pred(temp1, glob->st1, 36))
-        colmult(glob->pr1, glob->st1, table1a, 36);
-
-    for (x=0, y=glob->phase + 1; x < 8; buffer2[x++] = glob->history[(y++) % 8]);
-
-    co(10, 8, 20, buffer2, temp2, glob->st2a, glob->st2b, table2);
-
-    if (pred(temp2, glob->st2, 10))
-        colmult(glob->pr2, glob->st2, table2a, 10);
+    while (n--)
+        *tgt++ = *m1++ * *m2++;
 }
 
-/* Decode and produce output */
-static void decode(Real288_internal *glob, unsigned int input)
+static void convolve(float *tgt, const float *src, int len, int n)
 {
-    unsigned int x, y;
-    float f;
-    double sum, sumsum;
-    float *p1, *p2;
-    float buffer[5];
-    const float *table;
+    for (; n >= 0; n--)
+        tgt[n] = ff_dot_productf(src, src - n, len);
 
-    for (x=36; x--; glob->sb[x+5] = glob->sb[x]);
+}
 
-    for (x=5; x--;) {
-        p1 = glob->sb+x;
-        p2 = glob->pr1;
-        for (sum=0, y=36; y--; sum -= (*(++p1))*(*(p2++)));
+static void decode(RA288Context *ractx, float gain, int cb_coef)
+{
+    int i;
+    double sumsum;
+    float sum, buffer[5];
+    float *block = ractx->sp_hist + 70 + 36; // current block
+    float *gain_block = ractx->gain_hist + 28;
 
-        glob->sb[x] = sum;
-    }
+    memmove(ractx->sp_hist + 70, ractx->sp_hist + 75, 36*sizeof(*block));
 
-    f = amptable[input & 7];
-    table = codetable + (input >> 3) * 5;
+    /* block 46 of G.728 spec */
+    sum = 32.;
+    for (i=0; i < 10; i++)
+        sum -= gain_block[9-i] * ractx->gain_lpc[i];
 
-    /* convert log and do rms */
-    for (sum=32, x=10; x--; sum -= glob->pr2[x] * glob->lhist[x]);
+    /* block 47 of G.728 spec */
+    sum = av_clipf(sum, 0, 60);
 
-    if (sum < 0)
-        sum = 0;
-    else if (sum > 60)
-        sum = 60;
+    /* block 48 of G.728 spec */
+    /* exp(sum * 0.1151292546497) == pow(10.0,sum/20) */
+    sumsum = exp(sum * 0.1151292546497) * gain * (1.0/(1<<23));
 
-    sumsum = exp(sum * 0.1151292546497) * f;    /* pow(10.0,sum/20)*f */
+    for (i=0; i < 5; i++)
+        buffer[i] = codetable[cb_coef][i] * sumsum;
 
-    for (sum=0, x=5; x--;) {
-        buffer[x] = table[x] * sumsum;
-        sum += buffer[x] * buffer[x];
-    }
+    sum = ff_dot_productf(buffer, buffer, 5) * ((1<<24)/5.);
 
-    if ((sum /= 5) < 1)
-        sum = 1;
+    sum = FFMAX(sum, 1);
 
     /* shift and store */
-    for (x=10; --x; glob->lhist[x] = glob->lhist[x-1]);
+    memmove(gain_block, gain_block + 1, 9 * sizeof(*gain_block));
 
-    *glob->lhist = glob->history[glob->phase] = 10 * log10(sum) - 32;
+    gain_block[9] = 10 * log10(sum) - 32;
 
-    for (x=1; x < 5; x++)
-        for (y=x; y--; buffer[x] -= glob->pr1[x-y-1] * buffer[y]);
+    ff_celp_lp_synthesis_filterf(block, ractx->sp_lpc, buffer, 5, 36);
 
     /* output */
-    for (x=0; x < 5; x++) {
-        f = glob->sb[4-x] + buffer[x];
-
-        if (f > 4095)
-            f = 4095;
-        else if (f < -4095)
-            f = -4095;
-
-        glob->output[glob->phasep+x] = glob->sb[4-x] = f;
-    }
-}
-
-/* column multiply */
-static void colmult(float *tgt, float *m1, const float *m2, int n)
-{
-    while (n--)
-        *(tgt++) = (*(m1++)) * (*(m2++));
+    for (i=0; i < 5; i++)
+        block[i] = av_clipf(block[i], -4095./4096., 4095./4096.);
 }
 
-static int pred(float *in, float *tgt, int n)
+/**
+ * Hybrid window filtering, see blocks 36 and 49 of the G.728 specification.
+ *
+ * @param order   filter order
+ * @param n       input length
+ * @param non_rec number of non-recursive samples
+ * @param out     filter output
+ * @param hist    pointer to the input history of the filter
+ * @param out     pointer to the non-recursive part of the output
+ * @param out2    pointer to the recursive part of the output
+ * @param window  pointer to the windowing function table
+ */
+static void do_hybrid_window(int order, int n, int non_rec, float *out,
+                             float *hist, float *out2, const float *window)
 {
-    int x, y;
-    float *p1, *p2;
-    double f0, f1, f2;
-    float temp;
-
-    if (in[n] == 0)
-        return 0;
+    int i;
+    float buffer1[order + 1];
+    float buffer2[order + 1];
+    float work[order + n + non_rec];
 
-    if ((f0 = *in) <= 0)
-        return 0;
-
-    for (x=1 ; ; x++) {
-        if (n < x)
-            return 1;
-
-        p1 = in + x;
-        p2 = tgt;
-        f1 = *(p1--);
-        for (y=x; --y; f1 += (*(p1--))*(*(p2++)));
-
-        p1 = tgt + x - 1;
-        p2 = tgt;
-        *(p1--) = f2 = -f1/f0;
-        for (y=x >> 1; y--;) {
-            temp = *p2 + *p1 * f2;
-            *(p1--) += *p2 * f2;
-            *(p2++) = temp;
-        }
-        if ((f0 += f1*f2) < 0)
-            return 0;
-    }
-}
-
-static void co(int n, int i, int j, float *in, float *out, float *st1,
-               float *st2, const float *table)
-{
-    int a, b, c;
-    unsigned int x;
-    float *fp;
-    float buffer1[37];
-    float buffer2[37];
-    float work[111];
-
-    /* rotate and multiply */
-    c = (b = (a = n + i) + j) - i;
-    fp = st1 + i;
-    for (x=0; x < b; x++) {
-        if (x == c)
-            fp=in;
-        work[x] = *(table++) * (*(st1++) = *(fp++));
-    }
+    apply_window(work, window, hist, order + n + non_rec);
 
-    prodsum(buffer1, work + n, i, n);
-    prodsum(buffer2, work + a, j, n);
+    convolve(buffer1, work + order    , n      , order);
+    convolve(buffer2, work + order + n, non_rec, order);
 
-    for (x=0;x<=n;x++) {
-        *st2 = *st2 * (0.5625) + buffer1[x];
-        out[x] = *(st2++) + buffer2[x];
+    for (i=0; i <= order; i++) {
+        out2[i] = out2[i] * 0.5625 + buffer1[i];
+        out [i] = out2[i]          + buffer2[i];
     }
-    *out *= 1.00390625; /* to prevent clipping */
-}
 
-/* product sum (lsf) */
-static void prodsum(float *tgt, float *src, int len, int n)
-{
-    unsigned int x;
-    float *p1, *p2;
-    double sum;
-
-    while (n >= 0) {
-        p1 = (p2 = src) - n;
-        for (sum=0, x=len; x--; sum += (*p1++) * (*p2++));
-        tgt[n--] = sum;
-    }
+    /* Multiply by the white noise correcting factor (WNCF). */
+    *out *= 257./256.;
 }
 
-static void * decode_block(AVCodecContext * avctx, const unsigned char *in,
-                           signed short int *out, unsigned len)
+/**
+ * Backward synthesis filter, find the LPC coefficients from past speech data.
+ */
+static void backward_filter(float *hist, float *rec, const float *window,
+                            float *lpc, const float *tab,
+                            int order, int n, int non_rec, int move_size)
 {
-    int x, y;
-    Real288_internal *glob = avctx->priv_data;
-    unsigned short int buffer[len];
+    float temp[order+1];
 
-    unpack(buffer, in, len);
+    do_hybrid_window(order, n, non_rec, temp, hist, rec, window);
 
-    for (x=0; x < 32; x++) {
-        glob->phasep = (glob->phase = x & 7) * 5;
-        decode(glob, buffer[x]);
+    if (!compute_lpc_coefs(temp, order, lpc, 0, 1, 1))
+        apply_window(lpc, lpc, tab, order);
 
-        for (y=0; y<5; *(out++) = 8 * glob->output[glob->phasep+(y++)]);
-
-        if (glob->phase == 3)
-            update(glob);
-    }
-
-    return out;
+    memmove(hist, hist + n, move_size*sizeof(*hist));
 }
 
-/* Decode a block (celp) */
 static int ra288_decode_frame(AVCodecContext * avctx, void *data,
-                              int *data_size, const uint8_t * buf,
-                              int buf_size)
+                              int *data_size, AVPacket *avpkt)
 {
-    void *datao;
+    const uint8_t *buf = avpkt->data;
+    int buf_size = avpkt->size;
+    float *out = data;
+    int i, j;
+    RA288Context *ractx = avctx->priv_data;
+    GetBitContext gb;
 
     if (buf_size < avctx->block_align) {
         av_log(avctx, AV_LOG_ERROR,
-               "ffra288: Error! Input buffer is too small [%d<%d]\n",
+               "Error! Input buffer is too small [%d<%d]\n",
                buf_size, avctx->block_align);
         return 0;
     }
 
-    datao = data;
-    data = decode_block(avctx, buf, (signed short *)data, avctx->block_align);
+    if (*data_size < 32*5*4)
+        return -1;
+
+    init_get_bits(&gb, buf, avctx->block_align * 8);
+
+    for (i=0; i < 32; i++) {
+        float gain = amptable[get_bits(&gb, 3)];
+        int cb_coef = get_bits(&gb, 6 + (i&1));
+
+        decode(ractx, gain, cb_coef);
 
-    *data_size = (char *)data - (char *)datao;
+        for (j=0; j < 5; j++)
+            *(out++) = ractx->sp_hist[70 + 36 + j];
+
+        if ((i & 7) == 3) {
+            backward_filter(ractx->sp_hist, ractx->sp_rec, syn_window,
+                            ractx->sp_lpc, syn_bw_tab, 36, 40, 35, 70);
+
+            backward_filter(ractx->gain_hist, ractx->gain_rec, gain_window,
+                            ractx->gain_lpc, gain_bw_tab, 10, 8, 20, 28);
+        }
+    }
+
+    *data_size = (char *)out - (char *)data;
     return avctx->block_align;
 }
 
@@ -285,8 +208,8 @@ AVCodec ra_288_decoder =
     "real_288",
     CODEC_TYPE_AUDIO,
     CODEC_ID_RA_288,
-    sizeof(Real288_internal),
-    NULL,
+    sizeof(RA288Context),
+    ra288_decode_init,
     NULL,
     NULL,
     ra288_decode_frame,