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