vrshoot: e7ca128b8713 libs/vorbis/smallft.c

vrshoot

view libs/vorbis/smallft.c @ 1:e7ca128b8713

looks nice :)

author	John Tsiombikas <nuclear@member.fsf.org>
date	Sun, 02 Feb 2014 00:35:22 +0200
parents
children

line source

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 * *

9 * by the Xiph.Org Foundation http://www.xiph.org/ *

10 * *

11 ********************************************************************

13 function: *unnormalized* fft transform

14 last mod: $Id: smallft.c 16227 2009-07-08 06:58:46Z xiphmont $

16 ********************************************************************/

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 */

31 #include <stdlib.h>

32 #include <string.h>

33 #include <math.h>

34 #include "smallft.h"

35 #include "os.h"

36 #include "misc.h"

38 static void drfti1(int n, float *wa, int *ifac){

39 static int ntryh[4] = { 4,2,3,5 };

40 static float tpi = 6.28318530717958648f;

41 float arg,argh,argld,fi;

42 int ntry=0,i,j=-1;

43 int k1, l1, l2, ib;

44 int ld, ii, ip, is, nq, nr;

45 int ido, ipm, nfm1;

46 int nl=n;

47 int nf=0;

49 L101:

50 j++;

51 if (j < 4)

52 ntry=ntryh[j];

53 else

54 ntry+=2;

56 L104:

57 nq=nl/ntry;

58 nr=nl-ntry*nq;

59 if (nr!=0) goto L101;

61 nf++;

62 ifac[nf+1]=ntry;

63 nl=nq;

64 if(ntry!=2)goto L107;

65 if(nf==1)goto L107;

67 for (i=1;i<nf;i++){

68 ib=nf-i+1;

69 ifac[ib+1]=ifac[ib];

70 }

71 ifac[2] = 2;

73 L107:

74 if(nl!=1)goto L104;

75 ifac[0]=n;

76 ifac[1]=nf;

77 argh=tpi/n;

78 is=0;

79 nfm1=nf-1;

80 l1=1;

82 if(nfm1==0)return;

84 for (k1=0;k1<nfm1;k1++){

85 ip=ifac[k1+2];

86 ld=0;

87 l2=l1*ip;

88 ido=n/l2;

89 ipm=ip-1;

91 for (j=0;j<ipm;j++){

92 ld+=l1;

93 i=is;

94 argld=(float)ld*argh;

95 fi=0.f;

96 for (ii=2;ii<ido;ii+=2){

97 fi+=1.f;

98 arg=fi*argld;

99 wa[i++]=cos(arg);

100 wa[i++]=sin(arg);

101 }

102 is+=ido;

103 }

104 l1=l2;

105 }

106 }

107

108 static void fdrffti(int n, float *wsave, int *ifac){

109

110 if (n == 1) return;

111 drfti1(n, wsave+n, ifac);

112 }

113

114 static void dradf2(int ido,int l1,float *cc,float *ch,float *wa1){

115 int i,k;

116 float ti2,tr2;

117 int t0,t1,t2,t3,t4,t5,t6;

118

119 t1=0;

120 t0=(t2=l1*ido);

121 t3=ido<<1;

122 for(k=0;k<l1;k++){

123 ch[t1<<1]=cc[t1]+cc[t2];

124 ch[(t1<<1)+t3-1]=cc[t1]-cc[t2];

125 t1+=ido;

126 t2+=ido;

127 }

128

129 if(ido<2)return;

130 if(ido==2)goto L105;

131

132 t1=0;

133 t2=t0;

134 for(k=0;k<l1;k++){

135 t3=t2;

136 t4=(t1<<1)+(ido<<1);

137 t5=t1;

138 t6=t1+t1;

139 for(i=2;i<ido;i+=2){

140 t3+=2;

141 t4-=2;

142 t5+=2;

143 t6+=2;

144 tr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];

145 ti2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];

146 ch[t6]=cc[t5]+ti2;

147 ch[t4]=ti2-cc[t5];

148 ch[t6-1]=cc[t5-1]+tr2;

149 ch[t4-1]=cc[t5-1]-tr2;

150 }

151 t1+=ido;

152 t2+=ido;

153 }

154

155 if(ido%2==1)return;

156

157 L105:

158 t3=(t2=(t1=ido)-1);

159 t2+=t0;

160 for(k=0;k<l1;k++){

161 ch[t1]=-cc[t2];

162 ch[t1-1]=cc[t3];

163 t1+=ido<<1;

164 t2+=ido;

165 t3+=ido;

166 }

167 }

168

