Mercurial > hg > peckfft

annotate _peck_fft_guts.h @ 0:723f588b82ac

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
import
author Peter Meerwald <p.meerwald@bct-electronic.com>
date Fri, 16 Sep 2011 12:53:08 +0200
parents
children 2d6c49fcafcb
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
1 /*
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
2 Copyright (c) 2003-2010, Mark Borgerding
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
3
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
4 All rights reserved.
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
5
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
6 Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
7
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
8 * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
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.
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
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.
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
11
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
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.
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
13 */
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
14
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
15 /* peck_fft.h
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
16 defines peck_fft_scalar as either short or a float type
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
17 and defines
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
18 typedef struct { peck_fft_scalar r; peck_fft_scalar i; }peck_fft_cpx; */
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
19 #include "peck_fft.h"
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
20 #include <limits.h>
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
21
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
22 #define MAXFACTORS 32
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
23 /* e.g. an fft of length 128 has 4 factors
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
24 as far as kissfft is concerned
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
25 4*4*4*2
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
26 */
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
27
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
28 struct peck_fft_state{
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
29 int nfft;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
30 int inverse;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
31 int factors[2*MAXFACTORS];
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
32 peck_fft_cpx twiddles[1];
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
33 };
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
34
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
35 /*
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
36 Explanation of macros dealing with complex math:
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
37
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
38 C_MUL(m,a,b) : m = a*b
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
39 C_FIXDIV( c , div ) : if a fixed point impl., c /= div. noop otherwise
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
40 C_SUB( res, a,b) : res = a - b
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
41 C_SUBFROM( res , a) : res -= a
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
42 C_ADDTO( res , a) : res += a
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
43 * */
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
44 #ifdef FIXED_POINT
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
45 #if (FIXED_POINT==32)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
46 # define FRACBITS 31
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
47 # define SAMPPROD int64_t
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
48 #define SAMP_MAX 2147483647
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
49 #else
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
50 # define FRACBITS 15
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
51 # define SAMPPROD int32_t
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
52 #define SAMP_MAX 32767
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
53 #endif
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
54
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
55 #define SAMP_MIN -SAMP_MAX
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
56
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
57 #if defined(CHECK_OVERFLOW)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
58 # define CHECK_OVERFLOW_OP(a,op,b) \
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
59 if ( (SAMPPROD)(a) op (SAMPPROD)(b) > SAMP_MAX || (SAMPPROD)(a) op (SAMPPROD)(b) < SAMP_MIN ) { \
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
60 fprintf(stderr,"WARNING:overflow @ " __FILE__ "(%d): (%d " #op" %d) = %ld\n",__LINE__,(a),(b),(SAMPPROD)(a) op (SAMPPROD)(b) ); }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
61 #endif
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
62
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
63
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
64 # define smul(a,b) ( (SAMPPROD)(a)*(b) )
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
65 # define sround( x ) (peck_fft_scalar)( ( (x) + (1<<(FRACBITS-1)) ) >> FRACBITS )
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
66
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
67 # define S_MUL(a,b) sround( smul(a,b) )
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
68
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
69 # define C_MUL(m,a,b) \
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
70 do{ (m).r = sround( smul((a).r,(b).r) - smul((a).i,(b).i) ); \
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
71 (m).i = sround( smul((a).r,(b).i) + smul((a).i,(b).r) ); }while(0)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
72
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
73 # define DIVSCALAR(x,k) \
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
74 (x) = sround( smul( x, SAMP_MAX/k ) )
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
75
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
76 # define C_FIXDIV(c,div) \
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
77 do { DIVSCALAR( (c).r , div); \
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
78 DIVSCALAR( (c).i , div); }while (0)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
79
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
80 # define C_MULBYSCALAR( c, s ) \
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
81 do{ (c).r = sround( smul( (c).r , s ) ) ;\
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
82 (c).i = sround( smul( (c).i , s ) ) ; }while(0)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
83
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
84 #else /* not FIXED_POINT*/
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
85
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
86 # define S_MUL(a,b) ( (a)*(b) )
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
87 #define C_MUL(m,a,b) \
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
88 do{ (m).r = (a).r*(b).r - (a).i*(b).i;\
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
89 (m).i = (a).r*(b).i + (a).i*(b).r; }while(0)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
90 # define C_FIXDIV(c,div) /* NOOP */
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
91 # define C_MULBYSCALAR( c, s ) \
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
92 do{ (c).r *= (s);\
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
93 (c).i *= (s); }while(0)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
94 #endif
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
95
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
96 #ifndef CHECK_OVERFLOW_OP
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
97 # define CHECK_OVERFLOW_OP(a,op,b) /* noop */
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
98 #endif
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
99
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
100 #define C_ADD( res, a,b)\
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
101 do { \
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
102 CHECK_OVERFLOW_OP((a).r,+,(b).r)\
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
103 CHECK_OVERFLOW_OP((a).i,+,(b).i)\
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
104 (res).r=(a).r+(b).r; (res).i=(a).i+(b).i; \
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
105 }while(0)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
106 #define C_SUB( res, a,b)\
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
107 do { \
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
108 CHECK_OVERFLOW_OP((a).r,-,(b).r)\
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
109 CHECK_OVERFLOW_OP((a).i,-,(b).i)\
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
110 (res).r=(a).r-(b).r; (res).i=(a).i-(b).i; \
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
111 }while(0)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
112 #define C_ADDTO( res , a)\
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
113 do { \
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
114 CHECK_OVERFLOW_OP((res).r,+,(a).r)\
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
115 CHECK_OVERFLOW_OP((res).i,+,(a).i)\
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
116 (res).r += (a).r; (res).i += (a).i;\
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
117 }while(0)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
118
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
119 #define C_SUBFROM( res , a)\
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
120 do {\
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
121 CHECK_OVERFLOW_OP((res).r,-,(a).r)\
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
122 CHECK_OVERFLOW_OP((res).i,-,(a).i)\
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
123 (res).r -= (a).r; (res).i -= (a).i; \
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
124 }while(0)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
125
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
126
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
127 #ifdef FIXED_POINT
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
128 # define PECK_FFT_COS(phase) floor(.5+SAMP_MAX * cos (phase))
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
129 # define PECK_FFT_SIN(phase) floor(.5+SAMP_MAX * sin (phase))
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
130 # define HALF_OF(x) ((x)>>1)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
131 #elif defined(USE_SIMD)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
132 # define PECK_FFT_COS(phase) _mm_set1_ps( cos(phase) )
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
133 # define PECK_FFT_SIN(phase) _mm_set1_ps( sin(phase) )
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
134 # define HALF_OF(x) ((x)*_mm_set1_ps(.5))
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
135 #else
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
136 # define PECK_FFT_COS(phase) (peck_fft_scalar) cos(phase)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
137 # define PECK_FFT_SIN(phase) (peck_fft_scalar) sin(phase)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
138 # define HALF_OF(x) ((x)*.5)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
139 #endif
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
140
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
141 #define kf_cexp(x,phase) \
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
142 do{ \
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
143 (x)->r = PECK_FFT_COS(phase);\
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
144 (x)->i = PECK_FFT_SIN(phase);\
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
145 }while(0)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
146
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
147
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
148 /* a debugging function */
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
149 #define pcpx(c)\
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
150 fprintf(stderr,"%g + %gi\n",(double)((c)->r),(double)((c)->i) )
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
151
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
152
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
153 #ifdef PECK_FFT_USE_ALLOCA
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
154 // define this to allow use of alloca instead of malloc for temporary buffers
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
155 // Temporary buffers are used in two case:
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
156 // 1. FFT sizes that have "bad" factors. i.e. not 2,3 and 5
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
157 // 2. "in-place" FFTs. Notice the quotes, since kissfft does not really do an in-place transform.
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
158 #include <alloca.h>
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
159 #define PECK_FFT_TMP_ALLOC(nbytes) alloca(nbytes)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
160 #define PECK_FFT_TMP_FREE(ptr)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
161 #else
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
162 #define PECK_FFT_TMP_ALLOC(nbytes) PECK_FFT_MALLOC(nbytes)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
163 #define PECK_FFT_TMP_FREE(ptr) PECK_FFT_FREE(ptr)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
164 #endif

Repositories maintained by Peter Meerwald, pmeerw@pmeerw.net.