diff options
Diffstat (limited to 'src/_kiss_fft_guts.h')
| -rw-r--r-- | src/_kiss_fft_guts.h | 194 |
1 files changed, 194 insertions, 0 deletions
diff --git a/src/_kiss_fft_guts.h b/src/_kiss_fft_guts.h new file mode 100644 index 0000000..16305af --- /dev/null +++ b/src/_kiss_fft_guts.h @@ -0,0 +1,194 @@ +/* +Copyright (c) 2003-2010, Mark Borgerding + +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: + + * Redistributions of source code must retain the above copyright notice, +this list of conditions and the following disclaimer. + * 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. + * 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. + +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. +*/ + +/* kiss_fft.h + defines kiss_fft_scalar as either short or a float type + and defines + typedef struct { kiss_fft_scalar r; kiss_fft_scalar i; }kiss_fft_cpx; */ +#include <limits.h> + +#include "kiss_fft.h" + +#define MAXFACTORS 32 +/* e.g. an fft of length 128 has 4 factors + as far as kissfft is concerned + 4*4*4*2 + */ + +struct kiss_fft_state { + int nfft; + int inverse; + int factors[2 * MAXFACTORS]; + kiss_fft_cpx twiddles[1]; +}; + +/* + Explanation of macros dealing with complex math: + + C_MUL(m,a,b) : m = a*b + C_FIXDIV( c , div ) : if a fixed point impl., c /= div. noop otherwise + C_SUB( res, a,b) : res = a - b + C_SUBFROM( res , a) : res -= a + C_ADDTO( res , a) : res += a + * */ +#ifdef FIXED_POINT +#if (FIXED_POINT == 32) +#define FRACBITS 31 +#define SAMPPROD int64_t +#define SAMP_MAX 2147483647 +#else +#define FRACBITS 15 +#define SAMPPROD int32_t +#define SAMP_MAX 32767 +#endif + +#define SAMP_MIN -SAMP_MAX + +#if defined(CHECK_OVERFLOW) +#define CHECK_OVERFLOW_OP(a, op, b) \ + if ((SAMPPROD)(a)op(SAMPPROD)(b) > SAMP_MAX || \ + (SAMPPROD)(a)op(SAMPPROD)(b) < SAMP_MIN) { \ + fprintf(stderr, \ + "WARNING:overflow @ " __FILE__ "(%d): (%d " #op " %d) = %ld\n", \ + __LINE__, (a), (b), (SAMPPROD)(a)op(SAMPPROD)(b)); \ + } +#endif + +#define smul(a, b) ((SAMPPROD)(a) * (b)) +#define sround(x) (kiss_fft_scalar)(((x) + (1 << (FRACBITS - 1))) >> FRACBITS) + +#define S_MUL(a, b) sround(smul(a, b)) + +#define C_MUL(m, a, b) \ + do { \ + (m).r = sround(smul((a).r, (b).r) - smul((a).i, (b).i)); \ + (m).i = sround(smul((a).r, (b).i) + smul((a).i, (b).r)); \ + } while (0) + +#define DIVSCALAR(x, k) (x) = sround(smul(x, SAMP_MAX / k)) + +#define C_FIXDIV(c, div) \ + do { \ + DIVSCALAR((c).r, div); \ + DIVSCALAR((c).i, div); \ + } while (0) + +#define C_MULBYSCALAR(c, s) \ + do { \ + (c).r = sround(smul((c).r, s)); \ + (c).i = sround(smul((c).i, s)); \ + } while (0) + +#else /* not FIXED_POINT*/ + +#define S_MUL(a, b) ((a) * (b)) +#define C_MUL(m, a, b) \ + do { \ + (m).r = (a).r * (b).r - (a).i * (b).i; \ + (m).i = (a).r * (b).i + (a).i * (b).r; \ + } while (0) +#define C_FIXDIV(c, div) /* NOOP */ +#define C_MULBYSCALAR(c, s) \ + do { \ + (c).r *= (s); \ + (c).i *= (s); \ + } while (0) +#endif + +#ifndef CHECK_OVERFLOW_OP +#define CHECK_OVERFLOW_OP(a, op, b) /* noop */ +#endif + +#define C_ADD(res, a, b) \ + do { \ + CHECK_OVERFLOW_OP((a).r, +, (b).r) \ + CHECK_OVERFLOW_OP((a).i, +, (b).i) \ + (res).r = (a).r + (b).r; \ + (res).i = (a).i + (b).i; \ + } while (0) +#define C_SUB(res, a, b) \ + do { \ + CHECK_OVERFLOW_OP((a).r, -, (b).r) \ + CHECK_OVERFLOW_OP((a).i, -, (b).i) \ + (res).r = (a).r - (b).r; \ + (res).i = (a).i - (b).i; \ + } while (0) +#define C_ADDTO(res, a) \ + do { \ + CHECK_OVERFLOW_OP((res).r, +, (a).r) \ + CHECK_OVERFLOW_OP((res).i, +, (a).i) \ + (res).r += (a).r; \ + (res).i += (a).i; \ + } while (0) + +#define C_SUBFROM(res, a) \ + do { \ + CHECK_OVERFLOW_OP((res).r, -, (a).r) \ + CHECK_OVERFLOW_OP((res).i, -, (a).i) \ + (res).r -= (a).r; \ + (res).i -= (a).i; \ + } while (0) + +#ifdef FIXED_POINT +#define KISS_FFT_COS(phase) floorf(.5 + SAMP_MAX * cosf(phase)) +#define KISS_FFT_SIN(phase) floorf(.5 + SAMP_MAX * sinf(phase)) +#define HALF_OF(x) ((x) >> 1) +#elif defined(USE_SIMD) +#define KISS_FFT_COS(phase) _mm_set1_ps(cosf(phase)) +#define KISS_FFT_SIN(phase) _mm_set1_ps(sinf(phase)) +#define HALF_OF(x) ((x)*_mm_set1_ps(.5)) +#else +#define KISS_FFT_COS(phase) (kiss_fft_scalar) cosf(phase) +#define KISS_FFT_SIN(phase) (kiss_fft_scalar) sinf(phase) +#define HALF_OF(x) ((x)*.5) +#endif + +#define kf_cexp(x, phase) \ + do { \ + (x)->r = KISS_FFT_COS(phase); \ + (x)->i = KISS_FFT_SIN(phase); \ + } while (0) + +/* a debugging function */ +#define pcpx(c) \ + fprintf(stderr, "%g + %gi\n", (double)((c)->r), (double)((c)->i)) + +#ifdef KISS_FFT_USE_ALLOCA +// define this to allow use of alloca instead of malloc for temporary buffers +// Temporary buffers are used in two case: +// 1. FFT sizes that have "bad" factors. i.e. not 2,3 and 5 +// 2. "in-place" FFTs. Notice the quotes, since kissfft does not really do an +// in-place transform. +#include <alloca.h> +#define KISS_FFT_TMP_ALLOC(nbytes) alloca(nbytes) +#define KISS_FFT_TMP_FREE(ptr) +#else +#define KISS_FFT_TMP_ALLOC(nbytes) KISS_FFT_MALLOC(nbytes) +#define KISS_FFT_TMP_FREE(ptr) KISS_FFT_FREE(ptr) +#endif |