diff options
Diffstat (limited to 'apps/plugins/pdbox/PDa/src/d_mayer_fft.c')
-rw-r--r-- | apps/plugins/pdbox/PDa/src/d_mayer_fft.c | 418 |
1 files changed, 0 insertions, 418 deletions
diff --git a/apps/plugins/pdbox/PDa/src/d_mayer_fft.c b/apps/plugins/pdbox/PDa/src/d_mayer_fft.c index 7c29c6e7c8..95fb303e91 100644 --- a/apps/plugins/pdbox/PDa/src/d_mayer_fft.c +++ b/apps/plugins/pdbox/PDa/src/d_mayer_fft.c | |||
@@ -417,422 +417,4 @@ void mayer_realifft(int n, REAL *real) | |||
417 | } | 417 | } |
418 | mayer_fht(real,n); | 418 | mayer_fht(real,n); |
419 | } | 419 | } |
420 | /* | ||
421 | ** FFT and FHT routines | ||
422 | ** Copyright 1988, 1993; Ron Mayer | ||
423 | ** | ||
424 | ** mayer_fht(fz,n); | ||
425 | ** Does a hartley transform of "n" points in the array "fz". | ||
426 | ** mayer_fft(n,real,imag) | ||
427 | ** Does a fourier transform of "n" points of the "real" and | ||
428 | ** "imag" arrays. | ||
429 | ** mayer_ifft(n,real,imag) | ||
430 | ** Does an inverse fourier transform of "n" points of the "real" | ||
431 | ** and "imag" arrays. | ||
432 | ** mayer_realfft(n,real) | ||
433 | ** Does a real-valued fourier transform of "n" points of the | ||
434 | ** "real" array. The real part of the transform ends | ||
435 | ** up in the first half of the array and the imaginary part of the | ||
436 | ** transform ends up in the second half of the array. | ||
437 | ** mayer_realifft(n,real) | ||
438 | ** The inverse of the realfft() routine above. | ||
439 | ** | ||
440 | ** | ||
441 | ** NOTE: This routine uses at least 2 patented algorithms, and may be | ||
442 | ** under the restrictions of a bunch of different organizations. | ||
443 | ** Although I wrote it completely myself, it is kind of a derivative | ||
444 | ** of a routine I once authored and released under the GPL, so it | ||
445 | ** may fall under the free software foundation's restrictions; | ||
446 | ** it was worked on as a Stanford Univ project, so they claim | ||
447 | ** some rights to it; it was further optimized at work here, so | ||
448 | ** I think this company claims parts of it. The patents are | ||
449 | ** held by R. Bracewell (the FHT algorithm) and O. Buneman (the | ||
450 | ** trig generator), both at Stanford Univ. | ||
451 | ** If it were up to me, I'd say go do whatever you want with it; | ||
452 | ** but it would be polite to give credit to the following people | ||
453 | ** if you use this anywhere: | ||
454 | ** Euler - probable inventor of the fourier transform. | ||
455 | ** Gauss - probable inventor of the FFT. | ||
456 | ** Hartley - probable inventor of the hartley transform. | ||
457 | ** Buneman - for a really cool trig generator | ||
458 | ** Mayer(me) - for authoring this particular version and | ||
459 | ** including all the optimizations in one package. | ||
460 | ** Thanks, | ||
461 | ** Ron Mayer; mayer@acuson.com | ||
462 | ** | ||
463 | */ | ||
464 | |||
465 | /* This is a slightly modified version of Mayer's contribution; write | ||
466 | * msp@ucsd.edu for the original code. Kudos to Mayer for a fine piece | ||
467 | * of work. -msp | ||
468 | */ | ||
469 | |||
470 | #ifdef MSW | ||
471 | #pragma warning( disable : 4305 ) /* uncast const double to float */ | ||
472 | #pragma warning( disable : 4244 ) /* uncast double to float */ | ||
473 | #pragma warning( disable : 4101 ) /* unused local variables */ | ||
474 | #endif | ||
475 | |||
476 | #define REAL float | ||
477 | #define GOOD_TRIG | ||
478 | |||
479 | #ifdef GOOD_TRIG | ||
480 | #else | ||
481 | #define FAST_TRIG | ||
482 | #endif | ||
483 | |||
484 | #if defined(GOOD_TRIG) | ||
485 | #define FHT_SWAP(a,b,t) {(t)=(a);(a)=(b);(b)=(t);} | ||
486 | #define TRIG_VARS \ | ||
487 | int t_lam=0; | ||
488 | #define TRIG_INIT(k,c,s) \ | ||
489 | { \ | ||
490 | int i; \ | ||
491 | for (i=2 ; i<=k ; i++) \ | ||
492 | {coswrk[i]=costab[i];sinwrk[i]=sintab[i];} \ | ||
493 | t_lam = 0; \ | ||
494 | c = 1; \ | ||
495 | s = 0; \ | ||
496 | } | ||
497 | #define TRIG_NEXT(k,c,s) \ | ||
498 | { \ | ||
499 | int i,j; \ | ||
500 | (t_lam)++; \ | ||
501 | for (i=0 ; !((1<<i)&t_lam) ; i++); \ | ||
502 | i = k-i; \ | ||
503 | s = sinwrk[i]; \ | ||
504 | c = coswrk[i]; \ | ||
505 | if (i>1) \ | ||
506 | { \ | ||
507 | for (j=k-i+2 ; (1<<j)&t_lam ; j++); \ | ||
508 | j = k - j; \ | ||
509 | sinwrk[i] = halsec[i] * (sinwrk[i-1] + sinwrk[j]); \ | ||
510 | coswrk[i] = halsec[i] * (coswrk[i-1] + coswrk[j]); \ | ||
511 | } \ | ||
512 | } | ||
513 | #define TRIG_RESET(k,c,s) | ||
514 | #endif | ||
515 | |||
516 | #if defined(FAST_TRIG) | ||
517 | #define TRIG_VARS \ | ||
518 | REAL t_c,t_s; | ||
519 | #define TRIG_INIT(k,c,s) \ | ||
520 | { \ | ||
521 | t_c = costab[k]; \ | ||
522 | t_s = sintab[k]; \ | ||
523 | c = 1; \ | ||
524 | s = 0; \ | ||
525 | } | ||
526 | #define TRIG_NEXT(k,c,s) \ | ||
527 | { \ | ||
528 | REAL t = c; \ | ||
529 | c = t*t_c - s*t_s; \ | ||
530 | s = t*t_s + s*t_c; \ | ||
531 | } | ||
532 | #define TRIG_RESET(k,c,s) | ||
533 | #endif | ||
534 | |||
535 | static REAL halsec[20]= | ||
536 | { | ||
537 | 0, | ||
538 | 0, | ||
539 | .54119610014619698439972320536638942006107206337801, | ||
540 | .50979557910415916894193980398784391368261849190893, | ||
541 | .50241928618815570551167011928012092247859337193963, | ||
542 | .50060299823519630134550410676638239611758632599591, | ||
543 | .50015063602065098821477101271097658495974913010340, | ||
544 | .50003765191554772296778139077905492847503165398345, | ||
545 | .50000941253588775676512870469186533538523133757983, | ||
546 | .50000235310628608051401267171204408939326297376426, | ||
547 | .50000058827484117879868526730916804925780637276181, | ||
548 | .50000014706860214875463798283871198206179118093251, | ||
549 | .50000003676714377807315864400643020315103490883972, | ||
550 | .50000000919178552207366560348853455333939112569380, | ||
551 | .50000000229794635411562887767906868558991922348920, | ||
552 | .50000000057448658687873302235147272458812263401372 | ||
553 | }; | ||
554 | static REAL costab[20]= | ||
555 | { | ||
556 | .00000000000000000000000000000000000000000000000000, | ||
557 | .70710678118654752440084436210484903928483593768847, | ||
558 | .92387953251128675612818318939678828682241662586364, | ||
559 | .98078528040323044912618223613423903697393373089333, | ||
560 | .99518472667219688624483695310947992157547486872985, | ||
561 | .99879545620517239271477160475910069444320361470461, | ||
562 | .99969881869620422011576564966617219685006108125772, | ||
563 | .99992470183914454092164649119638322435060646880221, | ||
564 | .99998117528260114265699043772856771617391725094433, | ||
565 | .99999529380957617151158012570011989955298763362218, | ||
566 | .99999882345170190992902571017152601904826792288976, | ||
567 | .99999970586288221916022821773876567711626389934930, | ||
568 | .99999992646571785114473148070738785694820115568892, | ||
569 | .99999998161642929380834691540290971450507605124278, | ||
570 | .99999999540410731289097193313960614895889430318945, | ||
571 | .99999999885102682756267330779455410840053741619428 | ||
572 | }; | ||
573 | static REAL sintab[20]= | ||
574 | { | ||
575 | 1.0000000000000000000000000000000000000000000000000, | ||
576 | .70710678118654752440084436210484903928483593768846, | ||
577 | .38268343236508977172845998403039886676134456248561, | ||
578 | .19509032201612826784828486847702224092769161775195, | ||
579 | .09801714032956060199419556388864184586113667316749, | ||
580 | .04906767432741801425495497694268265831474536302574, | ||
581 | .02454122852291228803173452945928292506546611923944, | ||
582 | .01227153828571992607940826195100321214037231959176, | ||
583 | .00613588464915447535964023459037258091705788631738, | ||
584 | .00306795676296597627014536549091984251894461021344, | ||
585 | .00153398018628476561230369715026407907995486457522, | ||
586 | .00076699031874270452693856835794857664314091945205, | ||
587 | .00038349518757139558907246168118138126339502603495, | ||
588 | .00019174759731070330743990956198900093346887403385, | ||
589 | .00009587379909597734587051721097647635118706561284, | ||
590 | .00004793689960306688454900399049465887274686668768 | ||
591 | }; | ||
592 | static REAL coswrk[20]= | ||
593 | { | ||
594 | .00000000000000000000000000000000000000000000000000, | ||
595 | .70710678118654752440084436210484903928483593768847, | ||
596 | .92387953251128675612818318939678828682241662586364, | ||
597 | .98078528040323044912618223613423903697393373089333, | ||
598 | .99518472667219688624483695310947992157547486872985, | ||
599 | .99879545620517239271477160475910069444320361470461, | ||
600 | .99969881869620422011576564966617219685006108125772, | ||
601 | .99992470183914454092164649119638322435060646880221, | ||
602 | .99998117528260114265699043772856771617391725094433, | ||
603 | .99999529380957617151158012570011989955298763362218, | ||
604 | .99999882345170190992902571017152601904826792288976, | ||
605 | .99999970586288221916022821773876567711626389934930, | ||
606 | .99999992646571785114473148070738785694820115568892, | ||
607 | .99999998161642929380834691540290971450507605124278, | ||
608 | .99999999540410731289097193313960614895889430318945, | ||
609 | .99999999885102682756267330779455410840053741619428 | ||
610 | }; | ||
611 | static REAL sinwrk[20]= | ||
612 | { | ||
613 | 1.0000000000000000000000000000000000000000000000000, | ||
614 | .70710678118654752440084436210484903928483593768846, | ||
615 | .38268343236508977172845998403039886676134456248561, | ||
616 | .19509032201612826784828486847702224092769161775195, | ||
617 | .09801714032956060199419556388864184586113667316749, | ||
618 | .04906767432741801425495497694268265831474536302574, | ||
619 | .02454122852291228803173452945928292506546611923944, | ||
620 | .01227153828571992607940826195100321214037231959176, | ||
621 | .00613588464915447535964023459037258091705788631738, | ||
622 | .00306795676296597627014536549091984251894461021344, | ||
623 | .00153398018628476561230369715026407907995486457522, | ||
624 | .00076699031874270452693856835794857664314091945205, | ||
625 | .00038349518757139558907246168118138126339502603495, | ||
626 | .00019174759731070330743990956198900093346887403385, | ||
627 | .00009587379909597734587051721097647635118706561284, | ||
628 | .00004793689960306688454900399049465887274686668768 | ||
629 | }; | ||
630 | |||
631 | |||
632 | #define SQRT2_2 0.70710678118654752440084436210484 | ||
633 | #define SQRT2 2*0.70710678118654752440084436210484 | ||
634 | |||
635 | void mayer_fht(REAL *fz, int n) | ||
636 | { | ||
637 | /* REAL a,b; | ||
638 | REAL c1,s1,s2,c2,s3,c3,s4,c4; | ||
639 | REAL f0,g0,f1,g1,f2,g2,f3,g3; */ | ||
640 | int k,k1,k2,k3,k4,kx; | ||
641 | REAL *fi,*fn,*gi; | ||
642 | TRIG_VARS; | ||
643 | 420 | ||
644 | for (k1=1,k2=0;k1<n;k1++) | ||
645 | { | ||
646 | REAL aa; | ||
647 | for (k=n>>1; (!((k2^=k)&k)); k>>=1); | ||
648 | if (k1>k2) | ||
649 | { | ||
650 | aa=fz[k1];fz[k1]=fz[k2];fz[k2]=aa; | ||
651 | } | ||
652 | } | ||
653 | for ( k=0 ; (1<<k)<n ; k++ ); | ||
654 | k &= 1; | ||
655 | if (k==0) | ||
656 | { | ||
657 | for (fi=fz,fn=fz+n;fi<fn;fi+=4) | ||
658 | { | ||
659 | REAL f0,f1,f2,f3; | ||
660 | f1 = fi[0 ]-fi[1 ]; | ||
661 | f0 = fi[0 ]+fi[1 ]; | ||
662 | f3 = fi[2 ]-fi[3 ]; | ||
663 | f2 = fi[2 ]+fi[3 ]; | ||
664 | fi[2 ] = (f0-f2); | ||
665 | fi[0 ] = (f0+f2); | ||
666 | fi[3 ] = (f1-f3); | ||
667 | fi[1 ] = (f1+f3); | ||
668 | } | ||
669 | } | ||
670 | else | ||
671 | { | ||
672 | for (fi=fz,fn=fz+n,gi=fi+1;fi<fn;fi+=8,gi+=8) | ||
673 | { | ||
674 | REAL bs1,bc1,bs2,bc2,bs3,bc3,bs4,bc4, | ||
675 | bg0,bf0,bf1,bg1,bf2,bg2,bf3,bg3; | ||
676 | bc1 = fi[0 ] - gi[0 ]; | ||
677 | bs1 = fi[0 ] + gi[0 ]; | ||
678 | bc2 = fi[2 ] - gi[2 ]; | ||
679 | bs2 = fi[2 ] + gi[2 ]; | ||
680 | bc3 = fi[4 ] - gi[4 ]; | ||
681 | bs3 = fi[4 ] + gi[4 ]; | ||
682 | bc4 = fi[6 ] - gi[6 ]; | ||
683 | bs4 = fi[6 ] + gi[6 ]; | ||
684 | bf1 = (bs1 - bs2); | ||
685 | bf0 = (bs1 + bs2); | ||
686 | bg1 = (bc1 - bc2); | ||
687 | bg0 = (bc1 + bc2); | ||
688 | bf3 = (bs3 - bs4); | ||
689 | bf2 = (bs3 + bs4); | ||
690 | bg3 = SQRT2*bc4; | ||
691 | bg2 = SQRT2*bc3; | ||
692 | fi[4 ] = bf0 - bf2; | ||
693 | fi[0 ] = bf0 + bf2; | ||
694 | fi[6 ] = bf1 - bf3; | ||
695 | fi[2 ] = bf1 + bf3; | ||
696 | gi[4 ] = bg0 - bg2; | ||
697 | gi[0 ] = bg0 + bg2; | ||
698 | gi[6 ] = bg1 - bg3; | ||
699 | gi[2 ] = bg1 + bg3; | ||
700 | } | ||
701 | } | ||
702 | if (n<16) return; | ||
703 | |||
704 | do | ||
705 | { | ||
706 | REAL s1,c1; | ||
707 | int ii; | ||
708 | k += 2; | ||
709 | k1 = 1 << k; | ||
710 | k2 = k1 << 1; | ||
711 | k4 = k2 << 1; | ||
712 | k3 = k2 + k1; | ||
713 | kx = k1 >> 1; | ||
714 | fi = fz; | ||
715 | gi = fi + kx; | ||
716 | fn = fz + n; | ||
717 | do | ||
718 | { | ||
719 | REAL g0,f0,f1,g1,f2,g2,f3,g3; | ||
720 | f1 = fi[0 ] - fi[k1]; | ||
721 | f0 = fi[0 ] + fi[k1]; | ||
722 | f3 = fi[k2] - fi[k3]; | ||
723 | f2 = fi[k2] + fi[k3]; | ||
724 | fi[k2] = f0 - f2; | ||
725 | fi[0 ] = f0 + f2; | ||
726 | fi[k3] = f1 - f3; | ||
727 | fi[k1] = f1 + f3; | ||
728 | g1 = gi[0 ] - gi[k1]; | ||
729 | g0 = gi[0 ] + gi[k1]; | ||
730 | g3 = SQRT2 * gi[k3]; | ||
731 | g2 = SQRT2 * gi[k2]; | ||
732 | gi[k2] = g0 - g2; | ||
733 | gi[0 ] = g0 + g2; | ||
734 | gi[k3] = g1 - g3; | ||
735 | gi[k1] = g1 + g3; | ||
736 | gi += k4; | ||
737 | fi += k4; | ||
738 | } while (fi<fn); | ||
739 | TRIG_INIT(k,c1,s1); | ||
740 | for (ii=1;ii<kx;ii++) | ||
741 | { | ||
742 | REAL c2,s2; | ||
743 | TRIG_NEXT(k,c1,s1); | ||
744 | c2 = c1*c1 - s1*s1; | ||
745 | s2 = 2*(c1*s1); | ||
746 | fn = fz + n; | ||
747 | fi = fz +ii; | ||
748 | gi = fz +k1-ii; | ||
749 | do | ||
750 | { | ||
751 | REAL a,b,g0,f0,f1,g1,f2,g2,f3,g3; | ||
752 | b = s2*fi[k1] - c2*gi[k1]; | ||
753 | a = c2*fi[k1] + s2*gi[k1]; | ||
754 | f1 = fi[0 ] - a; | ||
755 | f0 = fi[0 ] + a; | ||
756 | g1 = gi[0 ] - b; | ||
757 | g0 = gi[0 ] + b; | ||
758 | b = s2*fi[k3] - c2*gi[k3]; | ||
759 | a = c2*fi[k3] + s2*gi[k3]; | ||
760 | f3 = fi[k2] - a; | ||
761 | f2 = fi[k2] + a; | ||
762 | g3 = gi[k2] - b; | ||
763 | g2 = gi[k2] + b; | ||
764 | b = s1*f2 - c1*g3; | ||
765 | a = c1*f2 + s1*g3; | ||
766 | fi[k2] = f0 - a; | ||
767 | fi[0 ] = f0 + a; | ||
768 | gi[k3] = g1 - b; | ||
769 | gi[k1] = g1 + b; | ||
770 | b = c1*g2 - s1*f3; | ||
771 | a = s1*g2 + c1*f3; | ||
772 | gi[k2] = g0 - a; | ||
773 | gi[0 ] = g0 + a; | ||
774 | fi[k3] = f1 - b; | ||
775 | fi[k1] = f1 + b; | ||
776 | gi += k4; | ||
777 | fi += k4; | ||
778 | } while (fi<fn); | ||
779 | } | ||
780 | TRIG_RESET(k,c1,s1); | ||
781 | } while (k4<n); | ||
782 | } | ||
783 | |||
784 | void mayer_fft(int n, REAL *real, REAL *imag) | ||
785 | { | ||
786 | REAL a,b,c,d; | ||
787 | REAL q,r,s,t; | ||
788 | int i,j,k; | ||
789 | for (i=1,j=n-1,k=n/2;i<k;i++,j--) { | ||
790 | a = real[i]; b = real[j]; q=a+b; r=a-b; | ||
791 | c = imag[i]; d = imag[j]; s=c+d; t=c-d; | ||
792 | real[i] = (q+t)*.5; real[j] = (q-t)*.5; | ||
793 | imag[i] = (s-r)*.5; imag[j] = (s+r)*.5; | ||
794 | } | ||
795 | mayer_fht(real,n); | ||
796 | mayer_fht(imag,n); | ||
797 | } | ||
798 | |||
799 | void mayer_ifft(int n, REAL *real, REAL *imag) | ||
800 | { | ||
801 | REAL a,b,c,d; | ||
802 | REAL q,r,s,t; | ||
803 | int i,j,k; | ||
804 | mayer_fht(real,n); | ||
805 | mayer_fht(imag,n); | ||
806 | for (i=1,j=n-1,k=n/2;i<k;i++,j--) { | ||
807 | a = real[i]; b = real[j]; q=a+b; r=a-b; | ||
808 | c = imag[i]; d = imag[j]; s=c+d; t=c-d; | ||
809 | imag[i] = (s+r)*0.5; imag[j] = (s-r)*0.5; | ||
810 | real[i] = (q-t)*0.5; real[j] = (q+t)*0.5; | ||
811 | } | ||
812 | } | ||
813 | |||
814 | void mayer_realfft(int n, REAL *real) | ||
815 | { | ||
816 | REAL a,b,c,d; | ||
817 | int i,j,k; | ||
818 | mayer_fht(real,n); | ||
819 | for (i=1,j=n-1,k=n/2;i<k;i++,j--) { | ||
820 | a = real[i]; | ||
821 | b = real[j]; | ||
822 | real[j] = (a-b)*0.5; | ||
823 | real[i] = (a+b)*0.5; | ||
824 | } | ||
825 | } | ||
826 | |||
827 | void mayer_realifft(int n, REAL *real) | ||
828 | { | ||
829 | REAL a,b,c,d; | ||
830 | int i,j,k; | ||
831 | for (i=1,j=n-1,k=n/2;i<k;i++,j--) { | ||
832 | a = real[i]; | ||
833 | b = real[j]; | ||
834 | real[j] = (a-b); | ||
835 | real[i] = (a+b); | ||
836 | } | ||
837 | mayer_fht(real,n); | ||
838 | } | ||