Line data    Source code

       1             : /******************************************************************************
       2             : *                        ETSI TS 103 634 V1.5.1                               *
       3             : *              Low Complexity Communication Codec Plus (LC3plus)              *
       4             : *                                                                             *
       5             : * Copyright licence is solely granted through ETSI Intellectual Property      *
       6             : * Rights Policy, 3rd April 2019. No patent licence is granted by implication, *
       7             : * estoppel or otherwise.                                                      *
       8             : ******************************************************************************/
       9             : 
      10             : #include "options.h"
      11             : #include "wmc_auto.h"
      12             : #include <assert.h>
      13             : #include <string.h> /* for mmove */
      14             : #include <stdio.h>
      15             : #include <stdlib.h>
      16             : #include "iisfft.h"
      17             : #include "cfft.h"
      18             : 
      19             : /*  the fixed length fft functions have been split into sevelral headers to
      20             :     have smaller files. to give the compiler more room to optimize the ffts
      21             :     can't be in separate compilation units. the header approach seemed to be
      22             :     the best compromise. to prevent them being included from anywhere else,
      23             :     they are guarded by the INCLUDED_FROM_IISFFT_C macro.
      24             : */
      25             : #define INCLUDED_FROM_IISFFT_C
      26             : #include "fft_2_9.h"
      27             : #include "fft_15_16.h"
      28             : #include "fft_32.h"
      29             : #include "fft_60_128.h"
      30             : #include "fft_240_480.h"
      31             : #include "fft_384_768.h"
      32             : #include "fft_generic.h"
      33             : 
      34             : 
      35           0 : void LC3_iisfft_apply(Iisfft* handle, LC3_FLOAT* x)
      36             : {
      37           0 :     if (handle->sign == -1) {
      38           0 :         if (!fft_n(x, handle->length))
      39             :         {
      40             :             LC3_FLOAT scratch[2*MAX_LEN];
      41           0 :             pfaDFT(x, handle->length, scratch, handle->num_factors, handle->factors, handle->scratch2,
      42           0 :                    handle->isPrime);
      43             :         }
      44             :     } else {
      45           0 :         assert(0);
      46             :     }
      47           0 : }
      48             : 
      49             : /* returns 1 if there is no specialized function for length or 1 if a scratch needs to be allocated.
      50             :    check the fft_n function */
      51           0 : static LC3_INT need_scratch(LC3_INT n)
      52             : {
      53           0 :     return n != 2 && n != 3 && n != 4 && n != 5 && n != 7 && n != 8 && n != 9 && n != 15 && n != 16 && n != 32 &&
      54           0 :            n != 60 && n != 64 && n != 128 && n != 192 && n != 240 && n != 256 && n != 384 && n != 480 && n != 512 && n != 768 &&
      55             :            n != 1024;
      56             : }
      57             : 
      58           0 : IIS_FFT_ERROR LC3_iisfft_plan(Iisfft* handle, LC3_INT length, LC3_INT sign)
      59             : {
      60           0 :     memset(handle, 0, sizeof(Iisfft));
      61           0 :     if (length < 2)
      62             :     {
      63           0 :         return IIS_FFT_LENGTH_ERROR;
      64             :     }
      65           0 :     handle->length = length;
      66           0 :     handle->sign   = sign;
      67           0 :     if (need_scratch(length)) {
      68             :         /* only needed for prime numbers bigger than BORDER_FOR_SECOND_SCRATCH */
      69           0 :         LC3_INT i                    = 0;
      70           0 :         LC3_INT lengthOfPrimeScratch = BORDER_FOR_SECOND_SCRATCH;
      71           0 :         if (!factorize(length, &handle->num_factors, handle->factors, handle->isPrime))
      72             :         {
      73           0 :             return IIS_FFT_LENGTH_ERROR;
      74             :         }
      75             :         /* create additional scratch for primeFFT() */
      76           0 :         for (i = 0; i < handle->num_factors; i++) {
      77           0 :             if (handle->isPrime[i] == 1 && handle->factors[i] > lengthOfPrimeScratch) {
      78           0 :                 lengthOfPrimeScratch = handle->factors[i];
      79             :             }
      80             :         }
      81           0 :         if (lengthOfPrimeScratch > BORDER_FOR_SECOND_SCRATCH) {
      82           0 :             handle->scratch2 = (LC3_INT*)malloc(sizeof(LC3_INT) * lengthOfPrimeScratch);
      83           0 :             if (!handle->scratch2)
      84             :             {
      85           0 :                 return IIS_FFT_MEMORY_ERROR;
      86             :             }
      87             :         }
      88             :     }
      89             : 
      90           0 :     return IIS_FFT_NO_ERROR;
      91             : }
      92             : 
      93           0 : void LC3_iisfft_free(Iisfft* handle)
      94             : {
      95           0 :     handle->length = 0;
      96           0 :     if (handle->scratch2)
      97             :     {
      98           0 :         free(handle->scratch2);
      99             :     }
     100           0 : }
     101             : 
     102             : 
     103             : 
     104             : /*  generate sine table needed by LC3_rfft_pre/rfft/post. the table must be freed with iisFree */
     105           0 : void LC3_create_sine_table(LC3_INT32 len, LC3_FLOAT *sine_table)
     106             : {
     107             :        LC3_INT32 i;
     108             : 
     109           0 :     for (i = 0; i < len / 2 + 1; i++) {
     110           0 :         sine_table[i] = (LC3_FLOAT)sin(2.0 * M_PIl * i / len);
     111             :     }
     112           0 : }
     113             : 
     114           0 : void LC3_rfft_post(const LC3_FLOAT* restrict sine_table, LC3_FLOAT* restrict buf, LC3_INT32 len)
     115             : {
     116             :     LC3_FLOAT tmp1, tmp2, tmp3, tmp4, s, c;
     117             :     LC3_INT32 i;
     118             : 
     119           0 :     tmp1 = buf[0] + buf[1];
     120           0 :     buf[1] = buf[0] - buf[1];
     121           0 :     buf[0] = tmp1;
     122             : 
     123           0 :     for (i = 1; i <= (len + 2) / 4; i++) {
     124           0 :         s = sine_table[i];           /* sin(pi*i/(len/2)) */
     125           0 :         c = sine_table[i + len / 4]; /* cos(pi*i/(len/2)) */
     126             : 
     127           0 :         tmp1 = buf[2 * i] - buf[len - 2 * i];
     128           0 :         tmp2 = buf[2 * i + 1] + buf[len - 2 * i + 1];
     129           0 :         tmp3 = s * tmp1 - c * tmp2; /* real part of j*W(k,N)*[T(k) - T'(N-k)] */
     130           0 :         tmp4 = c * tmp1 + s * tmp2; /* imag part of j*W(k,N)*[T(k) - T'(N-k)] */
     131           0 :         tmp1 = buf[2 * i] + buf[len - 2 * i];
     132           0 :         tmp2 = buf[2 * i + 1] - buf[len - 2 * i + 1];
     133             : 
     134           0 :         buf[2 * i] = 0.5f * (tmp1 - tmp3);
     135           0 :         buf[2 * i + 1] = 0.5f * (tmp2 - tmp4);
     136           0 :         buf[len - 2 * i] = 0.5f * (tmp1 + tmp3);
     137           0 :         buf[len - 2 * i + 1] = -0.5f * (tmp2 + tmp4);
     138             :     }
     139             : 
     140           0 :     return;
     141             : }
     142             : 
     143           0 : void LC3_rfft_pre(const LC3_FLOAT* restrict sine_table, LC3_FLOAT* restrict buf, LC3_INT32 len)
     144             : {
     145             :     LC3_FLOAT scale ;
     146             :     LC3_FLOAT tmp1, tmp2, tmp3, tmp4, s, c;
     147             :     LC3_INT32 i;
     148             :    
     149           0 :     scale = 1.0f / len; /* constant */
     150             : 
     151           0 :     tmp1 = buf[0] + buf[1];
     152           0 :     buf[1] = scale * (buf[0] - buf[1]);
     153           0 :     buf[0] = scale * tmp1;
     154             : 
     155           0 :     for (i = 1; i <= (len + 2) / 4; i++) {
     156           0 :         s = sine_table[i];           /* sin(pi*i/(len/2)) */
     157           0 :         c = sine_table[i + len / 4]; /* cos(pi*i/(len/2)) */
     158             : 
     159           0 :         tmp1 = buf[2 * i] - buf[len - 2 * i];
     160           0 :         tmp2 = buf[2 * i + 1] + buf[len - 2 * i + 1];
     161           0 :         tmp3 = s * tmp1 + c * tmp2;  /* real part of j*W(k,N)*[T(k) - T'(N-k)] */
     162           0 :         tmp4 = -c * tmp1 + s * tmp2; /* imag part of j*W(k,N)*[T(k) - T'(N-k)] */
     163           0 :         tmp1 = buf[2 * i] + buf[len - 2 * i];
     164           0 :         tmp2 = buf[2 * i + 1] - buf[len - 2 * i + 1];
     165             : 
     166           0 :         buf[2 * i] = scale * (tmp1 + tmp3);
     167           0 :         buf[2 * i + 1] = -scale * (tmp2 + tmp4);
     168           0 :         buf[len - 2 * i] = scale * (tmp1 - tmp3);
     169           0 :         buf[len - 2 * i + 1] = scale * (tmp2 - tmp4);
     170             :     }
     171           0 :      return;
     172             : }
     173             :