git.sesse.net Git - ffmpeg/blob - libavutil/tx.c

   1 /*
   2  * This file is part of FFmpeg.
   3  *
   4  * FFmpeg is free software; you can redistribute it and/or
   5  * modify it under the terms of the GNU Lesser General Public
   6  * License as published by the Free Software Foundation; either
   7  * version 2.1 of the License, or (at your option) any later version.
   8  *
   9  * FFmpeg is distributed in the hope that it will be useful,
  10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12  * Lesser General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU Lesser General Public
  15  * License along with FFmpeg; if not, write to the Free Software
  16  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
  17  */
  18
  19 #include "tx_priv.h"
  20
  21 int ff_tx_type_is_mdct(enum AVTXType type)
  22 {
  23     switch (type) {
  24     case AV_TX_FLOAT_MDCT:
  25     case AV_TX_DOUBLE_MDCT:
  26     case AV_TX_INT32_MDCT:
  27         return 1;
  28     default:
  29         return 0;
  30     }
  31 }
  32
  33 /* Calculates the modular multiplicative inverse */
  34 static av_always_inline int mulinv(int n, int m)
  35 {
  36     n = n % m;
  37     for (int x = 1; x < m; x++)
  38         if (((n * x) % m) == 1)
  39             return x;
  40     av_assert0(0); /* Never reached */
  41 }
  42
  43 /* Guaranteed to work for any n, m where gcd(n, m) == 1 */
  44 int ff_tx_gen_compound_mapping(AVTXContext *s)
  45 {
  46     int *in_map, *out_map;
  47     const int n     = s->n;
  48     const int m     = s->m;
  49     const int inv   = s->inv;
  50     const int len   = n*m;
  51     const int m_inv = mulinv(m, n);
  52     const int n_inv = mulinv(n, m);
  53     const int mdct  = ff_tx_type_is_mdct(s->type);
  54
  55     if (!(s->pfatab = av_malloc(2*len*sizeof(*s->pfatab))))
  56         return AVERROR(ENOMEM);
  57
  58     in_map  = s->pfatab;
  59     out_map = s->pfatab + n*m;
  60
  61     /* Ruritanian map for input, CRT map for output, can be swapped */
  62     for (int j = 0; j < m; j++) {
  63         for (int i = 0; i < n; i++) {
  64             /* Shifted by 1 to simplify MDCTs */
  65             in_map[j*n + i] = ((i*m + j*n) % len) << mdct;
  66             out_map[(i*m*m_inv + j*n*n_inv) % len] = i*m + j;
  67         }
  68     }
  69
  70     /* Change transform direction by reversing all ACs */
  71     if (inv) {
  72         for (int i = 0; i < m; i++) {
  73             int *in = &in_map[i*n + 1]; /* Skip the DC */
  74             for (int j = 0; j < ((n - 1) >> 1); j++)
  75                 FFSWAP(int, in[j], in[n - j - 2]);
  76         }
  77     }
  78
  79     /* Our 15-point transform is also a compound one, so embed its input map */
  80     if (n == 15) {
  81         for (int k = 0; k < m; k++) {
  82             int tmp[15];
  83             memcpy(tmp, &in_map[k*15], 15*sizeof(*tmp));
  84             for (int i = 0; i < 5; i++) {
  85                 for (int j = 0; j < 3; j++)
  86                     in_map[k*15 + i*3 + j] = tmp[(i*3 + j*5) % 15];
  87             }
  88         }
  89     }
  90
  91     return 0;
  92 }
  93
  94 static inline int split_radix_permutation(int i, int m, int inverse)
  95 {
  96     m >>= 1;
  97     if (m <= 1)
  98         return i & 1;
  99     if (!(i & m))
 100         return (split_radix_permutation(i, m, inverse) << 1);
 101     m >>= 1;
 102     return (split_radix_permutation(i, m, inverse) << 2) + 1 - 2*(!(i & m) ^ inverse);
 103 }
 104
 105 int ff_tx_gen_ptwo_revtab(AVTXContext *s, int invert_lookup)
 106 {
 107     const int m = s->m, inv = s->inv;
 108
 109     if (!(s->revtab = av_malloc(s->m*sizeof(*s->revtab))))
 110         return AVERROR(ENOMEM);
 111     if (!(s->revtab_c = av_malloc(m*sizeof(*s->revtab_c))))
 112         return AVERROR(ENOMEM);
 113
 114     /* Default */
 115     for (int i = 0; i < m; i++) {
 116         int k = -split_radix_permutation(i, m, inv) & (m - 1);
 117         if (invert_lookup)
 118             s->revtab[i] = s->revtab_c[i] = k;
 119         else
 120             s->revtab[i] = s->revtab_c[k] = i;
 121     }
 122
 123     return 0;
 124 }
 125
 126 int ff_tx_gen_ptwo_inplace_revtab_idx(AVTXContext *s, int *revtab)
 127 {
 128     int nb_inplace_idx = 0;
 129
 130     if (!(s->inplace_idx = av_malloc(s->m*sizeof(*s->inplace_idx))))
 131         return AVERROR(ENOMEM);
 132
 133     /* The first coefficient is always already in-place */
 134     for (int src = 1; src < s->m; src++) {
 135         int dst = revtab[src];
 136         int found = 0;
 137
 138         if (dst <= src)
 139             continue;
 140
 141         /* This just checks if a closed loop has been encountered before,
 142          * and if so, skips it, since to fully permute a loop we must only
 143          * enter it once. */
 144         do {
 145             for (int j = 0; j < nb_inplace_idx; j++) {
 146                 if (dst == s->inplace_idx[j]) {
 147                     found = 1;
 148                     break;
 149                 }
 150             }
 151             dst = revtab[dst];
 152         } while (dst != src && !found);
 153
 154         if (!found)
 155             s->inplace_idx[nb_inplace_idx++] = src;
 156     }
 157
 158     s->inplace_idx[nb_inplace_idx++] = 0;
 159
 160     return 0;
 161 }
 162
 163 static void parity_revtab_generator(int *revtab, int n, int inv, int offset,
 164                                     int is_dual, int dual_high, int len,
 165                                     int basis, int dual_stride)
 166 {
 167     len >>= 1;
 168
 169     if (len <= basis) {
 170         int k1, k2, *even, *odd, stride;
 171
 172         is_dual = is_dual && dual_stride;
 173         dual_high = is_dual & dual_high;
 174         stride = is_dual ? FFMIN(dual_stride, len) : 0;
 175
 176         even = &revtab[offset + dual_high*(stride - 2*len)];
 177         odd  = &even[len + (is_dual && !dual_high)*len + dual_high*len];
 178
 179         for (int i = 0; i < len; i++) {
 180             k1 = -split_radix_permutation(offset + i*2 + 0, n, inv) & (n - 1);
 181             k2 = -split_radix_permutation(offset + i*2 + 1, n, inv) & (n - 1);
 182             *even++ = k1;
 183             *odd++  = k2;
 184             if (stride && !((i + 1) % stride)) {
 185                 even += stride;
 186                 odd  += stride;
 187             }
 188         }
 189
 190         return;
 191     }
 192
 193     parity_revtab_generator(revtab, n, inv, offset,
 194                             0, 0, len >> 0, basis, dual_stride);
 195     parity_revtab_generator(revtab, n, inv, offset + (len >> 0),
 196                             1, 0, len >> 1, basis, dual_stride);
 197     parity_revtab_generator(revtab, n, inv, offset + (len >> 0) + (len >> 1),
 198                             1, 1, len >> 1, basis, dual_stride);
 199 }
 200
 201 void ff_tx_gen_split_radix_parity_revtab(int *revtab, int len, int inv,
 202                                          int basis, int dual_stride)
 203 {
 204     basis >>= 1;
 205     if (len < basis)
 206         return;
 207     av_assert0(!dual_stride || !(dual_stride & (dual_stride - 1)));
 208     av_assert0(dual_stride <= basis);
 209     parity_revtab_generator(revtab, len, inv, 0, 0, 0, len, basis, dual_stride);
 210 }
 211
 212 av_cold void av_tx_uninit(AVTXContext **ctx)
 213 {
 214     if (!(*ctx))
 215         return;
 216
 217     av_free((*ctx)->pfatab);
 218     av_free((*ctx)->exptab);
 219     av_free((*ctx)->revtab);
 220     av_free((*ctx)->revtab_c);
 221     av_free((*ctx)->inplace_idx);
 222     av_free((*ctx)->tmp);
 223
 224     av_freep(ctx);
 225 }
 226
 227 av_cold int av_tx_init(AVTXContext **ctx, av_tx_fn *tx, enum AVTXType type,
 228                        int inv, int len, const void *scale, uint64_t flags)
 229 {
 230     int err;
 231     AVTXContext *s = av_mallocz(sizeof(*s));
 232     if (!s)
 233         return AVERROR(ENOMEM);
 234
 235     switch (type) {
 236     case AV_TX_FLOAT_FFT:
 237     case AV_TX_FLOAT_MDCT:
 238         if ((err = ff_tx_init_mdct_fft_float(s, tx, type, inv, len, scale, flags)))
 239             goto fail;
 240         if (ARCH_X86)
 241             ff_tx_init_float_x86(s, tx);
 242         break;
 243     case AV_TX_DOUBLE_FFT:
 244     case AV_TX_DOUBLE_MDCT:
 245         if ((err = ff_tx_init_mdct_fft_double(s, tx, type, inv, len, scale, flags)))
 246             goto fail;
 247         break;
 248     case AV_TX_INT32_FFT:
 249     case AV_TX_INT32_MDCT:
 250         if ((err = ff_tx_init_mdct_fft_int32(s, tx, type, inv, len, scale, flags)))
 251             goto fail;
 252         break;
 253     default:
 254         err = AVERROR(EINVAL);
 255         goto fail;
 256     }
 257
 258     *ctx = s;
 259
 260     return 0;
 261
 262 fail:
 263     av_tx_uninit(&s);
 264     *tx = NULL;
 265     return err;
 266 }