LCOV - code coverage report
Current view: top level - libavcodec - mdct15.c (source / functions) Hit Total Coverage
Test: coverage.info Lines: 145 180 80.6 %
Date: 2017-12-16 01:21:47 Functions: 7 8 87.5 %

          Line data    Source code
       1             : /*
       2             :  * Copyright (c) 2013-2014 Mozilla Corporation
       3             :  * Copyright (c) 2017 Rostislav Pehlivanov <atomnuker@gmail.com>
       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
      24             :  * Celt non-power of 2 iMDCT
      25             :  */
      26             : 
      27             : #include <float.h>
      28             : #include <math.h>
      29             : #include <stddef.h>
      30             : 
      31             : #include "config.h"
      32             : 
      33             : #include "libavutil/attributes.h"
      34             : #include "libavutil/common.h"
      35             : 
      36             : #include "mdct15.h"
      37             : 
      38             : #define FFT_FLOAT 1
      39             : #include "fft-internal.h"
      40             : 
      41             : #define CMUL3(c, a, b) CMUL((c).re, (c).im, (a).re, (a).im, (b).re, (b).im)
      42             : 
      43         902 : av_cold void ff_mdct15_uninit(MDCT15Context **ps)
      44             : {
      45         902 :     MDCT15Context *s = *ps;
      46             : 
      47         902 :     if (!s)
      48           0 :         return;
      49             : 
      50         902 :     ff_fft_end(&s->ptwo_fft);
      51             : 
      52         902 :     av_freep(&s->pfa_prereindex);
      53         902 :     av_freep(&s->pfa_postreindex);
      54         902 :     av_freep(&s->twiddle_exptab);
      55         902 :     av_freep(&s->tmp);
      56             : 
      57         902 :     av_freep(ps);
      58             : }
      59             : 
      60         902 : static inline int init_pfa_reindex_tabs(MDCT15Context *s)
      61             : {
      62             :     int i, j;
      63         902 :     const int b_ptwo = s->ptwo_fft.nbits; /* Bits for the power of two FFTs */
      64         902 :     const int l_ptwo = 1 << b_ptwo; /* Total length for the power of two FFTs */
      65         902 :     const int inv_1 = l_ptwo << ((4 - b_ptwo) & 3); /* (2^b_ptwo)^-1 mod 15 */
      66         902 :     const int inv_2 = 0xeeeeeeef & ((1U << b_ptwo) - 1); /* 15^-1 mod 2^b_ptwo */
      67             : 
      68         902 :     s->pfa_prereindex = av_malloc_array(15 * l_ptwo, sizeof(*s->pfa_prereindex));
      69         902 :     if (!s->pfa_prereindex)
      70           0 :         return 1;
      71             : 
      72         902 :     s->pfa_postreindex = av_malloc_array(15 * l_ptwo, sizeof(*s->pfa_postreindex));
      73         902 :     if (!s->pfa_postreindex)
      74           0 :         return 1;
      75             : 
      76             :     /* Pre/Post-reindex */
      77       16154 :     for (i = 0; i < l_ptwo; i++) {
      78      244032 :         for (j = 0; j < 15; j++) {
      79      228780 :             const int q_pre = ((l_ptwo * j)/15 + i) >> b_ptwo;
      80      228780 :             const int q_post = (((j*inv_1)/15) + (i*inv_2)) >> b_ptwo;
      81      228780 :             const int k_pre = 15*i + (j - q_pre*15)*(1 << b_ptwo);
      82      228780 :             const int k_post = i*inv_2*15 + j*inv_1 - 15*q_post*l_ptwo;
      83      228780 :             s->pfa_prereindex[i*15 + j] = k_pre << 1;
      84      228780 :             s->pfa_postreindex[k_post] = l_ptwo*j + i;
      85             :         }
      86             :     }
      87             : 
      88         902 :     return 0;
      89             : }
      90             : 
      91             : /* Stride is hardcoded to 3 */
      92     2662464 : static inline void fft5(FFTComplex *out, FFTComplex *in, FFTComplex exptab[2])
      93             : {
      94             :     FFTComplex z0[4], t[6];
      95             : 
      96     2662464 :     t[0].re = in[3].re + in[12].re;
      97     2662464 :     t[0].im = in[3].im + in[12].im;
      98     2662464 :     t[1].im = in[3].re - in[12].re;
      99     2662464 :     t[1].re = in[3].im - in[12].im;
     100     2662464 :     t[2].re = in[6].re + in[ 9].re;
     101     2662464 :     t[2].im = in[6].im + in[ 9].im;
     102     2662464 :     t[3].im = in[6].re - in[ 9].re;
     103     2662464 :     t[3].re = in[6].im - in[ 9].im;
     104             : 
     105     2662464 :     out[0].re = in[0].re + in[3].re + in[6].re + in[9].re + in[12].re;
     106     2662464 :     out[0].im = in[0].im + in[3].im + in[6].im + in[9].im + in[12].im;
     107             : 
     108     2662464 :     t[4].re = exptab[0].re * t[2].re - exptab[1].re * t[0].re;
     109     2662464 :     t[4].im = exptab[0].re * t[2].im - exptab[1].re * t[0].im;
     110     2662464 :     t[0].re = exptab[0].re * t[0].re - exptab[1].re * t[2].re;
     111     2662464 :     t[0].im = exptab[0].re * t[0].im - exptab[1].re * t[2].im;
     112     2662464 :     t[5].re = exptab[0].im * t[3].re - exptab[1].im * t[1].re;
     113     2662464 :     t[5].im = exptab[0].im * t[3].im - exptab[1].im * t[1].im;
     114     2662464 :     t[1].re = exptab[0].im * t[1].re + exptab[1].im * t[3].re;
     115     2662464 :     t[1].im = exptab[0].im * t[1].im + exptab[1].im * t[3].im;
     116             : 
     117     2662464 :     z0[0].re = t[0].re - t[1].re;
     118     2662464 :     z0[0].im = t[0].im - t[1].im;
     119     2662464 :     z0[1].re = t[4].re + t[5].re;
     120     2662464 :     z0[1].im = t[4].im + t[5].im;
     121             : 
     122     2662464 :     z0[2].re = t[4].re - t[5].re;
     123     2662464 :     z0[2].im = t[4].im - t[5].im;
     124     2662464 :     z0[3].re = t[0].re + t[1].re;
     125     2662464 :     z0[3].im = t[0].im + t[1].im;
     126             : 
     127     2662464 :     out[1].re = in[0].re + z0[3].re;
     128     2662464 :     out[1].im = in[0].im + z0[0].im;
     129     2662464 :     out[2].re = in[0].re + z0[2].re;
     130     2662464 :     out[2].im = in[0].im + z0[1].im;
     131     2662464 :     out[3].re = in[0].re + z0[1].re;
     132     2662464 :     out[3].im = in[0].im + z0[2].im;
     133     2662464 :     out[4].re = in[0].re + z0[0].re;
     134     2662464 :     out[4].im = in[0].im + z0[3].im;
     135     2662464 : }
     136             : 
     137      887488 : static void fft15_c(FFTComplex *out, FFTComplex *in, FFTComplex *exptab, ptrdiff_t stride)
     138             : {
     139             :     int k;
     140             :     FFTComplex tmp1[5], tmp2[5], tmp3[5];
     141             : 
     142      887488 :     fft5(tmp1, in + 0, exptab + 19);
     143      887488 :     fft5(tmp2, in + 1, exptab + 19);
     144      887488 :     fft5(tmp3, in + 2, exptab + 19);
     145             : 
     146     5324928 :     for (k = 0; k < 5; k++) {
     147             :         FFTComplex t[2];
     148             : 
     149     4437440 :         CMUL3(t[0], tmp2[k], exptab[k]);
     150     4437440 :         CMUL3(t[1], tmp3[k], exptab[2 * k]);
     151     4437440 :         out[stride*k].re = tmp1[k].re + t[0].re + t[1].re;
     152     4437440 :         out[stride*k].im = tmp1[k].im + t[0].im + t[1].im;
     153             : 
     154     4437440 :         CMUL3(t[0], tmp2[k], exptab[k + 5]);
     155     4437440 :         CMUL3(t[1], tmp3[k], exptab[2 * (k + 5)]);
     156     4437440 :         out[stride*(k + 5)].re = tmp1[k].re + t[0].re + t[1].re;
     157     4437440 :         out[stride*(k + 5)].im = tmp1[k].im + t[0].im + t[1].im;
     158             : 
     159     4437440 :         CMUL3(t[0], tmp2[k], exptab[k + 10]);
     160     4437440 :         CMUL3(t[1], tmp3[k], exptab[2 * k + 5]);
     161     4437440 :         out[stride*(k + 10)].re = tmp1[k].re + t[0].re + t[1].re;
     162     4437440 :         out[stride*(k + 10)].im = tmp1[k].im + t[0].im + t[1].im;
     163             :     }
     164      887488 : }
     165             : 
     166           0 : static void mdct15(MDCT15Context *s, float *dst, const float *src, ptrdiff_t stride)
     167             : {
     168             :     int i, j;
     169           0 :     const int len4 = s->len4, len3 = len4 * 3, len8 = len4 >> 1;
     170           0 :     const int l_ptwo = 1 << s->ptwo_fft.nbits;
     171             :     FFTComplex fft15in[15];
     172             : 
     173             :     /* Folding and pre-reindexing */
     174           0 :     for (i = 0; i < l_ptwo; i++) {
     175           0 :         for (j = 0; j < 15; j++) {
     176           0 :             const int k = s->pfa_prereindex[i*15 + j];
     177           0 :             FFTComplex tmp, exp = s->twiddle_exptab[k >> 1];
     178           0 :             if (k < len4) {
     179           0 :                 tmp.re = -src[ len4 + k] + src[1*len4 - 1 - k];
     180           0 :                 tmp.im = -src[ len3 + k] - src[1*len3 - 1 - k];
     181             :             } else {
     182           0 :                 tmp.re = -src[ len4 + k] - src[5*len4 - 1 - k];
     183           0 :                 tmp.im =  src[-len4 + k] - src[1*len3 - 1 - k];
     184             :             }
     185           0 :             CMUL(fft15in[j].im, fft15in[j].re, tmp.re, tmp.im, exp.re, exp.im);
     186             :         }
     187           0 :         s->fft15(s->tmp + s->ptwo_fft.revtab[i], fft15in, s->exptab, l_ptwo);
     188             :     }
     189             : 
     190             :     /* Then a 15xN FFT (where N is a power of two) */
     191           0 :     for (i = 0; i < 15; i++)
     192           0 :         s->ptwo_fft.fft_calc(&s->ptwo_fft, s->tmp + l_ptwo*i);
     193             : 
     194             :     /* Reindex again, apply twiddles and output */
     195           0 :     for (i = 0; i < len8; i++) {
     196           0 :         const int i0 = len8 + i, i1 = len8 - i - 1;
     197           0 :         const int s0 = s->pfa_postreindex[i0], s1 = s->pfa_postreindex[i1];
     198             : 
     199           0 :         CMUL(dst[2*i1*stride + stride], dst[2*i0*stride], s->tmp[s0].re, s->tmp[s0].im,
     200             :              s->twiddle_exptab[i0].im, s->twiddle_exptab[i0].re);
     201           0 :         CMUL(dst[2*i0*stride + stride], dst[2*i1*stride], s->tmp[s1].re, s->tmp[s1].im,
     202             :              s->twiddle_exptab[i1].im, s->twiddle_exptab[i1].re);
     203             :     }
     204           0 : }
     205             : 
     206       99968 : static void imdct15_half(MDCT15Context *s, float *dst, const float *src,
     207             :                          ptrdiff_t stride)
     208             : {
     209             :     FFTComplex fft15in[15];
     210       99968 :     FFTComplex *z = (FFTComplex *)dst;
     211       99968 :     int i, j, len8 = s->len4 >> 1, l_ptwo = 1 << s->ptwo_fft.nbits;
     212       99968 :     const float *in1 = src, *in2 = src + (s->len2 - 1) * stride;
     213             : 
     214             :     /* Reindex input, putting it into a buffer and doing an Nx15 FFT */
     215      987456 :     for (i = 0; i < l_ptwo; i++) {
     216    14199808 :         for (j = 0; j < 15; j++) {
     217    13312320 :             const int k = s->pfa_prereindex[i*15 + j];
     218    13312320 :             FFTComplex tmp = { in2[-k*stride], in1[k*stride] };
     219    13312320 :             CMUL3(fft15in[j], tmp, s->twiddle_exptab[k >> 1]);
     220             :         }
     221      887488 :         s->fft15(s->tmp + s->ptwo_fft.revtab[i], fft15in, s->exptab, l_ptwo);
     222             :     }
     223             : 
     224             :     /* Then a 15xN FFT (where N is a power of two) */
     225     1599488 :     for (i = 0; i < 15; i++)
     226     1499520 :         s->ptwo_fft.fft_calc(&s->ptwo_fft, s->tmp + l_ptwo*i);
     227             : 
     228             :     /* Reindex again, apply twiddles and output */
     229       99968 :     s->postreindex(z, s->tmp, s->twiddle_exptab, s->pfa_postreindex, len8);
     230       99968 : }
     231             : 
     232       99968 : static void postrotate_c(FFTComplex *out, FFTComplex *in, FFTComplex *exp,
     233             :                          int *lut, ptrdiff_t len8)
     234             : {
     235             :     int i;
     236             : 
     237             :     /* Reindex again, apply twiddles and output */
     238     6756128 :     for (i = 0; i < len8; i++) {
     239     6656160 :         const int i0 = len8 + i, i1 = len8 - i - 1;
     240     6656160 :         const int s0 = lut[i0], s1 = lut[i1];
     241             : 
     242     6656160 :         CMUL(out[i1].re, out[i0].im, in[s1].im, in[s1].re, exp[i1].im, exp[i1].re);
     243     6656160 :         CMUL(out[i0].re, out[i1].im, in[s0].im, in[s0].re, exp[i0].im, exp[i0].re);
     244             :     }
     245       99968 : }
     246             : 
     247         902 : av_cold int ff_mdct15_init(MDCT15Context **ps, int inverse, int N, double scale)
     248             : {
     249             :     MDCT15Context *s;
     250             :     double alpha, theta;
     251         902 :     int len2 = 15 * (1 << N);
     252         902 :     int len  = 2 * len2;
     253             :     int i;
     254             : 
     255             :     /* Tested and verified to work on everything in between */
     256         902 :     if ((N < 2) || (N > 13))
     257           0 :         return AVERROR(EINVAL);
     258             : 
     259         902 :     s = av_mallocz(sizeof(*s));
     260         902 :     if (!s)
     261           0 :         return AVERROR(ENOMEM);
     262             : 
     263         902 :     s->fft_n       = N - 1;
     264         902 :     s->len4        = len2 / 2;
     265         902 :     s->len2        = len2;
     266         902 :     s->inverse     = inverse;
     267         902 :     s->fft15       = fft15_c;
     268         902 :     s->mdct        = mdct15;
     269         902 :     s->imdct_half  = imdct15_half;
     270         902 :     s->postreindex = postrotate_c;
     271             : 
     272         902 :     if (ff_fft_init(&s->ptwo_fft, N - 1, s->inverse) < 0)
     273           0 :         goto fail;
     274             : 
     275         902 :     if (init_pfa_reindex_tabs(s))
     276           0 :         goto fail;
     277             : 
     278         902 :     s->tmp  = av_malloc_array(len, 2 * sizeof(*s->tmp));
     279         902 :     if (!s->tmp)
     280           0 :         goto fail;
     281             : 
     282         902 :     s->twiddle_exptab = av_malloc_array(s->len4, sizeof(*s->twiddle_exptab));
     283         902 :     if (!s->twiddle_exptab)
     284           0 :         goto fail;
     285             : 
     286         902 :     theta = 0.125f + (scale < 0 ? s->len4 : 0);
     287         902 :     scale = sqrt(fabs(scale));
     288      229682 :     for (i = 0; i < s->len4; i++) {
     289      228780 :         alpha = 2 * M_PI * (i + theta) / len;
     290      228780 :         s->twiddle_exptab[i].re = cosf(alpha) * scale;
     291      228780 :         s->twiddle_exptab[i].im = sinf(alpha) * scale;
     292             :     }
     293             : 
     294             :     /* 15-point FFT exptab */
     295       18040 :     for (i = 0; i < 19; i++) {
     296       17138 :         if (i < 15) {
     297       13530 :             double theta = (2.0f * M_PI * i) / 15.0f;
     298       13530 :             if (!s->inverse)
     299           0 :                 theta *= -1;
     300       13530 :             s->exptab[i].re = cosf(theta);
     301       13530 :             s->exptab[i].im = sinf(theta);
     302             :         } else { /* Wrap around to simplify fft15 */
     303        3608 :             s->exptab[i] = s->exptab[i - 15];
     304             :         }
     305             :     }
     306             : 
     307             :     /* 5-point FFT exptab */
     308         902 :     s->exptab[19].re = cosf(2.0f * M_PI / 5.0f);
     309         902 :     s->exptab[19].im = sinf(2.0f * M_PI / 5.0f);
     310         902 :     s->exptab[20].re = cosf(1.0f * M_PI / 5.0f);
     311         902 :     s->exptab[20].im = sinf(1.0f * M_PI / 5.0f);
     312             : 
     313             :     /* Invert the phase for an inverse transform, do nothing for a forward transform */
     314         902 :     if (s->inverse) {
     315         902 :         s->exptab[19].im *= -1;
     316         902 :         s->exptab[20].im *= -1;
     317             :     }
     318             : 
     319             :     if (ARCH_X86)
     320         902 :         ff_mdct15_init_x86(s);
     321             : 
     322         902 :     *ps = s;
     323             : 
     324         902 :     return 0;
     325             : 
     326           0 : fail:
     327           0 :     ff_mdct15_uninit(&s);
     328           0 :     return AVERROR(ENOMEM);
     329             : }

Generated by: LCOV version 1.13