Mercurial > hg > peckfft

annotate peck_fft.c @ 2:3d08140650d8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
cleanup
author Peter Meerwald <p.meerwald@bct-electronic.com>
date Fri, 16 Sep 2011 13:02:34 +0200
parents cfec79393811
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
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
16 #include "_peck_fft_guts.h"
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
17 /* The guts header contains all the multiplication and addition macros that are defined for
1
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
18 * fixed or floating point complex numbers. It also delares the kf_ internal functions.
0
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
19 */
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
20
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
21 static void kf_bfly2(
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
22 peck_fft_cpx * Fout,
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
23 const size_t fstride,
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
24 const peck_fft_cfg st,
1
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
25 int m) {
0
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
26
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
27 //printf("kf_bfly2\n");
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
28
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
29 peck_fft_cpx * Fout2;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
30 peck_fft_cpx * tw1 = st->twiddles;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
31 peck_fft_cpx t;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
32 Fout2 = Fout + m;
1
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
33 do {
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
34 C_FIXDIV(*Fout, 2);
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
35 C_FIXDIV(*Fout2, 2);
0
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
36
1
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
37 C_MUL(t, *Fout2, *tw1);
0
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
38 tw1 += fstride;
1
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
39 C_SUB(*Fout2, *Fout, t);
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
40 C_ADDTO(*Fout, t);
0
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
41 ++Fout2;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
42 ++Fout;
1
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
43 } while (--m);
0
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
44 }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
45
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
46 static void kf_bfly4(
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
47 peck_fft_cpx * Fout,
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
48 const size_t fstride,
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
49 const peck_fft_cfg st,
1
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
50 const size_t m) {
0
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
51 peck_fft_cpx *tw1,*tw2,*tw3;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
52 peck_fft_cpx scratch[6];
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
53 size_t k=m;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
54 const size_t m2=2*m;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
55 const size_t m3=3*m;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
56
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
57 //printf("kf_bfly4\n");
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
58
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
59 tw3 = tw2 = tw1 = st->twiddles;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
60
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
61 do {
1
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
62 C_FIXDIV(*Fout, 4);
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
63 C_FIXDIV(Fout[m], 4);
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
64 C_FIXDIV(Fout[m2], 4);
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
65 C_FIXDIV(Fout[m3], 4);
0
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
66
1
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
67 C_MUL(scratch[0], Fout[m], *tw1);
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
68 C_MUL(scratch[1], Fout[m2], *tw2);
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
69 C_MUL(scratch[2], Fout[m3], *tw3);
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
70
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
71 C_SUB(scratch[5], *Fout, scratch[1]);
0
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
72 C_ADDTO(*Fout, scratch[1]);
1
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
73 C_ADD(scratch[3], scratch[0], scratch[2]);
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
74 C_SUB(scratch[4], scratch[0], scratch[2]);
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
75 C_SUB(Fout[m2], *Fout, scratch[3]);
0
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
76 tw1 += fstride;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
77 tw2 += fstride*2;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
78 tw3 += fstride*3;
1
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
79 C_ADDTO(*Fout, scratch[3]);
0
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
80
1
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
81 if (st->inverse) {
0
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
82 Fout[m].r = scratch[5].r - scratch[4].i;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
83 Fout[m].i = scratch[5].i + scratch[4].r;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
84 Fout[m3].r = scratch[5].r + scratch[4].i;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
85 Fout[m3].i = scratch[5].i - scratch[4].r;
1
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
86 } else {
0
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
87 Fout[m].r = scratch[5].r + scratch[4].i;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
88 Fout[m].i = scratch[5].i - scratch[4].r;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
89 Fout[m3].r = scratch[5].r - scratch[4].i;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
90 Fout[m3].i = scratch[5].i + scratch[4].r;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
91 }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
92 ++Fout;
1
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
93 } while (--k);
0
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
94 }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
95
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
96 static void kf_bfly3(
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
97 peck_fft_cpx * Fout,
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
98 const size_t fstride,
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
99 const peck_fft_cfg st,
1
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
100 size_t m) {
0
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
101 size_t k=m;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
102 const size_t m2 = 2*m;
1
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
103 peck_fft_cpx *tw1, *tw2;
0
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
104 peck_fft_cpx scratch[5];
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
105 peck_fft_cpx epi3;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
106 epi3 = st->twiddles[fstride*m];
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
107
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
108 printf("kf_bfly3\n");
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
109
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
110
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
111 tw1=tw2=st->twiddles;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
112
1
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
113 do {
0
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
114 C_FIXDIV(*Fout,3); C_FIXDIV(Fout[m],3); C_FIXDIV(Fout[m2],3);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
115
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
116 C_MUL(scratch[1],Fout[m] , *tw1);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
117 C_MUL(scratch[2],Fout[m2] , *tw2);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
118
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
119 C_ADD(scratch[3],scratch[1],scratch[2]);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
120 C_SUB(scratch[0],scratch[1],scratch[2]);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
121 tw1 += fstride;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
122 tw2 += fstride*2;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
123
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
124 Fout[m].r = Fout->r - HALF_OF(scratch[3].r);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
125 Fout[m].i = Fout->i - HALF_OF(scratch[3].i);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
126
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
127 C_MULBYSCALAR( scratch[0] , epi3.i );
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
128
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
129 C_ADDTO(*Fout,scratch[3]);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
130
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
131 Fout[m2].r = Fout[m].r + scratch[0].i;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
132 Fout[m2].i = Fout[m].i - scratch[0].r;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
133
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
134 Fout[m].r -= scratch[0].i;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
135 Fout[m].i += scratch[0].r;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
136
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
137 ++Fout;
1
cfec79393811 cleanup
Peter Meerwald <p.meerwald@bct-electronic.com>
parents: 0
diff changeset
138 } while (--k);
0
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
139 }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
140
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
141 static void kf_bfly5(
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
142 peck_fft_cpx * Fout,
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
143 const size_t fstride,
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
144 const peck_fft_cfg st,
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
145 int m
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 peck_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
149 int u;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
150 peck_fft_cpx scratch[13];
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
151 peck_fft_cpx * twiddles = st->twiddles;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
152 peck_fft_cpx *tw;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
153 peck_fft_cpx ya,yb;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
154 ya = twiddles[fstride*m];
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
155 yb = twiddles[fstride*2*m];
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
156
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
157 printf("kf_bfly5\n");
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
158
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
159
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
160 Fout0=Fout;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
161 Fout1=Fout0+m;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
162 Fout2=Fout0+2*m;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
163 Fout3=Fout0+3*m;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
164 Fout4=Fout0+4*m;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
165
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
166 tw=st->twiddles;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
167 for ( u=0; u<m; ++u ) {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
168 C_FIXDIV( *Fout0,5); C_FIXDIV( *Fout1,5); C_FIXDIV( *Fout2,5); C_FIXDIV( *Fout3,5); C_FIXDIV( *Fout4,5);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
169 scratch[0] = *Fout0;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
170
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
171 C_MUL(scratch[1] ,*Fout1, tw[u*fstride]);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
172 C_MUL(scratch[2] ,*Fout2, tw[2*u*fstride]);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
173 C_MUL(scratch[3] ,*Fout3, tw[3*u*fstride]);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
174 C_MUL(scratch[4] ,*Fout4, tw[4*u*fstride]);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
175
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
176 C_ADD( scratch[7],scratch[1],scratch[4]);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
177 C_SUB( scratch[10],scratch[1],scratch[4]);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
178 C_ADD( scratch[8],scratch[2],scratch[3]);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
179 C_SUB( scratch[9],scratch[2],scratch[3]);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
180
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
181 Fout0->r += scratch[7].r + scratch[8].r;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
182 Fout0->i += scratch[7].i + scratch[8].i;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
183
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
184 scratch[5].r = scratch[0].r + S_MUL(scratch[7].r,ya.r) + S_MUL(scratch[8].r,yb.r);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
185 scratch[5].i = scratch[0].i + S_MUL(scratch[7].i,ya.r) + S_MUL(scratch[8].i,yb.r);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
186
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
187 scratch[6].r = S_MUL(scratch[10].i,ya.i) + S_MUL(scratch[9].i,yb.i);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
188 scratch[6].i = -S_MUL(scratch[10].r,ya.i) - S_MUL(scratch[9].r,yb.i);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
189
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
190 C_SUB(*Fout1,scratch[5],scratch[6]);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
191 C_ADD(*Fout4,scratch[5],scratch[6]);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
192
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
193 scratch[11].r = scratch[0].r + S_MUL(scratch[7].r,yb.r) + S_MUL(scratch[8].r,ya.r);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
194 scratch[11].i = scratch[0].i + S_MUL(scratch[7].i,yb.r) + S_MUL(scratch[8].i,ya.r);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
195 scratch[12].r = - S_MUL(scratch[10].i,yb.i) + S_MUL(scratch[9].i,ya.i);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
196 scratch[12].i = S_MUL(scratch[10].r,yb.i) - S_MUL(scratch[9].r,ya.i);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
197
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
198 C_ADD(*Fout2,scratch[11],scratch[12]);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
199 C_SUB(*Fout3,scratch[11],scratch[12]);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
200
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
201 ++Fout0;++Fout1;++Fout2;++Fout3;++Fout4;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
202 }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
203 }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
204
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
205 /* perform the butterfly for one stage of a mixed radix FFT */
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
206 static void kf_bfly_generic(
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
207 peck_fft_cpx * Fout,
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
208 const size_t fstride,
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
209 const peck_fft_cfg st,
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
210 int m,
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
211 int p
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
212 )
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
213 {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
214 int u,k,q1,q;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
215 peck_fft_cpx * twiddles = st->twiddles;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
216 peck_fft_cpx t;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
217 int Norig = st->nfft;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
218
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
219 printf("kf_bfly_generic\n");
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
220
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
221
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
222 peck_fft_cpx * scratch = (peck_fft_cpx*)PECK_FFT_TMP_ALLOC(sizeof(peck_fft_cpx)*p);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
223
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
224 for ( u=0; u<m; ++u ) {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
225 k=u;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
226 for (q1 = 0; q1 < p; ++q1) {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
227 scratch[q1] = Fout[k];
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
228 C_FIXDIV(scratch[q1], p);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
229 k += m;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
230 }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
231
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
232 k=u;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
233 for ( q1=0 ; q1<p ; ++q1 ) {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
234 int twidx=0;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
235 Fout[ k ] = scratch[0];
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
236 for (q=1;q<p;++q ) {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
237 twidx += fstride * k;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
238 if (twidx>=Norig) twidx-=Norig;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
239 C_MUL(t,scratch[q] , twiddles[twidx] );
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
240 C_ADDTO( Fout[ k ] ,t);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
241 }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
242 k += m;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
243 }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
244 }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
245 PECK_FFT_TMP_FREE(scratch);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
246 }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
247
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
248 static void kf_work(
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
249 peck_fft_cpx * Fout,
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
250 const peck_fft_cpx * f,
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
251 const size_t fstride,
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
252 int *factors,
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
253 const peck_fft_cfg st) {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
254 peck_fft_cpx *Fout_beg = Fout;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
255 const int p = *factors++; /* the radix */
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
256 const int m = *factors++; /* stage's FFT length / p */
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
257 const peck_fft_cpx *Fout_end = Fout + p*m;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
258
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
259 // printf("kf_work\n");
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
260
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
261 if (m == 1) {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
262 do {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
263 *Fout = *f;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
264 f += fstride;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
265 } while (++Fout != Fout_end);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
266 } else {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
267 do {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
268 // recursive call:
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
269 // DFT of size m*p performed by doing
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
270 // p instances of smaller DFTs of size m,
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
271 // each one takes a decimated version of the input
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
272 kf_work(Fout, f, fstride*p, factors, st);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
273 f += fstride;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
274 } while ((Fout += m) != Fout_end);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
275 }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
276
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
277 Fout=Fout_beg;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
278
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
279 // recombine the p smaller DFTs
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
280 switch (p) {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
281 case 2: kf_bfly2(Fout, fstride, st, m); break;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
282 case 3: kf_bfly3(Fout, fstride, st, m); break;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
283 case 4: kf_bfly4(Fout, fstride, st, m); break;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
284 case 5: kf_bfly5(Fout, fstride, st, m); break;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
285 default: kf_bfly_generic(Fout, fstride, st, m, p); break;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
286 }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
287 }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
288
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
289 /*
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
290 * facbuf is populated by p1, m1, p2, m2, ...
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
291 * where
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
292 * p[i] * m[i] = m[i-1]
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
293 * m0 = n
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
294 */
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
295 static void kf_factor(int n, int * facbuf) {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
296 int p = 4;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
297 float floor_sqrt;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
298 floor_sqrt = floorf(sqrtf(n));
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
299
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
300 /* factor out powers of 4, powers of 2, then any remaining primes */
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
301 do {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
302 while (n % p) {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
303 switch (p) {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
304 case 4: p = 2; break;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
305 case 2: p = 3; break;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
306 default: p += 2; break;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
307 }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
308 if (p > floor_sqrt)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
309 p = n; /* no more factors, skip to end */
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
310 }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
311 n /= p;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
312 *facbuf++ = p;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
313 *facbuf++ = n;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
314 } while (n > 1);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
315 }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
316
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
317 /*
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
318 * User-callable function to allocate all necessary storage space for the fft.
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
319 * The return value is a contiguous block of memory, allocated with malloc. As such,
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
320 * it can be freed with free(), rather than a peck_fft-specific function.
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
321 */
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
322 peck_fft_cfg peck_fft_alloc(int nfft, int inverse_fft, void * mem, size_t * lenmem) {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
323 peck_fft_cfg st = NULL;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
324 size_t memneeded = sizeof(struct peck_fft_state)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
325 + sizeof(peck_fft_cpx)*(nfft-1); /* twiddle factors */
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
326
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
327 if (lenmem == NULL) {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
328 st = ( peck_fft_cfg)PECK_FFT_MALLOC(memneeded);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
329 } else {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
330 if (mem != NULL && *lenmem >= memneeded)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
331 st = (peck_fft_cfg)mem;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
332 *lenmem = memneeded;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
333 }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
334
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
335 if (st) {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
336 int i;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
337 st->nfft=nfft;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
338 st->inverse = inverse_fft;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
339
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
340 for (i = 0; i < nfft; ++i) {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
341 const float pi = 3.14159265359f;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
342 float phase = -2*pi*i / nfft;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
343 if (st->inverse)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
344 phase *= -1;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
345 kf_cexp(st->twiddles+i, phase);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
346 }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
347
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
348 kf_factor(nfft, st->factors);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
349 }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
350 return st;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
351 }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
352
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
353 void peck_fft(peck_fft_cfg cfg, const peck_fft_cpx *fin, peck_fft_cpx *fout) {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
354 kf_work(fout, fin, 1, cfg->factors, cfg);
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
355 }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
356
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
357 void peck_fft_cleanup(void) {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
358 /* nothing needed any more */
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
359 }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
360
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
361 int peck_fft_next_fast_size(int n) {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
362 while (1) {
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
363 int m = n;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
364 while ((m % 2) == 0) m /= 2;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
365 while ((m % 3) == 0) m /= 3;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
366 while ((m % 5) == 0) m /= 5;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
367 if (m <= 1)
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
368 break; /* n is completely factorable by twos, threes, and fives */
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
369 n++;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
370 }
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
371 return n;
Peter Meerwald <p.meerwald@bct-electronic.com>
parents:
diff changeset
372 }

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