169 static void dradf4(int ido,int l1,float *cc,float *ch,float *wa1,

170 float *wa2,float *wa3){

171 static float hsqt2 = .70710678118654752f;

172 int i,k,t0,t1,t2,t3,t4,t5,t6;

173 float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;

174 t0=l1*ido;

175

176 t1=t0;

177 t4=t1<<1;

178 t2=t1+(t1<<1);

179 t3=0;

180

181 for(k=0;k<l1;k++){

182 tr1=cc[t1]+cc[t2];

183 tr2=cc[t3]+cc[t4];

184

185 ch[t5=t3<<2]=tr1+tr2;

186 ch[(ido<<2)+t5-1]=tr2-tr1;

187 ch[(t5+=(ido<<1))-1]=cc[t3]-cc[t4];

188 ch[t5]=cc[t2]-cc[t1];

189

190 t1+=ido;

191 t2+=ido;

192 t3+=ido;

193 t4+=ido;

194 }

195

196 if(ido<2)return;

197 if(ido==2)goto L105;

198

199

200 t1=0;

201 for(k=0;k<l1;k++){

202 t2=t1;

203 t4=t1<<2;

204 t5=(t6=ido<<1)+t4;

205 for(i=2;i<ido;i+=2){

206 t3=(t2+=2);

207 t4+=2;

208 t5-=2;

209

210 t3+=t0;

211 cr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];

212 ci2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];

213 t3+=t0;

214 cr3=wa2[i-2]*cc[t3-1]+wa2[i-1]*cc[t3];

215 ci3=wa2[i-2]*cc[t3]-wa2[i-1]*cc[t3-1];

216 t3+=t0;

217 cr4=wa3[i-2]*cc[t3-1]+wa3[i-1]*cc[t3];

218 ci4=wa3[i-2]*cc[t3]-wa3[i-1]*cc[t3-1];

219

220 tr1=cr2+cr4;

221 tr4=cr4-cr2;

222 ti1=ci2+ci4;

223 ti4=ci2-ci4;

224

225 ti2=cc[t2]+ci3;

226 ti3=cc[t2]-ci3;

227 tr2=cc[t2-1]+cr3;

228 tr3=cc[t2-1]-cr3;

229

230 ch[t4-1]=tr1+tr2;

231 ch[t4]=ti1+ti2;

232

233 ch[t5-1]=tr3-ti4;

234 ch[t5]=tr4-ti3;

235

236 ch[t4+t6-1]=ti4+tr3;

237 ch[t4+t6]=tr4+ti3;

238

239 ch[t5+t6-1]=tr2-tr1;

240 ch[t5+t6]=ti1-ti2;

241 }

242 t1+=ido;

243 }

244 if(ido&1)return;

245

246 L105:

247

248 t2=(t1=t0+ido-1)+(t0<<1);

249 t3=ido<<2;

250 t4=ido;

251 t5=ido<<1;

252 t6=ido;

253

254 for(k=0;k<l1;k++){

255 ti1=-hsqt2*(cc[t1]+cc[t2]);

256 tr1=hsqt2*(cc[t1]-cc[t2]);

257

258 ch[t4-1]=tr1+cc[t6-1];

259 ch[t4+t5-1]=cc[t6-1]-tr1;

260

261 ch[t4]=ti1-cc[t1+t0];

262 ch[t4+t5]=ti1+cc[t1+t0];

263

264 t1+=ido;

265 t2+=ido;

266 t4+=t3;

267 t6+=ido;

268 }

269 }

270

271 static void dradfg(int ido,int ip,int l1,int idl1,float *cc,float *c1,

272 float *c2,float *ch,float *ch2,float *wa){

273

274 static float tpi=6.283185307179586f;

275 int idij,ipph,i,j,k,l,ic,ik,is;

276 int t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;

277 float dc2,ai1,ai2,ar1,ar2,ds2;

278 int nbd;

279 float dcp,arg,dsp,ar1h,ar2h;

280 int idp2,ipp2;

281

282 arg=tpi/(float)ip;

283 dcp=cos(arg);

284 dsp=sin(arg);

285 ipph=(ip+1)>>1;

286 ipp2=ip;

287 idp2=ido;

288 nbd=(ido-1)>>1;

289 t0=l1*ido;

290 t10=ip*ido;

291

292 if(ido==1)goto L119;

293 for(ik=0;ik<idl1;ik++)ch2[ik]=c2[ik];

294

295 t1=0;

296 for(j=1;j<ip;j++){

297 t1+=t0;

298 t2=t1;

299 for(k=0;k<l1;k++){

300 ch[t2]=c1[t2];

301 t2+=ido;

302 }

303 }

304

305 is=-ido;

306 t1=0;

307 if(nbd>l1){

308 for(j=1;j<ip;j++){

309 t1+=t0;

310 is+=ido;

311 t2= -ido+t1;

312 for(k=0;k<l1;k++){

313 idij=is-1;

314 t2+=ido;

315 t3=t2;

316 for(i=2;i<ido;i+=2){

317 idij+=2;

318 t3+=2;

319 ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3];

320 ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1];

321 }

322 }

323 }

