mirror of
https://github.com/mborgerding/kissfft.git
synced 2025-05-27 21:20:27 -04:00
5ebbc5e618
In order to use constants or trigonometric functions with a type other than
double, a suffix ('f' for float or 'l' for long double) has to be used in C.
This commit adds a preprocessor macro 'kiss_fft_suffix' which can be set to
either 'f' or 'l' and which will be added to floating point constants and to
the trigonometric functions (sin and cos).
Without this suffix, the code will use a too high precision for float and a
too low precision for long double.
206 lines
5.9 KiB
C
206 lines
5.9 KiB
C
/*
|
|
* Copyright (c) 2003-2010, Mark Borgerding. All rights reserved.
|
|
* This file is part of KISS FFT - https://github.com/mborgerding/kissfft
|
|
*
|
|
* SPDX-License-Identifier: BSD-3-Clause
|
|
* See COPYING file for more information.
|
|
*/
|
|
|
|
/* 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 "kiss_fft.h"
|
|
#include <limits.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
|
|
#include <stdint.h>
|
|
#if (FIXED_POINT==32)
|
|
# define FRACBITS 31
|
|
# define SAMPPROD int64_t
|
|
#define SAMP_MAX INT32_MAX
|
|
#define SAMP_MIN INT32_MIN
|
|
#else
|
|
# define FRACBITS 15
|
|
# define SAMPPROD int32_t
|
|
#define SAMP_MAX INT16_MAX
|
|
#define SAMP_MIN INT16_MIN
|
|
#endif
|
|
|
|
#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)
|
|
|
|
|
|
// Normally same as kiss_fft_scalar, but for USE_SIMD this is only a single number
|
|
#ifndef kiss_fft_scalar_one
|
|
#define kiss_fft_scalar_one kiss_fft_scalar
|
|
#endif
|
|
|
|
|
|
// If kiss_fft_suffix is not provided, determine it automatically based on kiss_fft_scalar_one
|
|
#if !defined (kiss_fft_suffix) && !defined (kiss_fft_suffix_empty)
|
|
#ifdef FIXED_POINT
|
|
#undef kiss_fft_suffix
|
|
#define kiss_fft_suffix_empty 1
|
|
#endif
|
|
|
|
#define KISS_X_float_X 1
|
|
#define KISS_X__X 2
|
|
#define KISS_X_long_X 3
|
|
#define KISS_concat2(a, b, c) a##b##c
|
|
#define KISS_concat(a, b, c) KISS_concat2(a, b, c)
|
|
|
|
#define double // This is because otherwise KISS_concat(KISS_X_, long double, _X) would produce two tokens which would break the preprocessor comparison
|
|
|
|
#if KISS_concat(KISS_X_, kiss_fft_scalar_one, _X) == KISS_X_float_X // float
|
|
#define kiss_fft_suffix f
|
|
#elif KISS_concat(KISS_X_, kiss_fft_scalar_one, _X) == KISS_X__X // double
|
|
#undef kiss_fft_suffix
|
|
#define kiss_fft_suffix_empty 1
|
|
#elif KISS_concat(KISS_X_, kiss_fft_scalar_one, _X) == KISS_X_long_X // long double
|
|
#define kiss_fft_suffix l
|
|
#else
|
|
#undef kiss_fft_suffix
|
|
#define kiss_fft_suffix_empty 1
|
|
#warning "Unknown kiss_fft_scalar type"
|
|
#endif
|
|
|
|
#undef double
|
|
#endif
|
|
|
|
|
|
#ifdef kiss_fft_suffix_empty
|
|
#define KISS_ADD_SUFFIX(x) x
|
|
#else
|
|
#define KISS_ADD_SUFFIX2(x, y) x##y
|
|
#define KISS_ADD_SUFFIX1(x, y) KISS_ADD_SUFFIX2(x, y)
|
|
#define KISS_ADD_SUFFIX(x) KISS_ADD_SUFFIX1(x, kiss_fft_suffix)
|
|
#endif
|
|
|
|
#ifdef FIXED_POINT
|
|
# define KISS_FFT_COS(phase) floor(.5+SAMP_MAX * cos (phase))
|
|
# define KISS_FFT_SIN(phase) floor(.5+SAMP_MAX * sin (phase))
|
|
# define HALF_OF(x) ((x)>>1)
|
|
#elif defined(USE_SIMD)
|
|
# define KISS_FFT_COS(phase) _mm_set1_ps( cos(phase) )
|
|
# define KISS_FFT_SIN(phase) _mm_set1_ps( sin(phase) )
|
|
# define HALF_OF(x) ((x)*_mm_set1_ps(.5))
|
|
#else
|
|
# define KISS_FFT_COS(phase) (kiss_fft_scalar) KISS_ADD_SUFFIX(cos)(phase)
|
|
# define KISS_FFT_SIN(phase) (kiss_fft_scalar) KISS_ADD_SUFFIX(sin)(phase)
|
|
# define HALF_OF(x) ((x)*KISS_ADD_SUFFIX(.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
|