src/third_party/openmax_dl/dl/sp/src/mips/mips_FFTInv_CCSToR_F32_complex.c

   1 /*
   2  *  Copyright (c) 2014 The WebRTC project authors. All Rights Reserved.
   3  *
   4  *  Use of this source code is governed by a BSD-style license
   5  *  that can be found in the LICENSE file in the root of the source
   6  *  tree. An additional intellectual property rights grant can be found
   7  *  in the file PATENTS.  All contributing project authors may
   8  *  be found in the AUTHORS file in the root of the source tree.
   9  *
  10  */
  11
  12 #include <stdint.h>
  13
  14 #include "dl/api/omxtypes.h"
  15 #include "dl/sp/api/mipsSP.h"
  16
  17 OMXResult mips_FFTInv_CCSToR_F32_complex(const OMX_F32* pSrc,
  18                                          OMX_F32* pDst,
  19                                          const MIPSFFTSpec_R_FC32* pFFTSpec) {
  20   OMX_U32 num_transforms, step;
  21   /* Quarter, half and three-quarters of transform size. */
  22   OMX_U32 n1_4, n1_2, n3_4;
  23   const OMX_F32* w_re_ptr;
  24   const OMX_F32* w_im_ptr;
  25   OMX_U32 fft_size = 1 << pFFTSpec->order;
  26   OMX_FC32* p_buf = (OMX_FC32*)pFFTSpec->pBuf;
  27   const OMX_FC32* p_src = (const OMX_FC32*)pSrc;
  28   OMX_FC32* p_dst;
  29   OMX_U16* p_bitrev = pFFTSpec->pBitRevInv;
  30   OMX_F32 tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7, tmp8, factor;
  31   OMX_F32 w_re, w_im;
  32
  33   w_re_ptr = pFFTSpec->pTwiddle + 1;
  34   w_im_ptr = pFFTSpec->pTwiddle + (OMX_U32)(1 << pFFTSpec->order - 2) - 1;
  35
  36   /*
  37    * Preliminary loop performing input adaptation due to computing real FFT
  38    * through complex FFT of half the size.
  39    */
  40   for (uint32_t n = 1; n < fft_size / 8; ++n) {
  41     tmp1 = p_src[n].Re;
  42     tmp2 = p_src[n].Im;
  43     tmp3 = p_src[fft_size / 2 - n].Re;
  44     tmp4 = p_src[fft_size / 2 - n].Im;
  45
  46     tmp5 = tmp1 + tmp3;
  47     tmp6 = tmp1 - tmp3;
  48     tmp7 = tmp2 + tmp4;
  49     tmp8 = tmp2 - tmp4;
  50
  51     w_re = *w_re_ptr;
  52     w_im = *w_im_ptr;
  53
  54     p_buf[p_bitrev[n]].Re = 0.5f * (tmp5 - w_re * tmp7 - w_im * tmp6);
  55     p_buf[p_bitrev[n]].Im = 0.5f * (tmp8 + w_re * tmp6 - w_im * tmp7);
  56     p_buf[p_bitrev[fft_size / 2 - n]].Re =
  57         0.5f * (tmp5 + w_re * tmp7 + w_im * tmp6);
  58     p_buf[p_bitrev[fft_size / 2 - n]].Im =
  59         0.5f * (-tmp8 + w_re * tmp6 - w_im * tmp7);
  60
  61     tmp1 = p_src[n + fft_size / 4].Re;
  62     tmp2 = p_src[n + fft_size / 4].Im;
  63     tmp3 = p_src[fft_size / 4 - n].Re;
  64     tmp4 = p_src[fft_size / 4 - n].Im;
  65
  66     tmp5 = tmp1 + tmp3;
  67     tmp6 = tmp1 - tmp3;
  68     tmp7 = tmp2 + tmp4;
  69     tmp8 = tmp2 - tmp4;
  70
  71     p_buf[p_bitrev[n + fft_size / 4]].Re =
  72         0.5f * (tmp5 + w_im * tmp7 - w_re * tmp6);
  73     p_buf[p_bitrev[n + fft_size / 4]].Im =
  74         0.5f * (tmp8 - w_im * tmp6 - w_re * tmp7);
  75     p_buf[p_bitrev[fft_size / 4 - n]].Re =
  76         0.5f * (tmp5 - w_im * tmp7 + w_re * tmp6);
  77     p_buf[p_bitrev[fft_size / 4 - n]].Im =
  78         0.5f * (-tmp8 - w_im * tmp6 - w_re * tmp7);
  79
  80     ++w_re_ptr;
  81     --w_im_ptr;
  82   }
  83   tmp1 = p_src[fft_size / 8].Re;
  84   tmp2 = p_src[fft_size / 8].Im;
  85   tmp3 = p_src[3 * fft_size / 8].Re;
  86   tmp4 = p_src[3 * fft_size / 8].Im;
  87
  88   tmp5 = tmp1 + tmp3;
  89   tmp6 = tmp1 - tmp3;
  90   tmp7 = tmp2 + tmp4;
  91   tmp8 = tmp2 - tmp4;
  92
  93   w_re = *w_re_ptr;
  94   w_im = *w_im_ptr;
  95
  96   p_buf[p_bitrev[fft_size / 8]].Re = 0.5f * (tmp5 - w_re * tmp7 - w_im * tmp6);
  97   p_buf[p_bitrev[fft_size / 8]].Im = 0.5f * (tmp8 + w_re * tmp6 - w_im * tmp7);
  98   p_buf[p_bitrev[3 * fft_size / 8]].Re =
  99       0.5f * (tmp5 + w_re * tmp7 + w_im * tmp6);
 100   p_buf[p_bitrev[3 * fft_size / 8]].Im =
 101       0.5f * (-tmp8 + w_re * tmp6 - w_im * tmp7);
 102
 103   p_buf[p_bitrev[0]].Re = 0.5f * (p_src[0].Re + p_src[fft_size / 2].Re);
 104   p_buf[p_bitrev[0]].Im = 0.5f * (p_src[0].Re - p_src[fft_size / 2].Re);
 105   p_buf[p_bitrev[fft_size / 4]].Re = p_src[fft_size / 4].Re;
 106   p_buf[p_bitrev[fft_size / 4]].Im = -p_src[fft_size / 4].Im;
 107
 108   /*
 109    * Loop performing sub-transforms of size 4,
 110    * which contain 2 butterfly operations.
 111    */
 112   num_transforms = (SUBTRANSFORM_CONST >> (17 - pFFTSpec->order)) | 1;
 113   for (uint32_t n = 0; n < num_transforms; ++n) {
 114     OMX_U32 offset = pFFTSpec->pOffset[n] << 2;
 115     OMX_FC32* p_tmp = p_buf + offset;
 116
 117     tmp1 = p_tmp[0].Re + p_tmp[1].Re;
 118     tmp2 = p_tmp[0].Im + p_tmp[1].Im;
 119     tmp3 = p_tmp[0].Re - p_tmp[1].Re;
 120     tmp4 = p_tmp[0].Im - p_tmp[1].Im;
 121     tmp5 = p_tmp[2].Re + p_tmp[3].Re;
 122     tmp6 = p_tmp[2].Im + p_tmp[3].Im;
 123     tmp8 = p_tmp[2].Im - p_tmp[3].Im;
 124     tmp7 = p_tmp[2].Re - p_tmp[3].Re;
 125
 126     p_tmp[0].Re = tmp1 + tmp5;
 127     p_tmp[2].Re = tmp1 - tmp5;
 128     p_tmp[0].Im = tmp2 + tmp6;
 129     p_tmp[2].Im = tmp2 - tmp6;
 130     p_tmp[1].Re = tmp3 - tmp8;
 131     p_tmp[3].Re = tmp3 + tmp8;
 132     p_tmp[1].Im = tmp4 + tmp7;
 133     p_tmp[3].Im = tmp4 - tmp7;
 134   }
 135
 136   num_transforms = (num_transforms >> 1) | 1;
 137   /*
 138    * Loop performing sub-transforms of size 8,
 139    * which contain four butterfly operations.
 140    */
 141   for (uint32_t n = 0; n < num_transforms; ++n) {
 142     OMX_U32 offset = pFFTSpec->pOffset[n] << 3;
 143     OMX_FC32* p_tmp = p_buf + offset;
 144
 145     tmp1 = p_tmp[4].Re + p_tmp[5].Re;
 146     tmp3 = p_tmp[6].Re + p_tmp[7].Re;
 147     tmp2 = p_tmp[4].Im + p_tmp[5].Im;
 148     tmp4 = p_tmp[6].Im + p_tmp[7].Im;
 149
 150     tmp5 = tmp1 + tmp3;
 151     tmp7 = tmp1 - tmp3;
 152     tmp6 = tmp2 + tmp4;
 153     tmp8 = tmp2 - tmp4;
 154
 155     tmp1 = p_tmp[4].Re - p_tmp[5].Re;
 156     tmp2 = p_tmp[4].Im - p_tmp[5].Im;
 157     tmp3 = p_tmp[6].Re - p_tmp[7].Re;
 158     tmp4 = p_tmp[6].Im - p_tmp[7].Im;
 159
 160     p_tmp[4].Re = p_tmp[0].Re - tmp5;
 161     p_tmp[0].Re = p_tmp[0].Re + tmp5;
 162     p_tmp[4].Im = p_tmp[0].Im - tmp6;
 163     p_tmp[0].Im = p_tmp[0].Im + tmp6;
 164     p_tmp[6].Re = p_tmp[2].Re + tmp8;
 165     p_tmp[2].Re = p_tmp[2].Re - tmp8;
 166     p_tmp[6].Im = p_tmp[2].Im - tmp7;
 167     p_tmp[2].Im = p_tmp[2].Im + tmp7;
 168
 169     tmp5 = SQRT1_2 * (tmp1 - tmp2);
 170     tmp7 = SQRT1_2 * (tmp3 + tmp4);
 171     tmp6 = SQRT1_2 * (tmp1 + tmp2);
 172     tmp8 = SQRT1_2 * (tmp4 - tmp3);
 173
 174     tmp1 = tmp5 + tmp7;
 175     tmp3 = tmp5 - tmp7;
 176     tmp2 = tmp6 + tmp8;
 177     tmp4 = tmp6 - tmp8;
 178
 179     p_tmp[5].Re = p_tmp[1].Re - tmp1;
 180     p_tmp[1].Re = p_tmp[1].Re + tmp1;
 181     p_tmp[5].Im = p_tmp[1].Im - tmp2;
 182     p_tmp[1].Im = p_tmp[1].Im + tmp2;
 183     p_tmp[7].Re = p_tmp[3].Re + tmp4;
 184     p_tmp[3].Re = p_tmp[3].Re - tmp4;
 185     p_tmp[7].Im = p_tmp[3].Im - tmp3;
 186     p_tmp[3].Im = p_tmp[3].Im + tmp3;
 187   }
 188
 189   step = 1 << (pFFTSpec->order - 4);
 190   n1_4 = 4;
 191   /* Outer loop that loops over FFT stages. */
 192   for (uint32_t fft_stage = 4; fft_stage <= pFFTSpec->order - 2; ++fft_stage) {
 193     n1_2 = 2 * n1_4;
 194     n3_4 = 3 * n1_4;
 195     num_transforms = (num_transforms >> 1) | 1;
 196     for (uint32_t n = 0; n < num_transforms; ++n) {
 197       OMX_U32 offset = pFFTSpec->pOffset[n] << fft_stage;
 198       OMX_FC32* p_tmp = p_buf + offset;
 199
 200       tmp1 = p_tmp[n1_2].Re + p_tmp[n3_4].Re;
 201       tmp2 = p_tmp[n1_2].Re - p_tmp[n3_4].Re;
 202       tmp3 = p_tmp[n1_2].Im + p_tmp[n3_4].Im;
 203       tmp4 = p_tmp[n1_2].Im - p_tmp[n3_4].Im;
 204
 205       p_tmp[n1_2].Re = p_tmp[0].Re - tmp1;
 206       p_tmp[0].Re = p_tmp[0].Re + tmp1;
 207       p_tmp[n1_2].Im = p_tmp[0].Im - tmp3;
 208       p_tmp[0].Im = p_tmp[0].Im + tmp3;
 209       p_tmp[n3_4].Re = p_tmp[n1_4].Re + tmp4;
 210       p_tmp[n1_4].Re = p_tmp[n1_4].Re - tmp4;
 211       p_tmp[n3_4].Im = p_tmp[n1_4].Im - tmp2;
 212       p_tmp[n1_4].Im = p_tmp[n1_4].Im + tmp2;
 213
 214       w_re_ptr = pFFTSpec->pTwiddle + step;
 215       w_im_ptr =
 216           pFFTSpec->pTwiddle + (OMX_U32)(1 << pFFTSpec->order - 2) - step;
 217
 218       /*
 219        * Loop performing split-radix butterfly operations for one sub-transform.
 220        */
 221       for (uint32_t i = 1; i < n1_4; ++i) {
 222         w_re = *w_re_ptr;
 223         w_im = *w_im_ptr;
 224
 225         tmp1 = w_re * p_tmp[n1_2 + i].Re - w_im * p_tmp[n1_2 + i].Im;
 226         tmp2 = w_re * p_tmp[n1_2 + i].Im + w_im * p_tmp[n1_2 + i].Re;
 227         tmp3 = w_re * p_tmp[n3_4 + i].Re + w_im * p_tmp[n3_4 + i].Im;
 228         tmp4 = w_re * p_tmp[n3_4 + i].Im - w_im * p_tmp[n3_4 + i].Re;
 229
 230         tmp5 = tmp1 + tmp3;
 231         tmp1 = tmp1 - tmp3;
 232         tmp6 = tmp2 + tmp4;
 233         tmp2 = tmp2 - tmp4;
 234
 235         p_tmp[n1_2 + i].Re = p_tmp[i].Re - tmp5;
 236         p_tmp[i].Re = p_tmp[i].Re + tmp5;
 237         p_tmp[n1_2 + i].Im = p_tmp[i].Im - tmp6;
 238         p_tmp[i].Im = p_tmp[i].Im + tmp6;
 239         p_tmp[n3_4 + i].Re = p_tmp[n1_4 + i].Re + tmp2;
 240         p_tmp[n1_4 + i].Re = p_tmp[n1_4 + i].Re - tmp2;
 241         p_tmp[n3_4 + i].Im = p_tmp[n1_4 + i].Im - tmp1;
 242         p_tmp[n1_4 + i].Im = p_tmp[n1_4 + i].Im + tmp1;
 243
 244         w_re_ptr += step;
 245         w_im_ptr -= step;
 246       }
 247     }
 248     step >>= 1;
 249     n1_4 <<= 1;
 250   }
 251
 252   /* Last FFT stage, write data to output buffer. */
 253   n1_2 = 2 * n1_4;
 254   n3_4 = 3 * n1_4;
 255   factor = (OMX_F32)2.0f / fft_size;
 256
 257   p_dst = (OMX_FC32*)pDst;
 258
 259   tmp1 = p_buf[n1_2].Re + p_buf[n3_4].Re;
 260   tmp2 = p_buf[n1_2].Re - p_buf[n3_4].Re;
 261   tmp3 = p_buf[n1_2].Im + p_buf[n3_4].Im;
 262   tmp4 = p_buf[n1_2].Im - p_buf[n3_4].Im;
 263
 264   p_dst[n1_2].Re = factor * (p_buf[0].Re - tmp1);
 265   p_dst[0].Re = factor * (p_buf[0].Re + tmp1);
 266   p_dst[n1_2].Im = factor * (p_buf[0].Im - tmp3);
 267   p_dst[0].Im = factor * (p_buf[0].Im + tmp3);
 268   p_dst[n3_4].Re = factor * (p_buf[n1_4].Re + tmp4);
 269   p_dst[n1_4].Re = factor * (p_buf[n1_4].Re - tmp4);
 270   p_dst[n3_4].Im = factor * (p_buf[n1_4].Im - tmp2);
 271   p_dst[n1_4].Im = factor * (p_buf[n1_4].Im + tmp2);
 272
 273   w_re_ptr = pFFTSpec->pTwiddle + step;
 274   w_im_ptr = pFFTSpec->pTwiddle + (OMX_U32)(1 << pFFTSpec->order - 2) - step;
 275
 276   for (uint32_t i = 1; i < n1_4; ++i) {
 277     w_re = *w_re_ptr;
 278     w_im = *w_im_ptr;
 279
 280     tmp1 = w_re * p_buf[n1_2 + i].Re - w_im * p_buf[n1_2 + i].Im;
 281     tmp2 = w_re * p_buf[n1_2 + i].Im + w_im * p_buf[n1_2 + i].Re;
 282     tmp3 = w_re * p_buf[n3_4 + i].Re + w_im * p_buf[n3_4 + i].Im;
 283     tmp4 = w_re * p_buf[n3_4 + i].Im - w_im * p_buf[n3_4 + i].Re;
 284
 285     tmp5 = tmp1 + tmp3;
 286     tmp1 = tmp1 - tmp3;
 287     tmp6 = tmp2 + tmp4;
 288     tmp2 = tmp2 - tmp4;
 289
 290     p_dst[n1_2 + i].Re = factor * (p_buf[i].Re - tmp5);
 291     p_dst[i].Re = factor * (p_buf[i].Re + tmp5);
 292     p_dst[n1_2 + i].Im = factor * (p_buf[i].Im - tmp6);
 293     p_dst[i].Im = factor * (p_buf[i].Im + tmp6);
 294     p_dst[n3_4 + i].Re = factor * (p_buf[n1_4 + i].Re + tmp2);
 295     p_dst[n1_4 + i].Re = factor * (p_buf[n1_4 + i].Re - tmp2);
 296     p_dst[n3_4 + i].Im = factor * (p_buf[n1_4 + i].Im - tmp1);
 297     p_dst[n1_4 + i].Im = factor * (p_buf[n1_4 + i].Im + tmp1);
 298
 299     w_re_ptr += step;
 300     w_im_ptr -= step;
 301   }
 302
 303   return OMX_Sts_NoErr;
 304 }