324 }else{

325

326 for(j=1;j<ip;j++){

327 is+=ido;

328 idij=is-1;

329 t1+=t0;

330 t2=t1;

331 for(i=2;i<ido;i+=2){

332 idij+=2;

333 t2+=2;

334 t3=t2;

335 for(k=0;k<l1;k++){

336 ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3];

337 ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1];

338 t3+=ido;

339 }

340 }

341 }

342 }

343

344 t1=0;

345 t2=ipp2*t0;

346 if(nbd<l1){

347 for(j=1;j<ipph;j++){

348 t1+=t0;

349 t2-=t0;

350 t3=t1;

351 t4=t2;

352 for(i=2;i<ido;i+=2){

353 t3+=2;

354 t4+=2;

355 t5=t3-ido;

356 t6=t4-ido;

357 for(k=0;k<l1;k++){

358 t5+=ido;

359 t6+=ido;

360 c1[t5-1]=ch[t5-1]+ch[t6-1];

361 c1[t6-1]=ch[t5]-ch[t6];

362 c1[t5]=ch[t5]+ch[t6];

363 c1[t6]=ch[t6-1]-ch[t5-1];

364 }

365 }

366 }

367 }else{

368 for(j=1;j<ipph;j++){

369 t1+=t0;

370 t2-=t0;

371 t3=t1;

372 t4=t2;

373 for(k=0;k<l1;k++){

374 t5=t3;

375 t6=t4;

376 for(i=2;i<ido;i+=2){

377 t5+=2;

378 t6+=2;

379 c1[t5-1]=ch[t5-1]+ch[t6-1];

380 c1[t6-1]=ch[t5]-ch[t6];

381 c1[t5]=ch[t5]+ch[t6];

382 c1[t6]=ch[t6-1]-ch[t5-1];

383 }

384 t3+=ido;

385 t4+=ido;

386 }

387 }

388 }

389

390 L119:

391 for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];

392

393 t1=0;

394 t2=ipp2*idl1;

395 for(j=1;j<ipph;j++){

396 t1+=t0;

397 t2-=t0;

398 t3=t1-ido;

399 t4=t2-ido;

400 for(k=0;k<l1;k++){

401 t3+=ido;

402 t4+=ido;

403 c1[t3]=ch[t3]+ch[t4];

404 c1[t4]=ch[t4]-ch[t3];

405 }

406 }

407

408 ar1=1.f;

409 ai1=0.f;

410 t1=0;

411 t2=ipp2*idl1;

412 t3=(ip-1)*idl1;

413 for(l=1;l<ipph;l++){

414 t1+=idl1;

415 t2-=idl1;

416 ar1h=dcp*ar1-dsp*ai1;

417 ai1=dcp*ai1+dsp*ar1;

418 ar1=ar1h;

419 t4=t1;

420 t5=t2;

421 t6=t3;

422 t7=idl1;

423

424 for(ik=0;ik<idl1;ik++){

425 ch2[t4++]=c2[ik]+ar1*c2[t7++];

426 ch2[t5++]=ai1*c2[t6++];

427 }

428

429 dc2=ar1;

430 ds2=ai1;

431 ar2=ar1;

432 ai2=ai1;

433

434 t4=idl1;

435 t5=(ipp2-1)*idl1;

436 for(j=2;j<ipph;j++){

437 t4+=idl1;

438 t5-=idl1;

439

440 ar2h=dc2*ar2-ds2*ai2;

441 ai2=dc2*ai2+ds2*ar2;

442 ar2=ar2h;

443

444 t6=t1;

445 t7=t2;

446 t8=t4;

447 t9=t5;

448 for(ik=0;ik<idl1;ik++){

449 ch2[t6++]+=ar2*c2[t8++];

450 ch2[t7++]+=ai2*c2[t9++];

451 }

452 }

453 }

454

455 t1=0;

456 for(j=1;j<ipph;j++){

457 t1+=idl1;

458 t2=t1;

459 for(ik=0;ik<idl1;ik++)ch2[ik]+=c2[t2++];

460 }

461

462 if(ido<l1)goto L132;

463

464 t1=0;

465 t2=0;

466 for(k=0;k<l1;k++){

467 t3=t1;

468 t4=t2;

469 for(i=0;i<ido;i++)cc[t4++]=ch[t3++];

470 t1+=ido;

471 t2+=t10;

472 }

473

474 goto L135;

475

476 L132:

477 for(i=0;i<ido;i++){

478 t1=i;

479 t2=i;

480 for(k=0;k<l1;k++){

481 cc[t2]=ch[t1];

482 t1+=ido;

483 t2+=t10;

484 }

485 }

486

487 L135:

488 t1=0;

489 t2=ido<<1;

490 t3=0;

491 t4=ipp2*t0;

