54 #include "mayer_fft.h"
64 #if defined(GOOD_TRIG)
65 #define FHT_SWAP(a,b,t) {(t)=(a);(a)=(b);(b)=(t);}
68 #define TRIG_INIT(k,c,s) \
71 for (i=2 ; i<=k ; i++) \
72 {coswrk[i]=costab[i];sinwrk[i]=sintab[i];} \
77 #define TRIG_NEXT(k,c,s) \
81 for (i=0 ; !((1<<i)&t_lam) ; i++); \
87 for (j=k-i+2 ; (1<<j)&t_lam ; j++); \
89 sinwrk[i] = halsec[i] * (sinwrk[i-1] + sinwrk[j]); \
90 coswrk[i] = halsec[i] * (coswrk[i-1] + coswrk[j]); \
93 #define TRIG_RESET(k,c,s)
96 #if defined(FAST_TRIG)
99 #define TRIG_INIT(k,c,s) \
106 #define TRIG_NEXT(k,c,s) \
112 #define TRIG_RESET(k,c,s)
115 static REAL halsec[20]=
119 .54119610014619698439972320536638942006107206337801,
120 .50979557910415916894193980398784391368261849190893,
121 .50241928618815570551167011928012092247859337193963,
122 .50060299823519630134550410676638239611758632599591,
123 .50015063602065098821477101271097658495974913010340,
124 .50003765191554772296778139077905492847503165398345,
125 .50000941253588775676512870469186533538523133757983,
126 .50000235310628608051401267171204408939326297376426,
127 .50000058827484117879868526730916804925780637276181,
128 .50000014706860214875463798283871198206179118093251,
129 .50000003676714377807315864400643020315103490883972,
130 .50000000919178552207366560348853455333939112569380,
131 .50000000229794635411562887767906868558991922348920,
132 .50000000057448658687873302235147272458812263401372
134 static REAL costab[20]=
136 .00000000000000000000000000000000000000000000000000,
137 .70710678118654752440084436210484903928483593768847,
138 .92387953251128675612818318939678828682241662586364,
139 .98078528040323044912618223613423903697393373089333,
140 .99518472667219688624483695310947992157547486872985,
141 .99879545620517239271477160475910069444320361470461,
142 .99969881869620422011576564966617219685006108125772,
143 .99992470183914454092164649119638322435060646880221,
144 .99998117528260114265699043772856771617391725094433,
145 .99999529380957617151158012570011989955298763362218,
146 .99999882345170190992902571017152601904826792288976,
147 .99999970586288221916022821773876567711626389934930,
148 .99999992646571785114473148070738785694820115568892,
149 .99999998161642929380834691540290971450507605124278,
150 .99999999540410731289097193313960614895889430318945,
151 .99999999885102682756267330779455410840053741619428
153 static REAL sintab[20]=
155 1.0000000000000000000000000000000000000000000000000,
156 .70710678118654752440084436210484903928483593768846,
157 .38268343236508977172845998403039886676134456248561,
158 .19509032201612826784828486847702224092769161775195,
159 .09801714032956060199419556388864184586113667316749,
160 .04906767432741801425495497694268265831474536302574,
161 .02454122852291228803173452945928292506546611923944,
162 .01227153828571992607940826195100321214037231959176,
163 .00613588464915447535964023459037258091705788631738,
164 .00306795676296597627014536549091984251894461021344,
165 .00153398018628476561230369715026407907995486457522,
166 .00076699031874270452693856835794857664314091945205,
167 .00038349518757139558907246168118138126339502603495,
168 .00019174759731070330743990956198900093346887403385,
169 .00009587379909597734587051721097647635118706561284,
170 .00004793689960306688454900399049465887274686668768
172 static REAL coswrk[20]=
174 .00000000000000000000000000000000000000000000000000,
175 .70710678118654752440084436210484903928483593768847,
176 .92387953251128675612818318939678828682241662586364,
177 .98078528040323044912618223613423903697393373089333,
178 .99518472667219688624483695310947992157547486872985,
179 .99879545620517239271477160475910069444320361470461,
180 .99969881869620422011576564966617219685006108125772,
181 .99992470183914454092164649119638322435060646880221,
182 .99998117528260114265699043772856771617391725094433,
183 .99999529380957617151158012570011989955298763362218,
184 .99999882345170190992902571017152601904826792288976,
185 .99999970586288221916022821773876567711626389934930,
186 .99999992646571785114473148070738785694820115568892,
187 .99999998161642929380834691540290971450507605124278,
188 .99999999540410731289097193313960614895889430318945,
189 .99999999885102682756267330779455410840053741619428
191 static REAL sinwrk[20]=
193 1.0000000000000000000000000000000000000000000000000,
194 .70710678118654752440084436210484903928483593768846,
195 .38268343236508977172845998403039886676134456248561,
196 .19509032201612826784828486847702224092769161775195,
197 .09801714032956060199419556388864184586113667316749,
198 .04906767432741801425495497694268265831474536302574,
199 .02454122852291228803173452945928292506546611923944,
200 .01227153828571992607940826195100321214037231959176,
201 .00613588464915447535964023459037258091705788631738,
202 .00306795676296597627014536549091984251894461021344,
203 .00153398018628476561230369715026407907995486457522,
204 .00076699031874270452693856835794857664314091945205,
205 .00038349518757139558907246168118138126339502603495,
206 .00019174759731070330743990956198900093346887403385,
207 .00009587379909597734587051721097647635118706561284,
208 .00004793689960306688454900399049465887274686668768
212 #define SQRT2_2 0.70710678118654752440084436210484
213 #define SQRT2 2*0.70710678118654752440084436210484
215 void mayer_fht(REAL *fz,
int n)
220 int k,k1,k2,k3,k4,kx;
224 for (k1=1,k2=0;k1<n;k1++)
227 for (k=n>>1; (!((k2^=k)&k)); k>>=1);
230 aa=fz[k1];fz[k1]=fz[k2];fz[k2]=aa;
233 for ( k=0 ; (1<<k)<n ; k++ );
237 for (fi=fz,fn=fz+n;fi<fn;fi+=4)
252 for (fi=fz,fn=fz+n,gi=fi+1;fi<fn;fi+=8,gi+=8)
254 REAL bs1,bc1,bs2,bc2,bs3,bc3,bs4,bc4,
255 bg0,bf0,bf1,bg1,bf2,bg2,bf3,bg3;
256 bc1 = fi[0 ] - gi[0 ];
257 bs1 = fi[0 ] + gi[0 ];
258 bc2 = fi[2 ] - gi[2 ];
259 bs2 = fi[2 ] + gi[2 ];
260 bc3 = fi[4 ] - gi[4 ];
261 bs3 = fi[4 ] + gi[4 ];
262 bc4 = fi[6 ] - gi[6 ];
263 bs4 = fi[6 ] + gi[6 ];
299 REAL g0,f0,f1,g1,f2,g2,f3,g3;
300 f1 = fi[0 ] - fi[k1];
301 f0 = fi[0 ] + fi[k1];
302 f3 = fi[k2] - fi[k3];
303 f2 = fi[k2] + fi[k3];
308 g1 = gi[0 ] - gi[k1];
309 g0 = gi[0 ] + gi[k1];
320 for (ii=1;ii<kx;ii++)
331 REAL a,b,g0,f0,f1,g1,f2,g2,f3,g3;
332 b = s2*fi[k1] - c2*gi[k1];
333 a = c2*fi[k1] + s2*gi[k1];
338 b = s2*fi[k3] - c2*gi[k3];
339 a = c2*fi[k3] + s2*gi[k3];
364 void mayer_fft(
int n, REAL *real, REAL *imag)
369 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
370 a = real[i]; b = real[j]; q=a+b; r=a-b;
371 c = imag[i]; d = imag[j]; s=c+d; t=c-d;
372 real[i] = (q+t)*.5; real[j] = (q-t)*.5;
373 imag[i] = (s-r)*.5; imag[j] = (s+r)*.5;
379 void mayer_ifft(
int n, REAL *real, REAL *imag)
386 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
387 a = real[i]; b = real[j]; q=a+b; r=a-b;
388 c = imag[i]; d = imag[j]; s=c+d; t=c-d;
389 imag[i] = (s+r)*0.5; imag[j] = (s-r)*0.5;
390 real[i] = (q-t)*0.5; real[j] = (q+t)*0.5;
394 void mayer_realfft(
int n, REAL *real)
400 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {
408 void mayer_realifft(
int n, REAL *real)
413 for (i=1,j=n-1,k=n/2;i<k;i++,j--) {