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