492 for(j=1;j<ipph;j++){

493

494 t1+=t2;

495 t3+=t0;

496 t4-=t0;

497

498 t5=t1;

499 t6=t3;

500 t7=t4;

501

502 for(k=0;k<l1;k++){

503 cc[t5-1]=ch[t6];

504 cc[t5]=ch[t7];

505 t5+=t10;

506 t6+=ido;

507 t7+=ido;

508 }

509 }

510

511 if(ido==1)return;

512 if(nbd<l1)goto L141;

513

514 t1=-ido;

515 t3=0;

516 t4=0;

517 t5=ipp2*t0;

518 for(j=1;j<ipph;j++){

519 t1+=t2;

520 t3+=t2;

521 t4+=t0;

522 t5-=t0;

523 t6=t1;

524 t7=t3;

525 t8=t4;

526 t9=t5;

527 for(k=0;k<l1;k++){

528 for(i=2;i<ido;i+=2){

529 ic=idp2-i;

530 cc[i+t7-1]=ch[i+t8-1]+ch[i+t9-1];

531 cc[ic+t6-1]=ch[i+t8-1]-ch[i+t9-1];

532 cc[i+t7]=ch[i+t8]+ch[i+t9];

533 cc[ic+t6]=ch[i+t9]-ch[i+t8];

534 }

535 t6+=t10;

536 t7+=t10;

537 t8+=ido;

538 t9+=ido;

539 }

540 }

541 return;

542

543 L141:

544

545 t1=-ido;

546 t3=0;

547 t4=0;

548 t5=ipp2*t0;

549 for(j=1;j<ipph;j++){

550 t1+=t2;

551 t3+=t2;

552 t4+=t0;

553 t5-=t0;

554 for(i=2;i<ido;i+=2){

555 t6=idp2+t1-i;

556 t7=i+t3;

557 t8=i+t4;

558 t9=i+t5;

559 for(k=0;k<l1;k++){

560 cc[t7-1]=ch[t8-1]+ch[t9-1];

561 cc[t6-1]=ch[t8-1]-ch[t9-1];

562 cc[t7]=ch[t8]+ch[t9];

563 cc[t6]=ch[t9]-ch[t8];

564 t6+=t10;

565 t7+=t10;

566 t8+=ido;

567 t9+=ido;

568 }

569 }

570 }

571 }

572

573 static void drftf1(int n,float *c,float *ch,float *wa,int *ifac){

574 int i,k1,l1,l2;

575 int na,kh,nf;

576 int ip,iw,ido,idl1,ix2,ix3;

577

578 nf=ifac[1];

579 na=1;

580 l2=n;

581 iw=n;

582

583 for(k1=0;k1<nf;k1++){

584 kh=nf-k1;

585 ip=ifac[kh+1];

586 l1=l2/ip;

587 ido=n/l2;

588 idl1=ido*l1;

589 iw-=(ip-1)*ido;

590 na=1-na;

591

592 if(ip!=4)goto L102;

593

594 ix2=iw+ido;

595 ix3=ix2+ido;

596 if(na!=0)

597 dradf4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);

598 else

599 dradf4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);

600 goto L110;

601

602 L102:

603 if(ip!=2)goto L104;

604 if(na!=0)goto L103;

605

606 dradf2(ido,l1,c,ch,wa+iw-1);

607 goto L110;

608

609 L103:

610 dradf2(ido,l1,ch,c,wa+iw-1);

611 goto L110;

612

613 L104:

614 if(ido==1)na=1-na;

615 if(na!=0)goto L109;

616

617 dradfg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1);

618 na=1;

619 goto L110;

620

621 L109:

622 dradfg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1);

623 na=0;

624

625 L110:

626 l2=l1;

627 }

628

629 if(na==1)return;

630

631 for(i=0;i<n;i++)c[i]=ch[i];

632 }

633

634 static void dradb2(int ido,int l1,float *cc,float *ch,float *wa1){

635 int i,k,t0,t1,t2,t3,t4,t5,t6;

636 float ti2,tr2;

637

638 t0=l1*ido;

639

640 t1=0;

641 t2=0;

642 t3=(ido<<1)-1;

643 for(k=0;k<l1;k++){

644 ch[t1]=cc[t2]+cc[t3+t2];

645 ch[t1+t0]=cc[t2]-cc[t3+t2];

646 t2=(t1+=ido)<<1;

647 }

648

649 if(ido<2)return;

650 if(ido==2)goto L105;

651

652 t1=0;

653 t2=0;

654 for(k=0;k<l1;k++){

655 t3=t1;

656 t5=(t4=t2)+(ido<<1);

657 t6=t0+t1;

658 for(i=2;i<ido;i+=2){

659 t3+=2;

660 t4+=2;

661 t5-=2;

662 t6+=2;

663 ch[t3-1]=cc[t4-1]+cc[t5-1];

664 tr2=cc[t4-1]-cc[t5-1];

665 ch[t3]=cc[t4]-cc[t5];

666 ti2=cc[t4]+cc[t5];

667 ch[t6-1]=wa1[i-2]*tr2-wa1[i-1]*ti2;

668 ch[t6]=wa1[i-2]*ti2+wa1[i-1]*tr2;

669 }

670 t2=(t1+=ido)<<1;

671 }

672

673 if(ido%2==1)return;

674

675 L105:

676 t1=ido-1;

677 t2=ido-1;

678 for(k=0;k<l1;k++){

679 ch[t1]=cc[t2]+cc[t2];

680 ch[t1+t0]=-(cc[t2+1]+cc[t2+1]);

681 t1+=ido;

682 t2+=ido<<1;

683 }

684 }

