diff options
author | Sean Bartell <wingedtachikoma@gmail.com> | 2011-06-25 21:32:25 -0400 |
---|---|---|
committer | Nils Wallménius <nils@rockbox.org> | 2012-04-25 22:13:20 +0200 |
commit | f40bfc9267b13b54e6379dfe7539447662879d24 (patch) | |
tree | 9b20069d5e62809ff434061ad730096836f916f2 /lib/rbcodec/codecs/libspeex/smallft.c | |
parent | a0009907de7a0107d49040d8a180f140e2eff299 (diff) | |
download | rockbox-f40bfc9267b13b54e6379dfe7539447662879d24.tar.gz rockbox-f40bfc9267b13b54e6379dfe7539447662879d24.zip |
Add codecs to librbcodec.
Change-Id: Id7f4717d51ed02d67cb9f9cb3c0ada4a81843f97
Reviewed-on: http://gerrit.rockbox.org/137
Reviewed-by: Nils Wallménius <nils@rockbox.org>
Tested-by: Nils Wallménius <nils@rockbox.org>
Diffstat (limited to 'lib/rbcodec/codecs/libspeex/smallft.c')
-rw-r--r-- | lib/rbcodec/codecs/libspeex/smallft.c | 1261 |
1 files changed, 1261 insertions, 0 deletions
diff --git a/lib/rbcodec/codecs/libspeex/smallft.c b/lib/rbcodec/codecs/libspeex/smallft.c new file mode 100644 index 0000000000..6e3a927b5f --- /dev/null +++ b/lib/rbcodec/codecs/libspeex/smallft.c | |||
@@ -0,0 +1,1261 @@ | |||
1 | /******************************************************************** | ||
2 | * * | ||
3 | * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. * | ||
4 | * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS * | ||
5 | * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE * | ||
6 | * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. * | ||
7 | * * | ||
8 | * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001 * | ||
9 | * by the XIPHOPHORUS Company http://www.xiph.org/ * | ||
10 | * * | ||
11 | ******************************************************************** | ||
12 | |||
13 | function: *unnormalized* fft transform | ||
14 | last mod: $Id$ | ||
15 | |||
16 | ********************************************************************/ | ||
17 | |||
18 | /* FFT implementation from OggSquish, minus cosine transforms, | ||
19 | * minus all but radix 2/4 case. In Vorbis we only need this | ||
20 | * cut-down version. | ||
21 | * | ||
22 | * To do more than just power-of-two sized vectors, see the full | ||
23 | * version I wrote for NetLib. | ||
24 | * | ||
25 | * Note that the packing is a little strange; rather than the FFT r/i | ||
26 | * packing following R_0, I_n, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1, | ||
27 | * it follows R_0, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1, I_n like the | ||
28 | * FORTRAN version | ||
29 | */ | ||
30 | |||
31 | #ifdef HAVE_CONFIG_H | ||
32 | #include "config-speex.h" | ||
33 | #endif | ||
34 | |||
35 | #include <math.h> | ||
36 | #include "smallft.h" | ||
37 | #include "arch.h" | ||
38 | #include "os_support.h" | ||
39 | |||
40 | static void drfti1(int n, float *wa, int *ifac){ | ||
41 | static int ntryh[4] = { 4,2,3,5 }; | ||
42 | static float tpi = 6.28318530717958648f; | ||
43 | float arg,argh,argld,fi; | ||
44 | int ntry=0,i,j=-1; | ||
45 | int k1, l1, l2, ib; | ||
46 | int ld, ii, ip, is, nq, nr; | ||
47 | int ido, ipm, nfm1; | ||
48 | int nl=n; | ||
49 | int nf=0; | ||
50 | |||
51 | L101: | ||
52 | j++; | ||
53 | if (j < 4) | ||
54 | ntry=ntryh[j]; | ||
55 | else | ||
56 | ntry+=2; | ||
57 | |||
58 | L104: | ||
59 | nq=nl/ntry; | ||
60 | nr=nl-ntry*nq; | ||
61 | if (nr!=0) goto L101; | ||
62 | |||
63 | nf++; | ||
64 | ifac[nf+1]=ntry; | ||
65 | nl=nq; | ||
66 | if(ntry!=2)goto L107; | ||
67 | if(nf==1)goto L107; | ||
68 | |||
69 | for (i=1;i<nf;i++){ | ||
70 | ib=nf-i+1; | ||
71 | ifac[ib+1]=ifac[ib]; | ||
72 | } | ||
73 | ifac[2] = 2; | ||
74 | |||
75 | L107: | ||
76 | if(nl!=1)goto L104; | ||
77 | ifac[0]=n; | ||
78 | ifac[1]=nf; | ||
79 | argh=tpi/n; | ||
80 | is=0; | ||
81 | nfm1=nf-1; | ||
82 | l1=1; | ||
83 | |||
84 | if(nfm1==0)return; | ||
85 | |||
86 | for (k1=0;k1<nfm1;k1++){ | ||
87 | ip=ifac[k1+2]; | ||
88 | ld=0; | ||
89 | l2=l1*ip; | ||
90 | ido=n/l2; | ||
91 | ipm=ip-1; | ||
92 | |||
93 | for (j=0;j<ipm;j++){ | ||
94 | ld+=l1; | ||
95 | i=is; | ||
96 | argld=(float)ld*argh; | ||
97 | fi=0.f; | ||
98 | for (ii=2;ii<ido;ii+=2){ | ||
99 | fi+=1.f; | ||
100 | arg=fi*argld; | ||
101 | wa[i++]=cos(arg); | ||
102 | wa[i++]=sin(arg); | ||
103 | } | ||
104 | is+=ido; | ||
105 | } | ||
106 | l1=l2; | ||
107 | } | ||
108 | } | ||
109 | |||
110 | static void fdrffti(int n, float *wsave, int *ifac){ | ||
111 | |||
112 | if (n == 1) return; | ||
113 | drfti1(n, wsave+n, ifac); | ||
114 | } | ||
115 | |||
116 | static void dradf2(int ido,int l1,float *cc,float *ch,float *wa1){ | ||
117 | int i,k; | ||
118 | float ti2,tr2; | ||
119 | int t0,t1,t2,t3,t4,t5,t6; | ||
120 | |||
121 | t1=0; | ||
122 | t0=(t2=l1*ido); | ||
123 | t3=ido<<1; | ||
124 | for(k=0;k<l1;k++){ | ||
125 | ch[t1<<1]=cc[t1]+cc[t2]; | ||
126 | ch[(t1<<1)+t3-1]=cc[t1]-cc[t2]; | ||
127 | t1+=ido; | ||
128 | t2+=ido; | ||
129 | } | ||
130 | |||
131 | if(ido<2)return; | ||
132 | if(ido==2)goto L105; | ||
133 | |||
134 | t1=0; | ||
135 | t2=t0; | ||
136 | for(k=0;k<l1;k++){ | ||
137 | t3=t2; | ||
138 | t4=(t1<<1)+(ido<<1); | ||
139 | t5=t1; | ||
140 | t6=t1+t1; | ||
141 | for(i=2;i<ido;i+=2){ | ||
142 | t3+=2; | ||
143 | t4-=2; | ||
144 | t5+=2; | ||
145 | t6+=2; | ||
146 | tr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3]; | ||
147 | ti2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1]; | ||
148 | ch[t6]=cc[t5]+ti2; | ||
149 | ch[t4]=ti2-cc[t5]; | ||
150 | ch[t6-1]=cc[t5-1]+tr2; | ||
151 | ch[t4-1]=cc[t5-1]-tr2; | ||
152 | } | ||
153 | t1+=ido; | ||
154 | t2+=ido; | ||
155 | } | ||
156 | |||
157 | if(ido%2==1)return; | ||
158 | |||
159 | L105: | ||
160 | t3=(t2=(t1=ido)-1); | ||
161 | t2+=t0; | ||
162 | for(k=0;k<l1;k++){ | ||
163 | ch[t1]=-cc[t2]; | ||
164 | ch[t1-1]=cc[t3]; | ||
165 | t1+=ido<<1; | ||
166 | t2+=ido; | ||
167 | t3+=ido; | ||
168 | } | ||
169 | } | ||
170 | |||
171 | static void dradf4(int ido,int l1,float *cc,float *ch,float *wa1, | ||
172 | float *wa2,float *wa3){ | ||
173 | static float hsqt2 = .70710678118654752f; | ||
174 | int i,k,t0,t1,t2,t3,t4,t5,t6; | ||
175 | float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4; | ||
176 | t0=l1*ido; | ||
177 | |||
178 | t1=t0; | ||
179 | t4=t1<<1; | ||
180 | t2=t1+(t1<<1); | ||
181 | t3=0; | ||
182 | |||
183 | for(k=0;k<l1;k++){ | ||
184 | tr1=cc[t1]+cc[t2]; | ||
185 | tr2=cc[t3]+cc[t4]; | ||
186 | |||
187 | ch[t5=t3<<2]=tr1+tr2; | ||
188 | ch[(ido<<2)+t5-1]=tr2-tr1; | ||
189 | ch[(t5+=(ido<<1))-1]=cc[t3]-cc[t4]; | ||
190 | ch[t5]=cc[t2]-cc[t1]; | ||
191 | |||
192 | t1+=ido; | ||
193 | t2+=ido; | ||
194 | t3+=ido; | ||
195 | t4+=ido; | ||
196 | } | ||
197 | |||
198 | if(ido<2)return; | ||
199 | if(ido==2)goto L105; | ||
200 | |||
201 | |||
202 | t1=0; | ||
203 | for(k=0;k<l1;k++){ | ||
204 | t2=t1; | ||
205 | t4=t1<<2; | ||
206 | t5=(t6=ido<<1)+t4; | ||
207 | for(i=2;i<ido;i+=2){ | ||
208 | t3=(t2+=2); | ||
209 | t4+=2; | ||
210 | t5-=2; | ||
211 | |||
212 | t3+=t0; | ||
213 | cr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3]; | ||
214 | ci2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1]; | ||
215 | t3+=t0; | ||
216 | cr3=wa2[i-2]*cc[t3-1]+wa2[i-1]*cc[t3]; | ||
217 | ci3=wa2[i-2]*cc[t3]-wa2[i-1]*cc[t3-1]; | ||
218 | t3+=t0; | ||
219 | cr4=wa3[i-2]*cc[t3-1]+wa3[i-1]*cc[t3]; | ||
220 | ci4=wa3[i-2]*cc[t3]-wa3[i-1]*cc[t3-1]; | ||
221 | |||
222 | tr1=cr2+cr4; | ||
223 | tr4=cr4-cr2; | ||
224 | ti1=ci2+ci4; | ||
225 | ti4=ci2-ci4; | ||
226 | |||
227 | ti2=cc[t2]+ci3; | ||
228 | ti3=cc[t2]-ci3; | ||
229 | tr2=cc[t2-1]+cr3; | ||
230 | tr3=cc[t2-1]-cr3; | ||
231 | |||
232 | ch[t4-1]=tr1+tr2; | ||
233 | ch[t4]=ti1+ti2; | ||
234 | |||
235 | ch[t5-1]=tr3-ti4; | ||
236 | ch[t5]=tr4-ti3; | ||
237 | |||
238 | ch[t4+t6-1]=ti4+tr3; | ||
239 | ch[t4+t6]=tr4+ti3; | ||
240 | |||
241 | ch[t5+t6-1]=tr2-tr1; | ||
242 | ch[t5+t6]=ti1-ti2; | ||
243 | } | ||
244 | t1+=ido; | ||
245 | } | ||
246 | if(ido&1)return; | ||
247 | |||
248 | L105: | ||
249 | |||
250 | t2=(t1=t0+ido-1)+(t0<<1); | ||
251 | t3=ido<<2; | ||
252 | t4=ido; | ||
253 | t5=ido<<1; | ||
254 | t6=ido; | ||
255 | |||
256 | for(k=0;k<l1;k++){ | ||
257 | ti1=-hsqt2*(cc[t1]+cc[t2]); | ||
258 | tr1=hsqt2*(cc[t1]-cc[t2]); | ||
259 | |||
260 | ch[t4-1]=tr1+cc[t6-1]; | ||
261 | ch[t4+t5-1]=cc[t6-1]-tr1; | ||
262 | |||
263 | ch[t4]=ti1-cc[t1+t0]; | ||
264 | ch[t4+t5]=ti1+cc[t1+t0]; | ||
265 | |||
266 | t1+=ido; | ||
267 | t2+=ido; | ||
268 | t4+=t3; | ||
269 | t6+=ido; | ||
270 | } | ||
271 | } | ||
272 | |||
273 | static void dradfg(int ido,int ip,int l1,int idl1,float *cc,float *c1, | ||
274 | float *c2,float *ch,float *ch2,float *wa){ | ||
275 | |||
276 | static float tpi=6.283185307179586f; | ||
277 | int idij,ipph,i,j,k,l,ic,ik,is; | ||
278 | int t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10; | ||
279 | float dc2,ai1,ai2,ar1,ar2,ds2; | ||
280 | int nbd; | ||
281 | float dcp,arg,dsp,ar1h,ar2h; | ||
282 | int idp2,ipp2; | ||
283 | |||
284 | arg=tpi/(float)ip; | ||
285 | dcp=cos(arg); | ||
286 | dsp=sin(arg); | ||
287 | ipph=(ip+1)>>1; | ||
288 | ipp2=ip; | ||
289 | idp2=ido; | ||
290 | nbd=(ido-1)>>1; | ||
291 | t0=l1*ido; | ||
292 | t10=ip*ido; | ||
293 | |||
294 | if(ido==1)goto L119; | ||
295 | for(ik=0;ik<idl1;ik++)ch2[ik]=c2[ik]; | ||
296 | |||
297 | t1=0; | ||
298 | for(j=1;j<ip;j++){ | ||
299 | t1+=t0; | ||
300 | t2=t1; | ||
301 | for(k=0;k<l1;k++){ | ||
302 | ch[t2]=c1[t2]; | ||
303 | t2+=ido; | ||
304 | } | ||
305 | } | ||
306 | |||
307 | is=-ido; | ||
308 | t1=0; | ||
309 | if(nbd>l1){ | ||
310 | for(j=1;j<ip;j++){ | ||
311 | t1+=t0; | ||
312 | is+=ido; | ||
313 | t2= -ido+t1; | ||
314 | for(k=0;k<l1;k++){ | ||
315 | idij=is-1; | ||
316 | t2+=ido; | ||
317 | t3=t2; | ||
318 | for(i=2;i<ido;i+=2){ | ||
319 | idij+=2; | ||
320 | t3+=2; | ||
321 | ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3]; | ||
322 | ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1]; | ||
323 | } | ||
324 | } | ||
325 | } | ||
326 | }else{ | ||
327 | |||
328 | for(j=1;j<ip;j++){ | ||
329 | is+=ido; | ||
330 | idij=is-1; | ||
331 | t1+=t0; | ||
332 | t2=t1; | ||
333 | for(i=2;i<ido;i+=2){ | ||
334 | idij+=2; | ||
335 | t2+=2; | ||
336 | t3=t2; | ||
337 | for(k=0;k<l1;k++){ | ||
338 | ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3]; | ||
339 | ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1]; | ||
340 | t3+=ido; | ||
341 | } | ||
342 | } | ||
343 | } | ||
344 | } | ||
345 | |||
346 | t1=0; | ||
347 | t2=ipp2*t0; | ||
348 | if(nbd<l1){ | ||
349 | for(j=1;j<ipph;j++){ | ||
350 | t1+=t0; | ||
351 | t2-=t0; | ||
352 | t3=t1; | ||
353 | t4=t2; | ||
354 | for(i=2;i<ido;i+=2){ | ||
355 | t3+=2; | ||
356 | t4+=2; | ||
357 | t5=t3-ido; | ||
358 | t6=t4-ido; | ||
359 | for(k=0;k<l1;k++){ | ||
360 | t5+=ido; | ||
361 | t6+=ido; | ||
362 | c1[t5-1]=ch[t5-1]+ch[t6-1]; | ||
363 | c1[t6-1]=ch[t5]-ch[t6]; | ||
364 | c1[t5]=ch[t5]+ch[t6]; | ||
365 | c1[t6]=ch[t6-1]-ch[t5-1]; | ||
366 | } | ||
367 | } | ||
368 | } | ||
369 | }else{ | ||
370 | for(j=1;j<ipph;j++){ | ||
371 | t1+=t0; | ||
372 | t2-=t0; | ||
373 | t3=t1; | ||
374 | t4=t2; | ||
375 | for(k=0;k<l1;k++){ | ||
376 | t5=t3; | ||
377 | t6=t4; | ||
378 | for(i=2;i<ido;i+=2){ | ||
379 | t5+=2; | ||
380 | t6+=2; | ||
381 | c1[t5-1]=ch[t5-1]+ch[t6-1]; | ||
382 | c1[t6-1]=ch[t5]-ch[t6]; | ||
383 | c1[t5]=ch[t5]+ch[t6]; | ||
384 | c1[t6]=ch[t6-1]-ch[t5-1]; | ||
385 | } | ||
386 | t3+=ido; | ||
387 | t4+=ido; | ||
388 | } | ||
389 | } | ||
390 | } | ||
391 | |||
392 | L119: | ||
393 | for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik]; | ||
394 | |||
395 | t1=0; | ||
396 | t2=ipp2*idl1; | ||
397 | for(j=1;j<ipph;j++){ | ||
398 | t1+=t0; | ||
399 | t2-=t0; | ||
400 | t3=t1-ido; | ||
401 | t4=t2-ido; | ||
402 | for(k=0;k<l1;k++){ | ||
403 | t3+=ido; | ||
404 | t4+=ido; | ||
405 | c1[t3]=ch[t3]+ch[t4]; | ||
406 | c1[t4]=ch[t4]-ch[t3]; | ||
407 | } | ||
408 | } | ||
409 | |||
410 | ar1=1.f; | ||
411 | ai1=0.f; | ||
412 | t1=0; | ||
413 | t2=ipp2*idl1; | ||
414 | t3=(ip-1)*idl1; | ||
415 | for(l=1;l<ipph;l++){ | ||
416 | t1+=idl1; | ||
417 | t2-=idl1; | ||
418 | ar1h=dcp*ar1-dsp*ai1; | ||
419 | ai1=dcp*ai1+dsp*ar1; | ||
420 | ar1=ar1h; | ||
421 | t4=t1; | ||
422 | t5=t2; | ||
423 | t6=t3; | ||
424 | t7=idl1; | ||
425 | |||
426 | for(ik=0;ik<idl1;ik++){ | ||
427 | ch2[t4++]=c2[ik]+ar1*c2[t7++]; | ||
428 | ch2[t5++]=ai1*c2[t6++]; | ||
429 | } | ||
430 | |||
431 | dc2=ar1; | ||
432 | ds2=ai1; | ||
433 | ar2=ar1; | ||
434 | ai2=ai1; | ||
435 | |||
436 | t4=idl1; | ||
437 | t5=(ipp2-1)*idl1; | ||
438 | for(j=2;j<ipph;j++){ | ||
439 | t4+=idl1; | ||
440 | t5-=idl1; | ||
441 | |||
442 | ar2h=dc2*ar2-ds2*ai2; | ||
443 | ai2=dc2*ai2+ds2*ar2; | ||
444 | ar2=ar2h; | ||
445 | |||
446 | t6=t1; | ||
447 | t7=t2; | ||
448 | t8=t4; | ||
449 | t9=t5; | ||
450 | for(ik=0;ik<idl1;ik++){ | ||
451 | ch2[t6++]+=ar2*c2[t8++]; | ||
452 | ch2[t7++]+=ai2*c2[t9++]; | ||
453 | } | ||
454 | } | ||
455 | } | ||
456 | |||
457 | t1=0; | ||
458 | for(j=1;j<ipph;j++){ | ||
459 | t1+=idl1; | ||
460 | t2=t1; | ||
461 | for(ik=0;ik<idl1;ik++)ch2[ik]+=c2[t2++]; | ||
462 | } | ||
463 | |||
464 | if(ido<l1)goto L132; | ||
465 | |||
466 | t1=0; | ||
467 | t2=0; | ||
468 | for(k=0;k<l1;k++){ | ||
469 | t3=t1; | ||
470 | t4=t2; | ||
471 | for(i=0;i<ido;i++)cc[t4++]=ch[t3++]; | ||
472 | t1+=ido; | ||
473 | t2+=t10; | ||
474 | } | ||
475 | |||
476 | goto L135; | ||
477 | |||
478 | L132: | ||
479 | for(i=0;i<ido;i++){ | ||
480 | t1=i; | ||
481 | t2=i; | ||
482 | for(k=0;k<l1;k++){ | ||
483 | cc[t2]=ch[t1]; | ||
484 | t1+=ido; | ||
485 | t2+=t10; | ||
486 | } | ||
487 | } | ||
488 | |||
489 | L135: | ||
490 | t1=0; | ||
491 | t2=ido<<1; | ||
492 | t3=0; | ||
493 | t4=ipp2*t0; | ||
494 | for(j=1;j<ipph;j++){ | ||
495 | |||
496 | t1+=t2; | ||
497 | t3+=t0; | ||
498 | t4-=t0; | ||
499 | |||
500 | t5=t1; | ||
501 | t6=t3; | ||
502 | t7=t4; | ||
503 | |||
504 | for(k=0;k<l1;k++){ | ||
505 | cc[t5-1]=ch[t6]; | ||
506 | cc[t5]=ch[t7]; | ||
507 | t5+=t10; | ||
508 | t6+=ido; | ||
509 | t7+=ido; | ||
510 | } | ||
511 | } | ||
512 | |||
513 | if(ido==1)return; | ||
514 | if(nbd<l1)goto L141; | ||
515 | |||
516 | t1=-ido; | ||
517 | t3=0; | ||
518 | t4=0; | ||
519 | t5=ipp2*t0; | ||
520 | for(j=1;j<ipph;j++){ | ||
521 | t1+=t2; | ||
522 | t3+=t2; | ||
523 | t4+=t0; | ||
524 | t5-=t0; | ||
525 | t6=t1; | ||
526 | t7=t3; | ||
527 | t8=t4; | ||
528 | t9=t5; | ||
529 | for(k=0;k<l1;k++){ | ||
530 | for(i=2;i<ido;i+=2){ | ||
531 | ic=idp2-i; | ||
532 | cc[i+t7-1]=ch[i+t8-1]+ch[i+t9-1]; | ||
533 | cc[ic+t6-1]=ch[i+t8-1]-ch[i+t9-1]; | ||
534 | cc[i+t7]=ch[i+t8]+ch[i+t9]; | ||
535 | cc[ic+t6]=ch[i+t9]-ch[i+t8]; | ||
536 | } | ||
537 | t6+=t10; | ||
538 | t7+=t10; | ||
539 | t8+=ido; | ||
540 | t9+=ido; | ||
541 | } | ||
542 | } | ||
543 | return; | ||
544 | |||
545 | L141: | ||
546 | |||
547 | t1=-ido; | ||
548 | t3=0; | ||
549 | t4=0; | ||
550 | t5=ipp2*t0; | ||
551 | for(j=1;j<ipph;j++){ | ||
552 | t1+=t2; | ||
553 | t3+=t2; | ||
554 | t4+=t0; | ||
555 | t5-=t0; | ||
556 | for(i=2;i<ido;i+=2){ | ||
557 | t6=idp2+t1-i; | ||
558 | t7=i+t3; | ||
559 | t8=i+t4; | ||
560 | t9=i+t5; | ||
561 | for(k=0;k<l1;k++){ | ||
562 | cc[t7-1]=ch[t8-1]+ch[t9-1]; | ||
563 | cc[t6-1]=ch[t8-1]-ch[t9-1]; | ||
564 | cc[t7]=ch[t8]+ch[t9]; | ||
565 | cc[t6]=ch[t9]-ch[t8]; | ||
566 | t6+=t10; | ||
567 | t7+=t10; | ||
568 | t8+=ido; | ||
569 | t9+=ido; | ||
570 | } | ||
571 | } | ||
572 | } | ||
573 | } | ||
574 | |||
575 | static void drftf1(int n,float *c,float *ch,float *wa,int *ifac){ | ||
576 | int i,k1,l1,l2; | ||
577 | int na,kh,nf; | ||
578 | int ip,iw,ido,idl1,ix2,ix3; | ||
579 | |||
580 | nf=ifac[1]; | ||
581 | na=1; | ||
582 | l2=n; | ||
583 | iw=n; | ||
584 | |||
585 | for(k1=0;k1<nf;k1++){ | ||
586 | kh=nf-k1; | ||
587 | ip=ifac[kh+1]; | ||
588 | l1=l2/ip; | ||
589 | ido=n/l2; | ||
590 | idl1=ido*l1; | ||
591 | iw-=(ip-1)*ido; | ||
592 | na=1-na; | ||
593 | |||
594 | if(ip!=4)goto L102; | ||
595 | |||
596 | ix2=iw+ido; | ||
597 | ix3=ix2+ido; | ||
598 | if(na!=0) | ||
599 | dradf4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1); | ||
600 | else | ||
601 | dradf4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1); | ||
602 | goto L110; | ||
603 | |||
604 | L102: | ||
605 | if(ip!=2)goto L104; | ||
606 | if(na!=0)goto L103; | ||
607 | |||
608 | dradf2(ido,l1,c,ch,wa+iw-1); | ||
609 | goto L110; | ||
610 | |||
611 | L103: | ||
612 | dradf2(ido,l1,ch,c,wa+iw-1); | ||
613 | goto L110; | ||
614 | |||
615 | L104: | ||
616 | if(ido==1)na=1-na; | ||
617 | if(na!=0)goto L109; | ||
618 | |||
619 | dradfg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1); | ||
620 | na=1; | ||
621 | goto L110; | ||
622 | |||
623 | L109: | ||
624 | dradfg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1); | ||
625 | na=0; | ||
626 | |||
627 | L110: | ||
628 | l2=l1; | ||
629 | } | ||
630 | |||
631 | if(na==1)return; | ||
632 | |||
633 | for(i=0;i<n;i++)c[i]=ch[i]; | ||
634 | } | ||
635 | |||
636 | static void dradb2(int ido,int l1,float *cc,float *ch,float *wa1){ | ||
637 | int i,k,t0,t1,t2,t3,t4,t5,t6; | ||
638 | float ti2,tr2; | ||
639 | |||
640 | t0=l1*ido; | ||
641 | |||
642 | t1=0; | ||
643 | t2=0; | ||
644 | t3=(ido<<1)-1; | ||
645 | for(k=0;k<l1;k++){ | ||
646 | ch[t1]=cc[t2]+cc[t3+t2]; | ||
647 | ch[t1+t0]=cc[t2]-cc[t3+t2]; | ||
648 | t2=(t1+=ido)<<1; | ||
649 | } | ||
650 | |||
651 | if(ido<2)return; | ||
652 | if(ido==2)goto L105; | ||
653 | |||
654 | t1=0; | ||
655 | t2=0; | ||
656 | for(k=0;k<l1;k++){ | ||
657 | t3=t1; | ||
658 | t5=(t4=t2)+(ido<<1); | ||
659 | t6=t0+t1; | ||
660 | for(i=2;i<ido;i+=2){ | ||
661 | t3+=2; | ||
662 | t4+=2; | ||
663 | t5-=2; | ||
664 | t6+=2; | ||
665 | ch[t3-1]=cc[t4-1]+cc[t5-1]; | ||
666 | tr2=cc[t4-1]-cc[t5-1]; | ||
667 | ch[t3]=cc[t4]-cc[t5]; | ||
668 | ti2=cc[t4]+cc[t5]; | ||
669 | ch[t6-1]=wa1[i-2]*tr2-wa1[i-1]*ti2; | ||
670 | ch[t6]=wa1[i-2]*ti2+wa1[i-1]*tr2; | ||
671 | } | ||
672 | t2=(t1+=ido)<<1; | ||
673 | } | ||
674 | |||
675 | if(ido%2==1)return; | ||
676 | |||
677 | L105: | ||
678 | t1=ido-1; | ||
679 | t2=ido-1; | ||
680 | for(k=0;k<l1;k++){ | ||
681 | ch[t1]=cc[t2]+cc[t2]; | ||
682 | ch[t1+t0]=-(cc[t2+1]+cc[t2+1]); | ||
683 | t1+=ido; | ||
684 | t2+=ido<<1; | ||
685 | } | ||
686 | } | ||
687 | |||
688 | static void dradb3(int ido,int l1,float *cc,float *ch,float *wa1, | ||
689 | float *wa2){ | ||
690 | static float taur = -.5f; | ||
691 | static float taui = .8660254037844386f; | ||
692 | int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10; | ||
693 | float ci2,ci3,di2,di3,cr2,cr3,dr2,dr3,ti2,tr2; | ||
694 | t0=l1*ido; | ||
695 | |||
696 | t1=0; | ||
697 | t2=t0<<1; | ||
698 | t3=ido<<1; | ||
699 | t4=ido+(ido<<1); | ||
700 | t5=0; | ||
701 | for(k=0;k<l1;k++){ | ||
702 | tr2=cc[t3-1]+cc[t3-1]; | ||
703 | cr2=cc[t5]+(taur*tr2); | ||
704 | ch[t1]=cc[t5]+tr2; | ||
705 | ci3=taui*(cc[t3]+cc[t3]); | ||
706 | ch[t1+t0]=cr2-ci3; | ||
707 | ch[t1+t2]=cr2+ci3; | ||
708 | t1+=ido; | ||
709 | t3+=t4; | ||
710 | t5+=t4; | ||
711 | } | ||
712 | |||
713 | if(ido==1)return; | ||
714 | |||
715 | t1=0; | ||
716 | t3=ido<<1; | ||
717 | for(k=0;k<l1;k++){ | ||
718 | t7=t1+(t1<<1); | ||
719 | t6=(t5=t7+t3); | ||
720 | t8=t1; | ||
721 | t10=(t9=t1+t0)+t0; | ||
722 | |||
723 | for(i=2;i<ido;i+=2){ | ||
724 | t5+=2; | ||
725 | t6-=2; | ||
726 | t7+=2; | ||
727 | t8+=2; | ||
728 | t9+=2; | ||
729 | t10+=2; | ||
730 | tr2=cc[t5-1]+cc[t6-1]; | ||
731 | cr2=cc[t7-1]+(taur*tr2); | ||
732 | ch[t8-1]=cc[t7-1]+tr2; | ||
733 | ti2=cc[t5]-cc[t6]; | ||
734 | ci2=cc[t7]+(taur*ti2); | ||
735 | ch[t8]=cc[t7]+ti2; | ||
736 | cr3=taui*(cc[t5-1]-cc[t6-1]); | ||
737 | ci3=taui*(cc[t5]+cc[t6]); | ||
738 | dr2=cr2-ci3; | ||
739 | dr3=cr2+ci3; | ||
740 | di2=ci2+cr3; | ||
741 | di3=ci2-cr3; | ||
742 | ch[t9-1]=wa1[i-2]*dr2-wa1[i-1]*di2; | ||
743 | ch[t9]=wa1[i-2]*di2+wa1[i-1]*dr2; | ||
744 | ch[t10-1]=wa2[i-2]*dr3-wa2[i-1]*di3; | ||
745 | ch[t10]=wa2[i-2]*di3+wa2[i-1]*dr3; | ||
746 | } | ||
747 | t1+=ido; | ||
748 | } | ||
749 | } | ||
750 | |||
751 | static void dradb4(int ido,int l1,float *cc,float *ch,float *wa1, | ||
752 | float *wa2,float *wa3){ | ||
753 | static float sqrt2=1.414213562373095f; | ||
754 | int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8; | ||
755 | float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4; | ||
756 | t0=l1*ido; | ||
757 | |||
758 | t1=0; | ||
759 | t2=ido<<2; | ||
760 | t3=0; | ||
761 | t6=ido<<1; | ||
762 | for(k=0;k<l1;k++){ | ||
763 | t4=t3+t6; | ||
764 | t5=t1; | ||
765 | tr3=cc[t4-1]+cc[t4-1]; | ||
766 | tr4=cc[t4]+cc[t4]; | ||
767 | tr1=cc[t3]-cc[(t4+=t6)-1]; | ||
768 | tr2=cc[t3]+cc[t4-1]; | ||
769 | ch[t5]=tr2+tr3; | ||
770 | ch[t5+=t0]=tr1-tr4; | ||
771 | ch[t5+=t0]=tr2-tr3; | ||
772 | ch[t5+=t0]=tr1+tr4; | ||
773 | t1+=ido; | ||
774 | t3+=t2; | ||
775 | } | ||
776 | |||
777 | if(ido<2)return; | ||
778 | if(ido==2)goto L105; | ||
779 | |||
780 | t1=0; | ||
781 | for(k=0;k<l1;k++){ | ||
782 | t5=(t4=(t3=(t2=t1<<2)+t6))+t6; | ||
783 | t7=t1; | ||
784 | for(i=2;i<ido;i+=2){ | ||
785 | t2+=2; | ||
786 | t3+=2; | ||
787 | t4-=2; | ||
788 | t5-=2; | ||
789 | t7+=2; | ||
790 | ti1=cc[t2]+cc[t5]; | ||
791 | ti2=cc[t2]-cc[t5]; | ||
792 | ti3=cc[t3]-cc[t4]; | ||
793 | tr4=cc[t3]+cc[t4]; | ||
794 | tr1=cc[t2-1]-cc[t5-1]; | ||
795 | tr2=cc[t2-1]+cc[t5-1]; | ||
796 | ti4=cc[t3-1]-cc[t4-1]; | ||
797 | tr3=cc[t3-1]+cc[t4-1]; | ||
798 | ch[t7-1]=tr2+tr3; | ||
799 | cr3=tr2-tr3; | ||
800 | ch[t7]=ti2+ti3; | ||
801 | ci3=ti2-ti3; | ||
802 | cr2=tr1-tr4; | ||
803 | cr4=tr1+tr4; | ||
804 | ci2=ti1+ti4; | ||
805 | ci4=ti1-ti4; | ||
806 | |||
807 | ch[(t8=t7+t0)-1]=wa1[i-2]*cr2-wa1[i-1]*ci2; | ||
808 | ch[t8]=wa1[i-2]*ci2+wa1[i-1]*cr2; | ||
809 | ch[(t8+=t0)-1]=wa2[i-2]*cr3-wa2[i-1]*ci3; | ||
810 | ch[t8]=wa2[i-2]*ci3+wa2[i-1]*cr3; | ||
811 | ch[(t8+=t0)-1]=wa3[i-2]*cr4-wa3[i-1]*ci4; | ||
812 | ch[t8]=wa3[i-2]*ci4+wa3[i-1]*cr4; | ||
813 | } | ||
814 | t1+=ido; | ||
815 | } | ||
816 | |||
817 | if(ido%2 == 1)return; | ||
818 | |||
819 | L105: | ||
820 | |||
821 | t1=ido; | ||
822 | t2=ido<<2; | ||
823 | t3=ido-1; | ||
824 | t4=ido+(ido<<1); | ||
825 | for(k=0;k<l1;k++){ | ||
826 | t5=t3; | ||
827 | ti1=cc[t1]+cc[t4]; | ||
828 | ti2=cc[t4]-cc[t1]; | ||
829 | tr1=cc[t1-1]-cc[t4-1]; | ||
830 | tr2=cc[t1-1]+cc[t4-1]; | ||
831 | ch[t5]=tr2+tr2; | ||
832 | ch[t5+=t0]=sqrt2*(tr1-ti1); | ||
833 | ch[t5+=t0]=ti2+ti2; | ||
834 | ch[t5+=t0]=-sqrt2*(tr1+ti1); | ||
835 | |||
836 | t3+=ido; | ||
837 | t1+=t2; | ||
838 | t4+=t2; | ||
839 | } | ||
840 | } | ||
841 | |||
842 | static void dradbg(int ido,int ip,int l1,int idl1,float *cc,float *c1, | ||
843 | float *c2,float *ch,float *ch2,float *wa){ | ||
844 | static float tpi=6.283185307179586f; | ||
845 | int idij,ipph,i,j,k,l,ik,is,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10, | ||
846 | t11,t12; | ||
847 | float dc2,ai1,ai2,ar1,ar2,ds2; | ||
848 | int nbd; | ||
849 | float dcp,arg,dsp,ar1h,ar2h; | ||
850 | int ipp2; | ||
851 | |||
852 | t10=ip*ido; | ||
853 | t0=l1*ido; | ||
854 | arg=tpi/(float)ip; | ||
855 | dcp=cos(arg); | ||
856 | dsp=sin(arg); | ||
857 | nbd=(ido-1)>>1; | ||
858 | ipp2=ip; | ||
859 | ipph=(ip+1)>>1; | ||
860 | if(ido<l1)goto L103; | ||
861 | |||
862 | t1=0; | ||
863 | t2=0; | ||
864 | for(k=0;k<l1;k++){ | ||
865 | t3=t1; | ||
866 | t4=t2; | ||
867 | for(i=0;i<ido;i++){ | ||
868 | ch[t3]=cc[t4]; | ||
869 | t3++; | ||
870 | t4++; | ||
871 | } | ||
872 | t1+=ido; | ||
873 | t2+=t10; | ||
874 | } | ||
875 | goto L106; | ||
876 | |||
877 | L103: | ||
878 | t1=0; | ||
879 | for(i=0;i<ido;i++){ | ||
880 | t2=t1; | ||
881 | t3=t1; | ||
882 | for(k=0;k<l1;k++){ | ||
883 | ch[t2]=cc[t3]; | ||
884 | t2+=ido; | ||
885 | t3+=t10; | ||
886 | } | ||
887 | t1++; | ||
888 | } | ||
889 | |||
890 | L106: | ||
891 | t1=0; | ||
892 | t2=ipp2*t0; | ||
893 | t7=(t5=ido<<1); | ||
894 | for(j=1;j<ipph;j++){ | ||
895 | t1+=t0; | ||
896 | t2-=t0; | ||
897 | t3=t1; | ||
898 | t4=t2; | ||
899 | t6=t5; | ||
900 | for(k=0;k<l1;k++){ | ||
901 | ch[t3]=cc[t6-1]+cc[t6-1]; | ||
902 | ch[t4]=cc[t6]+cc[t6]; | ||
903 | t3+=ido; | ||
904 | t4+=ido; | ||
905 | t6+=t10; | ||
906 | } | ||
907 | t5+=t7; | ||
908 | } | ||
909 | |||
910 | if (ido == 1)goto L116; | ||
911 | if(nbd<l1)goto L112; | ||
912 | |||
913 | t1=0; | ||
914 | t2=ipp2*t0; | ||
915 | t7=0; | ||
916 | for(j=1;j<ipph;j++){ | ||
917 | t1+=t0; | ||
918 | t2-=t0; | ||
919 | t3=t1; | ||
920 | t4=t2; | ||
921 | |||
922 | t7+=(ido<<1); | ||
923 | t8=t7; | ||
924 | for(k=0;k<l1;k++){ | ||
925 | t5=t3; | ||
926 | t6=t4; | ||
927 | t9=t8; | ||
928 | t11=t8; | ||
929 | for(i=2;i<ido;i+=2){ | ||
930 | t5+=2; | ||
931 | t6+=2; | ||
932 | t9+=2; | ||
933 | t11-=2; | ||
934 | ch[t5-1]=cc[t9-1]+cc[t11-1]; | ||
935 | ch[t6-1]=cc[t9-1]-cc[t11-1]; | ||
936 | ch[t5]=cc[t9]-cc[t11]; | ||
937 | ch[t6]=cc[t9]+cc[t11]; | ||
938 | } | ||
939 | t3+=ido; | ||
940 | t4+=ido; | ||
941 | t8+=t10; | ||
942 | } | ||
943 | } | ||
944 | goto L116; | ||
945 | |||
946 | L112: | ||
947 | t1=0; | ||
948 | t2=ipp2*t0; | ||
949 | t7=0; | ||
950 | for(j=1;j<ipph;j++){ | ||
951 | t1+=t0; | ||
952 | t2-=t0; | ||
953 | t3=t1; | ||
954 | t4=t2; | ||
955 | t7+=(ido<<1); | ||
956 | t8=t7; | ||
957 | t9=t7; | ||
958 | for(i=2;i<ido;i+=2){ | ||
959 | t3+=2; | ||
960 | t4+=2; | ||
961 | t8+=2; | ||
962 | t9-=2; | ||
963 | t5=t3; | ||
964 | t6=t4; | ||
965 | t11=t8; | ||
966 | t12=t9; | ||
967 | for(k=0;k<l1;k++){ | ||
968 | ch[t5-1]=cc[t11-1]+cc[t12-1]; | ||
969 | ch[t6-1]=cc[t11-1]-cc[t12-1]; | ||
970 | ch[t5]=cc[t11]-cc[t12]; | ||
971 | ch[t6]=cc[t11]+cc[t12]; | ||
972 | t5+=ido; | ||
973 | t6+=ido; | ||
974 | t11+=t10; | ||
975 | t12+=t10; | ||
976 | } | ||
977 | } | ||
978 | } | ||
979 | |||
980 | L116: | ||
981 | ar1=1.f; | ||
982 | ai1=0.f; | ||
983 | t1=0; | ||
984 | t9=(t2=ipp2*idl1); | ||
985 | t3=(ip-1)*idl1; | ||
986 | for(l=1;l<ipph;l++){ | ||
987 | t1+=idl1; | ||
988 | t2-=idl1; | ||
989 | |||
990 | ar1h=dcp*ar1-dsp*ai1; | ||
991 | ai1=dcp*ai1+dsp*ar1; | ||
992 | ar1=ar1h; | ||
993 | t4=t1; | ||
994 | t5=t2; | ||
995 | t6=0; | ||
996 | t7=idl1; | ||
997 | t8=t3; | ||
998 | for(ik=0;ik<idl1;ik++){ | ||
999 | c2[t4++]=ch2[t6++]+ar1*ch2[t7++]; | ||
1000 | c2[t5++]=ai1*ch2[t8++]; | ||
1001 | } | ||
1002 | dc2=ar1; | ||
1003 | ds2=ai1; | ||
1004 | ar2=ar1; | ||
1005 | ai2=ai1; | ||
1006 | |||
1007 | t6=idl1; | ||
1008 | t7=t9-idl1; | ||
1009 | for(j=2;j<ipph;j++){ | ||
1010 | t6+=idl1; | ||
1011 | t7-=idl1; | ||
1012 | ar2h=dc2*ar2-ds2*ai2; | ||
1013 | ai2=dc2*ai2+ds2*ar2; | ||
1014 | ar2=ar2h; | ||
1015 | t4=t1; | ||
1016 | t5=t2; | ||
1017 | t11=t6; | ||
1018 | t12=t7; | ||
1019 | for(ik=0;ik<idl1;ik++){ | ||
1020 | c2[t4++]+=ar2*ch2[t11++]; | ||
1021 | c2[t5++]+=ai2*ch2[t12++]; | ||
1022 | } | ||
1023 | } | ||
1024 | } | ||
1025 | |||
1026 | t1=0; | ||
1027 | for(j=1;j<ipph;j++){ | ||
1028 | t1+=idl1; | ||
1029 | t2=t1; | ||
1030 | for(ik=0;ik<idl1;ik++)ch2[ik]+=ch2[t2++]; | ||
1031 | } | ||
1032 | |||
1033 | t1=0; | ||
1034 | t2=ipp2*t0; | ||
1035 | for(j=1;j<ipph;j++){ | ||
1036 | t1+=t0; | ||
1037 | t2-=t0; | ||
1038 | t3=t1; | ||
1039 | t4=t2; | ||
1040 | for(k=0;k<l1;k++){ | ||
1041 | ch[t3]=c1[t3]-c1[t4]; | ||
1042 | ch[t4]=c1[t3]+c1[t4]; | ||
1043 | t3+=ido; | ||
1044 | t4+=ido; | ||
1045 | } | ||
1046 | } | ||
1047 | |||
1048 | if(ido==1)goto L132; | ||
1049 | if(nbd<l1)goto L128; | ||
1050 | |||
1051 | t1=0; | ||
1052 | t2=ipp2*t0; | ||
1053 | for(j=1;j<ipph;j++){ | ||
1054 | t1+=t0; | ||
1055 | t2-=t0; | ||
1056 | t3=t1; | ||
1057 | t4=t2; | ||
1058 | for(k=0;k<l1;k++){ | ||
1059 | t5=t3; | ||
1060 | t6=t4; | ||
1061 | for(i=2;i<ido;i+=2){ | ||
1062 | t5+=2; | ||
1063 | t6+=2; | ||
1064 | ch[t5-1]=c1[t5-1]-c1[t6]; | ||
1065 | ch[t6-1]=c1[t5-1]+c1[t6]; | ||
1066 | ch[t5]=c1[t5]+c1[t6-1]; | ||
1067 | ch[t6]=c1[t5]-c1[t6-1]; | ||
1068 | } | ||
1069 | t3+=ido; | ||
1070 | t4+=ido; | ||
1071 | } | ||
1072 | } | ||
1073 | goto L132; | ||
1074 | |||
1075 | L128: | ||
1076 | t1=0; | ||
1077 | t2=ipp2*t0; | ||
1078 | for(j=1;j<ipph;j++){ | ||
1079 | t1+=t0; | ||
1080 | t2-=t0; | ||
1081 | t3=t1; | ||
1082 | t4=t2; | ||
1083 | for(i=2;i<ido;i+=2){ | ||
1084 | t3+=2; | ||
1085 | t4+=2; | ||
1086 | t5=t3; | ||
1087 | t6=t4; | ||
1088 | for(k=0;k<l1;k++){ | ||
1089 | ch[t5-1]=c1[t5-1]-c1[t6]; | ||
1090 | ch[t6-1]=c1[t5-1]+c1[t6]; | ||
1091 | ch[t5]=c1[t5]+c1[t6-1]; | ||
1092 | ch[t6]=c1[t5]-c1[t6-1]; | ||
1093 | t5+=ido; | ||
1094 | t6+=ido; | ||
1095 | } | ||
1096 | } | ||
1097 | } | ||
1098 | |||
1099 | L132: | ||
1100 | if(ido==1)return; | ||
1101 | |||
1102 | for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik]; | ||
1103 | |||
1104 | t1=0; | ||
1105 | for(j=1;j<ip;j++){ | ||
1106 | t2=(t1+=t0); | ||
1107 | for(k=0;k<l1;k++){ | ||
1108 | c1[t2]=ch[t2]; | ||
1109 | t2+=ido; | ||
1110 | } | ||
1111 | } | ||
1112 | |||
1113 | if(nbd>l1)goto L139; | ||
1114 | |||
1115 | is= -ido-1; | ||
1116 | t1=0; | ||
1117 | for(j=1;j<ip;j++){ | ||
1118 | is+=ido; | ||
1119 | t1+=t0; | ||
1120 | idij=is; | ||
1121 | t2=t1; | ||
1122 | for(i=2;i<ido;i+=2){ | ||
1123 | t2+=2; | ||
1124 | idij+=2; | ||
1125 | t3=t2; | ||
1126 | for(k=0;k<l1;k++){ | ||
1127 | c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3]; | ||
1128 | c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1]; | ||
1129 | t3+=ido; | ||
1130 | } | ||
1131 | } | ||
1132 | } | ||
1133 | return; | ||
1134 | |||
1135 | L139: | ||
1136 | is= -ido-1; | ||
1137 | t1=0; | ||
1138 | for(j=1;j<ip;j++){ | ||
1139 | is+=ido; | ||
1140 | t1+=t0; | ||
1141 | t2=t1; | ||
1142 | for(k=0;k<l1;k++){ | ||
1143 | idij=is; | ||
1144 | t3=t2; | ||
1145 | for(i=2;i<ido;i+=2){ | ||
1146 | idij+=2; | ||
1147 | t3+=2; | ||
1148 | c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3]; | ||
1149 | c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1]; | ||
1150 | } | ||
1151 | t2+=ido; | ||
1152 | } | ||
1153 | } | ||
1154 | } | ||
1155 | |||
1156 | static void drftb1(int n, float *c, float *ch, float *wa, int *ifac){ | ||
1157 | int i,k1,l1,l2; | ||
1158 | int na; | ||
1159 | int nf,ip,iw,ix2,ix3,ido,idl1; | ||
1160 | |||
1161 | nf=ifac[1]; | ||
1162 | na=0; | ||
1163 | l1=1; | ||
1164 | iw=1; | ||
1165 | |||
1166 | for(k1=0;k1<nf;k1++){ | ||
1167 | ip=ifac[k1 + 2]; | ||
1168 | l2=ip*l1; | ||
1169 | ido=n/l2; | ||
1170 | idl1=ido*l1; | ||
1171 | if(ip!=4)goto L103; | ||
1172 | ix2=iw+ido; | ||
1173 | ix3=ix2+ido; | ||
1174 | |||
1175 | if(na!=0) | ||
1176 | dradb4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1); | ||
1177 | else | ||
1178 | dradb4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1); | ||
1179 | na=1-na; | ||
1180 | goto L115; | ||
1181 | |||
1182 | L103: | ||
1183 | if(ip!=2)goto L106; | ||
1184 | |||
1185 | if(na!=0) | ||
1186 | dradb2(ido,l1,ch,c,wa+iw-1); | ||
1187 | else | ||
1188 | dradb2(ido,l1,c,ch,wa+iw-1); | ||
1189 | na=1-na; | ||
1190 | goto L115; | ||
1191 | |||
1192 | L106: | ||
1193 | if(ip!=3)goto L109; | ||
1194 | |||
1195 | ix2=iw+ido; | ||
1196 | if(na!=0) | ||
1197 | dradb3(ido,l1,ch,c,wa+iw-1,wa+ix2-1); | ||
1198 | else | ||
1199 | dradb3(ido,l1,c,ch,wa+iw-1,wa+ix2-1); | ||
1200 | na=1-na; | ||
1201 | goto L115; | ||
1202 | |||
1203 | L109: | ||
1204 | /* The radix five case can be translated later..... */ | ||
1205 | /* if(ip!=5)goto L112; | ||
1206 | |||
1207 | ix2=iw+ido; | ||
1208 | ix3=ix2+ido; | ||
1209 | ix4=ix3+ido; | ||
1210 | if(na!=0) | ||
1211 | dradb5(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1); | ||
1212 | else | ||
1213 | dradb5(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1); | ||
1214 | na=1-na; | ||
1215 | goto L115; | ||
1216 | |||
1217 | L112:*/ | ||
1218 | if(na!=0) | ||
1219 | dradbg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1); | ||
1220 | else | ||
1221 | dradbg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1); | ||
1222 | if(ido==1)na=1-na; | ||
1223 | |||
1224 | L115: | ||
1225 | l1=l2; | ||
1226 | iw+=(ip-1)*ido; | ||
1227 | } | ||
1228 | |||
1229 | if(na==0)return; | ||
1230 | |||
1231 | for(i=0;i<n;i++)c[i]=ch[i]; | ||
1232 | } | ||
1233 | |||
1234 | void spx_drft_forward(struct drft_lookup *l,float *data){ | ||
1235 | if(l->n==1)return; | ||
1236 | drftf1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache); | ||
1237 | } | ||
1238 | |||
1239 | void spx_drft_backward(struct drft_lookup *l,float *data){ | ||
1240 | if (l->n==1)return; | ||
1241 | drftb1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache); | ||
1242 | } | ||
1243 | |||
1244 | void spx_drft_init(struct drft_lookup *l,int n) | ||
1245 | { | ||
1246 | l->n=n; | ||
1247 | l->trigcache=(float*)speex_alloc(3*n*sizeof(*l->trigcache)); | ||
1248 | l->splitcache=(int*)speex_alloc(32*sizeof(*l->splitcache)); | ||
1249 | fdrffti(n, l->trigcache, l->splitcache); | ||
1250 | } | ||
1251 | |||
1252 | void spx_drft_clear(struct drft_lookup *l) | ||
1253 | { | ||
1254 | if(l) | ||
1255 | { | ||
1256 | if(l->trigcache) | ||
1257 | speex_free(l->trigcache); | ||
1258 | if(l->splitcache) | ||
1259 | speex_free(l->splitcache); | ||
1260 | } | ||
1261 | } | ||