1 /*

     2 Copyright (c) 2003-2004, Mark Borgerding

     3 

     4 All rights reserved.

     5 

     6 Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

     7 

     8     * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.

     9     * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.

    10     * Neither the author nor the names of any contributors may be used to endorse or promote products derived from this software without specific prior written permission.

    11 

    12 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

    13 */

    14 

    15 #include "kiss_fftr_f64.h"

    16 #include "_kiss_fft_guts_f64.h"

    17 

    18 struct kiss_fftr_f64_state

    19 {

    20   kiss_fft_f64_cfg substate;

    21   kiss_fft_f64_cpx *tmpbuf;

    22   kiss_fft_f64_cpx *super_twiddles;

    23 #ifdef USE_SIMD

    24   long pad;

    25 #endif

    26 };

    27 #ifdef __SYMBIAN32__

    28 EXPORT_C

    29 #endif

    30 

    31 

    32 kiss_fftr_f64_cfg

    33 kiss_fftr_f64_alloc (int nfft, int inverse_fft, void *mem, size_t * lenmem)

    34 {

    35   int i;

    36   kiss_fftr_f64_cfg st = NULL;

    37   size_t subsize, memneeded;

    38 

    39   if (nfft & 1) {

    40     fprintf (stderr, "Real FFT optimization must be even.\n");

    41     return NULL;

    42   }

    43   nfft >>= 1;

    44 

    45   kiss_fft_f64_alloc (nfft, inverse_fft, NULL, &subsize);

    46   memneeded =

    47       sizeof (struct kiss_fftr_f64_state) + subsize +

    48       sizeof (kiss_fft_f64_cpx) * (nfft * 2);

    49 

    50   if (lenmem == NULL) {

    51     st = (kiss_fftr_f64_cfg) KISS_FFT_F64_MALLOC (memneeded);

    52   } else {

    53     if (*lenmem >= memneeded)

    54       st = (kiss_fftr_f64_cfg) mem;

    55     *lenmem = memneeded;

    56   }

    57   if (!st)

    58     return NULL;

    59 

    60   st->substate = (kiss_fft_f64_cfg) (st + 1);   /*just beyond kiss_fftr_f64_state struct */

    61   st->tmpbuf = (kiss_fft_f64_cpx *) (((char *) st->substate) + subsize);

    62   st->super_twiddles = st->tmpbuf + nfft;

    63   kiss_fft_f64_alloc (nfft, inverse_fft, st->substate, &subsize);

    64 

    65   for (i = 0; i < nfft; ++i) {

    66     double phase = -3.14159265358979323846264338327 * ((double) i / nfft + .5);

    67 

    68     if (inverse_fft)

    69       phase *= -1;

    70     kf_cexp (st->super_twiddles + i, phase);

    71   }

    72   return st;

    73 }

    74 #ifdef __SYMBIAN32__

    75 EXPORT_C

    76 #endif

    77 

    78 

    79 void

    80 kiss_fftr_f64 (kiss_fftr_f64_cfg st, const kiss_fft_f64_scalar * timedata,

    81     kiss_fft_f64_cpx * freqdata)

    82 {

    83   /* input buffer timedata is stored row-wise */

    84   int k, ncfft;

    85   kiss_fft_f64_cpx fpnk, fpk, f1k, f2k, tw, tdc;

    86 

    87   if (st->substate->inverse) {

    88     fprintf (stderr, "kiss fft usage error: improper alloc\n");

    89     exit (1);

    90   }

    91 

    92   ncfft = st->substate->nfft;

    93 

    94   /*perform the parallel fft of two real signals packed in real,imag */

    95   kiss_fft_f64 (st->substate, (const kiss_fft_f64_cpx *) timedata, st->tmpbuf);

    96   /* The real part of the DC element of the frequency spectrum in st->tmpbuf

    97    * contains the sum of the even-numbered elements of the input time sequence

    98    * The imag part is the sum of the odd-numbered elements

    99    *

   100    * The sum of tdc.r and tdc.i is the sum of the input time sequence. 

   101    *      yielding DC of input time sequence

   102    * The difference of tdc.r - tdc.i is the sum of the input (dot product) [1,-1,1,-1... 

   103    *      yielding Nyquist bin of input time sequence

   104    */

   105 

   106   tdc.r = st->tmpbuf[0].r;

   107   tdc.i = st->tmpbuf[0].i;

   108   C_FIXDIV (tdc, 2);

   109   CHECK_OVERFLOW_OP (tdc.r, +, tdc.i);

   110   CHECK_OVERFLOW_OP (tdc.r, -, tdc.i);

   111   freqdata[0].r = tdc.r + tdc.i;

   112   freqdata[ncfft].r = tdc.r - tdc.i;

   113 #ifdef USE_SIMD

   114   freqdata[ncfft].i = freqdata[0].i = _mm_set1_ps (0);

   115 #else

   116   freqdata[ncfft].i = freqdata[0].i = 0;

   117 #endif

   118 

   119   for (k = 1; k <= ncfft / 2; ++k) {

   120     fpk = st->tmpbuf[k];

   121     fpnk.r = st->tmpbuf[ncfft - k].r;

   122     fpnk.i = -st->tmpbuf[ncfft - k].i;

   123     C_FIXDIV (fpk, 2);

   124     C_FIXDIV (fpnk, 2);

   125 

   126     C_ADD (f1k, fpk, fpnk);

   127     C_SUB (f2k, fpk, fpnk);

   128     C_MUL (tw, f2k, st->super_twiddles[k]);

   129 

   130     freqdata[k].r = HALF_OF (f1k.r + tw.r);

   131     freqdata[k].i = HALF_OF (f1k.i + tw.i);

   132     freqdata[ncfft - k].r = HALF_OF (f1k.r - tw.r);

   133     freqdata[ncfft - k].i = HALF_OF (tw.i - f1k.i);

   134   }

   135 }

   136 #ifdef __SYMBIAN32__

   137 EXPORT_C

   138 #endif

   139 

   140 

   141 void

   142 kiss_fftri_f64 (kiss_fftr_f64_cfg st, const kiss_fft_f64_cpx * freqdata,

   143     kiss_fft_f64_scalar * timedata)

   144 {

   145   /* input buffer timedata is stored row-wise */

   146   int k, ncfft;

   147 

   148   if (st->substate->inverse == 0) {

   149     fprintf (stderr, "kiss fft usage error: improper alloc\n");

   150     exit (1);

   151   }

   152 

   153   ncfft = st->substate->nfft;

   154 

   155   st->tmpbuf[0].r = freqdata[0].r + freqdata[ncfft].r;

   156   st->tmpbuf[0].i = freqdata[0].r - freqdata[ncfft].r;

   157   C_FIXDIV (st->tmpbuf[0], 2);

   158 

   159   for (k = 1; k <= ncfft / 2; ++k) {

   160     kiss_fft_f64_cpx fk, fnkc, fek, fok, tmp;

   161 

   162     fk = freqdata[k];

   163     fnkc.r = freqdata[ncfft - k].r;

   164     fnkc.i = -freqdata[ncfft - k].i;

   165     C_FIXDIV (fk, 2);

   166     C_FIXDIV (fnkc, 2);

   167 

   168     C_ADD (fek, fk, fnkc);

   169     C_SUB (tmp, fk, fnkc);

   170     C_MUL (fok, tmp, st->super_twiddles[k]);

   171     C_ADD (st->tmpbuf[k], fek, fok);

   172     C_SUB (st->tmpbuf[ncfft - k], fek, fok);

   173 #ifdef USE_SIMD

   174     st->tmpbuf[ncfft - k].i *= _mm_set1_ps (-1.0);

   175 #else

   176     st->tmpbuf[ncfft - k].i *= -1;

   177 #endif

   178   }

   179   kiss_fft_f64 (st->substate, st->tmpbuf, (kiss_fft_f64_cpx *) timedata);

   180 }