685

686 static void dradb3(int ido,int l1,float *cc,float *ch,float *wa1,

687 float *wa2){

688 static float taur = -.5f;

689 static float taui = .8660254037844386f;

690 int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;

691 float ci2,ci3,di2,di3,cr2,cr3,dr2,dr3,ti2,tr2;

692 t0=l1*ido;

693

694 t1=0;

695 t2=t0<<1;

696 t3=ido<<1;

697 t4=ido+(ido<<1);

698 t5=0;

699 for(k=0;k<l1;k++){

700 tr2=cc[t3-1]+cc[t3-1];

701 cr2=cc[t5]+(taur*tr2);

702 ch[t1]=cc[t5]+tr2;

703 ci3=taui*(cc[t3]+cc[t3]);

704 ch[t1+t0]=cr2-ci3;

705 ch[t1+t2]=cr2+ci3;

706 t1+=ido;

707 t3+=t4;

708 t5+=t4;

709 }

710

711 if(ido==1)return;

712

713 t1=0;

714 t3=ido<<1;

715 for(k=0;k<l1;k++){

716 t7=t1+(t1<<1);

717 t6=(t5=t7+t3);

718 t8=t1;

719 t10=(t9=t1+t0)+t0;

720

721 for(i=2;i<ido;i+=2){

722 t5+=2;

723 t6-=2;

724 t7+=2;

725 t8+=2;

726 t9+=2;

727 t10+=2;

728 tr2=cc[t5-1]+cc[t6-1];

729 cr2=cc[t7-1]+(taur*tr2);

730 ch[t8-1]=cc[t7-1]+tr2;

731 ti2=cc[t5]-cc[t6];

732 ci2=cc[t7]+(taur*ti2);

733 ch[t8]=cc[t7]+ti2;

734 cr3=taui*(cc[t5-1]-cc[t6-1]);

735 ci3=taui*(cc[t5]+cc[t6]);

736 dr2=cr2-ci3;

737 dr3=cr2+ci3;

738 di2=ci2+cr3;

739 di3=ci2-cr3;

740 ch[t9-1]=wa1[i-2]*dr2-wa1[i-1]*di2;

741 ch[t9]=wa1[i-2]*di2+wa1[i-1]*dr2;

742 ch[t10-1]=wa2[i-2]*dr3-wa2[i-1]*di3;

743 ch[t10]=wa2[i-2]*di3+wa2[i-1]*dr3;

744 }

745 t1+=ido;

746 }

747 }

748

