1 /*
2 * This file is part of libfftpack.
3 *
4 * libfftpack is free software; you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation; either version 2 of the License, or
7 * (at your option) any later version.
8 *
9 * libfftpack is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with libfftpack; if not, write to the Free Software
16 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
17 */
18
19 /*
20 * libfftpack is being developed at the Max-Planck-Institut fuer Astrophysik
21 * and financially supported by the Deutsches Zentrum fuer Luft- und Raumfahrt
22 * (DLR).
23 */
24
25 /*
26 * Copyright (C) 2005, 2006, 2007, 2008 Max-Planck-Society
27 * \author Martin Reinecke
28 */
29
30 #include <math.h>
31 #include <stdlib.h>
32 #include "fftpack.h"
33 #include "bluestein.h"
34
35 /* returns the sum of all prime factors of n */
prime_factor_sum(size_t n)36 size_t prime_factor_sum (size_t n)
37 {
38 size_t result=0,x,limit,tmp;
39 while (((tmp=(n>>1))<<1)==n)
40 { result+=2; n=tmp; }
41
42 limit=(size_t)sqrt(n+0.01);
43 for (x=3; x<=limit; x+=2)
44 while ((tmp=(n/x))*x==n)
45 {
46 result+=x;
47 n=tmp;
48 limit=(size_t)sqrt(n+0.01);
49 }
50 if (n>1) result+=n;
51
52 return result;
53 }
54
55 /* returns the smallest composite of 2, 3 and 5 which is >= n */
good_size(size_t n)56 static size_t good_size(size_t n)
57 {
58 size_t f2, f23, f235, bestfac=2*n;
59 if (n<=6) return n;
60
61 for (f2=1; f2<bestfac; f2*=2)
62 for (f23=f2; f23<bestfac; f23*=3)
63 for (f235=f23; f235<bestfac; f235*=5)
64 if (f235>=n) bestfac=f235;
65 return bestfac;
66 }
67
bluestein_i(size_t n,double ** tstorage,size_t * worksize)68 void bluestein_i (size_t n, double **tstorage, size_t *worksize)
69 {
70 static const double pi=3.14159265358979323846;
71 size_t n2=good_size(n*2-1);
72 size_t m, coeff;
73 double angle, xn2;
74 double *bk, *bkf, *work;
75 double pibyn=pi/n;
76 *worksize=2+2*n+8*n2+16;
77 *tstorage = RALLOC(double,2+2*n+8*n2+16);
78 ((size_t *)(*tstorage))[0]=n2;
79 bk = *tstorage+2;
80 bkf = *tstorage+2+2*n;
81 work= *tstorage+2+2*(n+n2);
82
83 /* initialize b_k */
84 bk[0] = 1;
85 bk[1] = 0;
86
87 coeff=0;
88 for (m=1; m<n; ++m)
89 {
90 coeff+=2*m-1;
91 if (coeff>=2*n) coeff-=2*n;
92 angle = pibyn*coeff;
93 bk[2*m] = cos(angle);
94 bk[2*m+1] = sin(angle);
95 }
96
97 /* initialize the zero-padded, Fourier transformed b_k. Add normalisation. */
98 xn2 = 1./n2;
99 bkf[0] = bk[0]*xn2;
100 bkf[1] = bk[1]*xn2;
101 for (m=2; m<2*n; m+=2)
102 {
103 bkf[m] = bkf[2*n2-m] = bk[m] *xn2;
104 bkf[m+1] = bkf[2*n2-m+1] = bk[m+1] *xn2;
105 }
106 for (m=2*n;m<=(2*n2-2*n+1);++m)
107 bkf[m]=0.;
108 cffti (n2,work);
109 cfftf (n2,bkf,work);
110 }
111
bluestein(size_t n,double * data,double * tstorage,int isign)112 void bluestein (size_t n, double *data, double *tstorage, int isign)
113 {
114 size_t n2=*((size_t *)tstorage);
115 size_t m;
116 double *bk, *bkf, *akf, *work;
117 bk = tstorage+2;
118 bkf = tstorage+2+2*n;
119 work= tstorage+2+2*(n+n2);
120 akf = tstorage+2+2*n+6*n2+16;
121
122 /* initialize a_k and FFT it */
123 if (isign>0)
124 for (m=0; m<2*n; m+=2)
125 {
126 akf[m] = data[m]*bk[m] - data[m+1]*bk[m+1];
127 akf[m+1] = data[m]*bk[m+1] + data[m+1]*bk[m];
128 }
129 else
130 for (m=0; m<2*n; m+=2)
131 {
132 akf[m] = data[m]*bk[m] + data[m+1]*bk[m+1];
133 akf[m+1] =-data[m]*bk[m+1] + data[m+1]*bk[m];
134 }
135 for (m=2*n; m<2*n2; ++m)
136 akf[m]=0;
137
138 cfftf (n2,akf,work);
139
140 /* do the convolution */
141 if (isign>0)
142 for (m=0; m<2*n2; m+=2)
143 {
144 double im = -akf[m]*bkf[m+1] + akf[m+1]*bkf[m];
145 akf[m ] = akf[m]*bkf[m] + akf[m+1]*bkf[m+1];
146 akf[m+1] = im;
147 }
148 else
149 for (m=0; m<2*n2; m+=2)
150 {
151 double im = akf[m]*bkf[m+1] + akf[m+1]*bkf[m];
152 akf[m ] = akf[m]*bkf[m] - akf[m+1]*bkf[m+1];
153 akf[m+1] = im;
154 }
155
156
157 /* inverse FFT */
158 cfftb (n2,akf,work);
159
160 /* multiply by b_k* */
161 if (isign>0)
162 for (m=0; m<2*n; m+=2)
163 {
164 data[m] = bk[m] *akf[m] - bk[m+1]*akf[m+1];
165 data[m+1] = bk[m+1]*akf[m] + bk[m] *akf[m+1];
166 }
167 else
168 for (m=0; m<2*n; m+=2)
169 {
170 data[m] = bk[m] *akf[m] + bk[m+1]*akf[m+1];
171 data[m+1] =-bk[m+1]*akf[m] + bk[m] *akf[m+1];
172 }
173 }
174