749 static void dradb4(int ido,int l1,float *cc,float *ch,float *wa1,

750 float *wa2,float *wa3){

751 static float sqrt2=1.414213562373095f;

752 int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8;

753 float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;

754 t0=l1*ido;

755

756 t1=0;

757 t2=ido<<2;

758 t3=0;

759 t6=ido<<1;

760 for(k=0;k<l1;k++){

761 t4=t3+t6;

762 t5=t1;

763 tr3=cc[t4-1]+cc[t4-1];

764 tr4=cc[t4]+cc[t4];

765 tr1=cc[t3]-cc[(t4+=t6)-1];

766 tr2=cc[t3]+cc[t4-1];

767 ch[t5]=tr2+tr3;

768 ch[t5+=t0]=tr1-tr4;

769 ch[t5+=t0]=tr2-tr3;

770 ch[t5+=t0]=tr1+tr4;

771 t1+=ido;

772 t3+=t2;

773 }

774

775 if(ido<2)return;

776 if(ido==2)goto L105;

777

778 t1=0;

779 for(k=0;k<l1;k++){

780 t5=(t4=(t3=(t2=t1<<2)+t6))+t6;

781 t7=t1;

782 for(i=2;i<ido;i+=2){

783 t2+=2;

784 t3+=2;

785 t4-=2;

786 t5-=2;

787 t7+=2;

788 ti1=cc[t2]+cc[t5];

789 ti2=cc[t2]-cc[t5];

790 ti3=cc[t3]-cc[t4];

791 tr4=cc[t3]+cc[t4];

792 tr1=cc[t2-1]-cc[t5-1];

793 tr2=cc[t2-1]+cc[t5-1];

794 ti4=cc[t3-1]-cc[t4-1];

795 tr3=cc[t3-1]+cc[t4-1];

796 ch[t7-1]=tr2+tr3;

797 cr3=tr2-tr3;

798 ch[t7]=ti2+ti3;

799 ci3=ti2-ti3;

800 cr2=tr1-tr4;

801 cr4=tr1+tr4;

802 ci2=ti1+ti4;

803 ci4=ti1-ti4;

804

805 ch[(t8=t7+t0)-1]=wa1[i-2]*cr2-wa1[i-1]*ci2;

806 ch[t8]=wa1[i-2]*ci2+wa1[i-1]*cr2;

807 ch[(t8+=t0)-1]=wa2[i-2]*cr3-wa2[i-1]*ci3;

808 ch[t8]=wa2[i-2]*ci3+wa2[i-1]*cr3;

809 ch[(t8+=t0)-1]=wa3[i-2]*cr4-wa3[i-1]*ci4;

810 ch[t8]=wa3[i-2]*ci4+wa3[i-1]*cr4;

811 }

812 t1+=ido;

813 }

814

815 if(ido%2 == 1)return;

816

817 L105:

818

819 t1=ido;

820 t2=ido<<2;

821 t3=ido-1;

822 t4=ido+(ido<<1);

823 for(k=0;k<l1;k++){

824 t5=t3;

825 ti1=cc[t1]+cc[t4];

826 ti2=cc[t4]-cc[t1];

827 tr1=cc[t1-1]-cc[t4-1];

828 tr2=cc[t1-1]+cc[t4-1];

829 ch[t5]=tr2+tr2;

830 ch[t5+=t0]=sqrt2*(tr1-ti1);

831 ch[t5+=t0]=ti2+ti2;

832 ch[t5+=t0]=-sqrt2*(tr1+ti1);

833

834 t3+=ido;

835 t1+=t2;

836 t4+=t2;

837 }

838 }

839

840 static void dradbg(int ido,int ip,int l1,int idl1,float *cc,float *c1,

841 float *c2,float *ch,float *ch2,float *wa){

842 static float tpi=6.283185307179586f;

843 int idij,ipph,i,j,k,l,ik,is,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,

844 t11,t12;

845 float dc2,ai1,ai2,ar1,ar2,ds2;

846 int nbd;

847 float dcp,arg,dsp,ar1h,ar2h;

848 int ipp2;

849

850 t10=ip*ido;

851 t0=l1*ido;

852 arg=tpi/(float)ip;

853 dcp=cos(arg);

854 dsp=sin(arg);

855 nbd=(ido-1)>>1;

856 ipp2=ip;

857 ipph=(ip+1)>>1;

858 if(ido<l1)goto L103;

859

860 t1=0;

861 t2=0;

862 for(k=0;k<l1;k++){

863 t3=t1;

864 t4=t2;

865 for(i=0;i<ido;i++){

866 ch[t3]=cc[t4];

867 t3++;

868 t4++;

869 }

870 t1+=ido;

871 t2+=t10;

872 }

873 goto L106;

874

875 L103:

876 t1=0;

877 for(i=0;i<ido;i++){

878 t2=t1;

879 t3=t1;

880 for(k=0;k<l1;k++){

881 ch[t2]=cc[t3];

882 t2+=ido;

883 t3+=t10;

884 }

885 t1++;

886 }

887

888 L106:

889 t1=0;

890 t2=ipp2*t0;

891 t7=(t5=ido<<1);

892 for(j=1;j<ipph;j++){

893 t1+=t0;

894 t2-=t0;

895 t3=t1;

896 t4=t2;

897 t6=t5;

898 for(k=0;k<l1;k++){

899 ch[t3]=cc[t6-1]+cc[t6-1];

900 ch[t4]=cc[t6]+cc[t6];

901 t3+=ido;

902 t4+=ido;

903 t6+=t10;

904 }

905 t5+=t7;

906 }

907

908 if (ido == 1)goto L116;

909 if(nbd<l1)goto L112;

910

911 t1=0;

912 t2=ipp2*t0;

913 t7=0;

914 for(j=1;j<ipph;j++){

915 t1+=t0;

916 t2-=t0;

917 t3=t1;

918 t4=t2;

919

920 t7+=(ido<<1);

921 t8=t7;

922 for(k=0;k<l1;k++){

923 t5=t3;

924 t6=t4;

925 t9=t8;

926 t11=t8;

927 for(i=2;i<ido;i+=2){

928 t5+=2;

929 t6+=2;

930 t9+=2;

931 t11-=2;

932 ch[t5-1]=cc[t9-1]+cc[t11-1];

933 ch[t6-1]=cc[t9-1]-cc[t11-1];

934 ch[t5]=cc[t9]-cc[t11];

935 ch[t6]=cc[t9]+cc[t11];

936 }

937 t3+=ido;

938 t4+=ido;

939 t8+=t10;

940 }

941 }

942 goto L116;

943

944 L112:

945 t1=0;

946 t2=ipp2*t0;

947 t7=0;

948 for(j=1;j<ipph;j++){

949 t1+=t0;

950 t2-=t0;

951 t3=t1;

952 t4=t2;

953 t7+=(ido<<1);

954 t8=t7;

955 t9=t7;

956 for(i=2;i<ido;i+=2){

957 t3+=2;

958 t4+=2;

959 t8+=2;

960 t9-=2;

961 t5=t3;

962 t6=t4;

963 t11=t8;

964 t12=t9;

965 for(k=0;k<l1;k++){

966 ch[t5-1]=cc[t11-1]+cc[t12-1];

967 ch[t6-1]=cc[t11-1]-cc[t12-1];

968 ch[t5]=cc[t11]-cc[t12];

969 ch[t6]=cc[t11]+cc[t12];

970 t5+=ido;

971 t6+=ido;

972 t11+=t10;

973 t12+=t10;

974 }

975 }

976 }

977

978 L116:

979 ar1=1.f;

980 ai1=0.f;

981 t1=0;

982 t9=(t2=ipp2*idl1);

983 t3=(ip-1)*idl1;

984 for(l=1;l<ipph;l++){

985 t1+=idl1;

986 t2-=idl1;

987

988 ar1h=dcp*ar1-dsp*ai1;

989 ai1=dcp*ai1+dsp*ar1;

990 ar1=ar1h;

991 t4=t1;

992 t5=t2;

993 t6=0;

994 t7=idl1;

995 t8=t3;

996 for(ik=0;ik<idl1;ik++){

997 c2[t4++]=ch2[t6++]+ar1*ch2[t7++];

998 c2[t5++]=ai1*ch2[t8++];

999 }

1000 dc2=ar1;

1001 ds2=ai1;

1002 ar2=ar1;

1003 ai2=ai1;

1004

1005 t6=idl1;

1006 t7=t9-idl1;

1007 for(j=2;j<ipph;j++){

1008 t6+=idl1;

1009 t7-=idl1;

1010 ar2h=dc2*ar2-ds2*ai2;

1011 ai2=dc2*ai2+ds2*ar2;

1012 ar2=ar2h;

1013 t4=t1;

1014 t5=t2;

1015 t11=t6;

1016 t12=t7;

1017 for(ik=0;ik<idl1;ik++){

1018 c2[t4++]+=ar2*ch2[t11++];

1019 c2[t5++]+=ai2*ch2[t12++];

1020 }

1021 }

1022 }

1023

1024 t1=0;

1025 for(j=1;j<ipph;j++){

1026 t1+=idl1;

1027 t2=t1;

1028 for(ik=0;ik<idl1;ik++)ch2[ik]+=ch2[t2++];

1029 }

1030

1031 t1=0;

1032 t2=ipp2*t0;

1033 for(j=1;j<ipph;j++){

1034 t1+=t0;

1035 t2-=t0;

1036 t3=t1;

1037 t4=t2;

1038 for(k=0;k<l1;k++){

1039 ch[t3]=c1[t3]-c1[t4];

1040 ch[t4]=c1[t3]+c1[t4];

1041 t3+=ido;

1042 t4+=ido;

1043 }

1044 }

1045

1046 if(ido==1)goto L132;

1047 if(nbd<l1)goto L128;

1048

1049 t1=0;

1050 t2=ipp2*t0;

1051 for(j=1;j<ipph;j++){

1052 t1+=t0;

1053 t2-=t0;

1054 t3=t1;

1055 t4=t2;

1056 for(k=0;k<l1;k++){

1057 t5=t3;

1058 t6=t4;

1059 for(i=2;i<ido;i+=2){

1060 t5+=2;

1061 t6+=2;

1062 ch[t5-1]=c1[t5-1]-c1[t6];

1063 ch[t6-1]=c1[t5-1]+c1[t6];

1064 ch[t5]=c1[t5]+c1[t6-1];

1065 ch[t6]=c1[t5]-c1[t6-1];

1066 }

1067 t3+=ido;

1068 t4+=ido;

1069 }

1070 }

1071 goto L132;

1072

1073 L128:

1074 t1=0;

1075 t2=ipp2*t0;

1076 for(j=1;j<ipph;j++){

1077 t1+=t0;

1078 t2-=t0;

1079 t3=t1;

1080 t4=t2;

1081 for(i=2;i<ido;i+=2){

1082 t3+=2;

1083 t4+=2;

1084 t5=t3;

1085 t6=t4;

1086 for(k=0;k<l1;k++){

1087 ch[t5-1]=c1[t5-1]-c1[t6];

1088 ch[t6-1]=c1[t5-1]+c1[t6];

1089 ch[t5]=c1[t5]+c1[t6-1];

1090 ch[t6]=c1[t5]-c1[t6-1];

1091 t5+=ido;

1092 t6+=ido;

1093 }

1094 }

1095 }

1096

1097 L132:

1098 if(ido==1)return;

1099

1100 for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];

1101

1102 t1=0;

1103 for(j=1;j<ip;j++){

1104 t2=(t1+=t0);

1105 for(k=0;k<l1;k++){

1106 c1[t2]=ch[t2];

1107 t2+=ido;

1108 }

1109 }

1110

1111 if(nbd>l1)goto L139;

1112

1113 is= -ido-1;

1114 t1=0;

1115 for(j=1;j<ip;j++){

1116 is+=ido;

1117 t1+=t0;

1118 idij=is;

1119 t2=t1;

1120 for(i=2;i<ido;i+=2){

1121 t2+=2;

1122 idij+=2;

1123 t3=t2;

1124 for(k=0;k<l1;k++){

1125 c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3];

1126 c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1];

1127 t3+=ido;

1128 }

1129 }

1130 }

1131 return;

1132

1133 L139:

1134 is= -ido-1;

1135 t1=0;

1136 for(j=1;j<ip;j++){

1137 is+=ido;

1138 t1+=t0;

1139 t2=t1;

1140 for(k=0;k<l1;k++){

1141 idij=is;

1142 t3=t2;

1143 for(i=2;i<ido;i+=2){

1144 idij+=2;

1145 t3+=2;

1146 c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3];

1147 c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1];

1148 }

1149 t2+=ido;

1150 }

1151 }

1152 }

1153

1154 static void drftb1(int n, float *c, float *ch, float *wa, int *ifac){

1155 int i,k1,l1,l2;

1156 int na;

1157 int nf,ip,iw,ix2,ix3,ido,idl1;

1158

1159 nf=ifac[1];

1160 na=0;

1161 l1=1;

1162 iw=1;

1163

1164 for(k1=0;k1<nf;k1++){

1165 ip=ifac[k1 + 2];

1166 l2=ip*l1;

1167 ido=n/l2;

1168 idl1=ido*l1;

1169 if(ip!=4)goto L103;

1170 ix2=iw+ido;

1171 ix3=ix2+ido;

1172

1173 if(na!=0)

1174 dradb4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);

1175 else

1176 dradb4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);

1177 na=1-na;

1178 goto L115;

1179

1180 L103:

1181 if(ip!=2)goto L106;

1182

1183 if(na!=0)

1184 dradb2(ido,l1,ch,c,wa+iw-1);

1185 else

1186 dradb2(ido,l1,c,ch,wa+iw-1);

1187 na=1-na;

1188 goto L115;

1189

1190 L106:

1191 if(ip!=3)goto L109;

1192

1193 ix2=iw+ido;

1194 if(na!=0)

1195 dradb3(ido,l1,ch,c,wa+iw-1,wa+ix2-1);

1196 else

1197 dradb3(ido,l1,c,ch,wa+iw-1,wa+ix2-1);

1198 na=1-na;

1199 goto L115;

1200

1201 L109:

1202 /* The radix five case can be translated later..... */

1203 /* if(ip!=5)goto L112;

1204

1205 ix2=iw+ido;

1206 ix3=ix2+ido;

1207 ix4=ix3+ido;

1208 if(na!=0)

1209 dradb5(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);

1210 else

1211 dradb5(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);

1212 na=1-na;

1213 goto L115;

1214

1215 L112:*/

1216 if(na!=0)

1217 dradbg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1);

1218 else

1219 dradbg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1);

1220 if(ido==1)na=1-na;

1221

1222 L115:

1223 l1=l2;

1224 iw+=(ip-1)*ido;

1225 }

1226

1227 if(na==0)return;

1228

1229 for(i=0;i<n;i++)c[i]=ch[i];

1230 }

1231

1232 void drft_forward(drft_lookup *l,float *data){

1233 if(l->n==1)return;

1234 drftf1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);

1235 }

1236

1237 void drft_backward(drft_lookup *l,float *data){

1238 if (l->n==1)return;

1239 drftb1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);

1240 }

1241

1242 void drft_init(drft_lookup *l,int n){

1243 l->n=n;

1244 l->trigcache=_ogg_calloc(3*n,sizeof(*l->trigcache));

1245 l->splitcache=_ogg_calloc(32,sizeof(*l->splitcache));

1246 fdrffti(n, l->trigcache, l->splitcache);

1247 }

1248

1249 void drft_clear(drft_lookup *l){

1250 if(l){

1251 if(l->trigcache)_ogg_free(l->trigcache);

1252 if(l->splitcache)_ogg_free(l->splitcache);

1253 memset(l,0,sizeof(*l));

1254 }

1255 }