| 1 | /*
 | 
|---|
| 2 | ** Please only distribute this with it's associated FHT routine.
 | 
|---|
| 3 | ** This algorithm is apparently patented(!) and the code copyrighted. 
 | 
|---|
| 4 | ** See the comment with the fht routine for more info.
 | 
|---|
| 5 | **   -Thanks,
 | 
|---|
| 6 | **   Ron Mayer
 | 
|---|
| 7 | */
 | 
|---|
| 8 | 
 | 
|---|
| 9 | #ifdef REAL
 | 
|---|
| 10 | #else
 | 
|---|
| 11 | #define REAL double
 | 
|---|
| 12 | #endif
 | 
|---|
| 13 | 
 | 
|---|
| 14 | #ifdef GOOD_TRIG
 | 
|---|
| 15 | #else
 | 
|---|
| 16 | #define FAST_TRIG
 | 
|---|
| 17 | #endif
 | 
|---|
| 18 | 
 | 
|---|
| 19 | #if defined(GOOD_TRIG)
 | 
|---|
| 20 | #define FHT_SWAP(a,b,t) {(t)=(a);(a)=(b);(b)=(t);}
 | 
|---|
| 21 | #define TRIG_VARS                                                \
 | 
|---|
| 22 |       int t_lam=0;
 | 
|---|
| 23 | #define TRIG_INIT(k,c,s)                                         \
 | 
|---|
| 24 |      {                                                           \
 | 
|---|
| 25 |       int i;                                                     \
 | 
|---|
| 26 |       for (i=2 ; i<=k ; i++)                                     \
 | 
|---|
| 27 |           {coswrk[i]=costab[i];sinwrk[i]=sintab[i];}             \
 | 
|---|
| 28 |       t_lam = 0;                                                 \
 | 
|---|
| 29 |       c = 1;                                                     \
 | 
|---|
| 30 |       s = 0;                                                     \
 | 
|---|
| 31 |      }
 | 
|---|
| 32 | #define TRIG_NEXT(k,c,s)                                         \
 | 
|---|
| 33 |      {                                                           \
 | 
|---|
| 34 |          int i,j;                                                \
 | 
|---|
| 35 |          (t_lam)++;                                              \
 | 
|---|
| 36 |          for (i=0 ; !((1<<i)&t_lam) ; i++);                      \
 | 
|---|
| 37 |          i = k-i;                                                \
 | 
|---|
| 38 |          s = sinwrk[i];                                          \
 | 
|---|
| 39 |          c = coswrk[i];                                          \
 | 
|---|
| 40 |          if (i>1)                                                \
 | 
|---|
| 41 |             {                                                    \
 | 
|---|
| 42 |              for (j=k-i+2 ; (1<<j)&t_lam ; j++);                 \
 | 
|---|
| 43 |              j         = k - j;                                  \
 | 
|---|
| 44 |              sinwrk[i] = halsec[i] * (sinwrk[i-1] + sinwrk[j]);  \
 | 
|---|
| 45 |              coswrk[i] = halsec[i] * (coswrk[i-1] + coswrk[j]);  \
 | 
|---|
| 46 |             }                                                    \
 | 
|---|
| 47 |      }
 | 
|---|
| 48 | #define TRIG_RESET(k,c,s)
 | 
|---|
| 49 | #endif
 | 
|---|
| 50 | 
 | 
|---|
| 51 | #if defined(FAST_TRIG)
 | 
|---|
| 52 | #define TRIG_VARS                                        \
 | 
|---|
| 53 |       REAL t_c,t_s;
 | 
|---|
| 54 | #define TRIG_INIT(k,c,s)                                 \
 | 
|---|
| 55 |     {                                                    \
 | 
|---|
| 56 |      t_c  = costab[k];                                   \
 | 
|---|
| 57 |      t_s  = sintab[k];                                   \
 | 
|---|
| 58 |      c    = 1;                                           \
 | 
|---|
| 59 |      s    = 0;                                           \
 | 
|---|
| 60 |     }
 | 
|---|
| 61 | #define TRIG_NEXT(k,c,s)                                 \
 | 
|---|
| 62 |     {                                                    \
 | 
|---|
| 63 |      REAL t = c;                                         \
 | 
|---|
| 64 |      c   = t*t_c - s*t_s;                                \
 | 
|---|
| 65 |      s   = t*t_s + s*t_c;                                \
 | 
|---|
| 66 |     }
 | 
|---|
| 67 | #define TRIG_RESET(k,c,s)
 | 
|---|
| 68 | #endif
 | 
|---|
| 69 | 
 | 
|---|
| 70 | static REAL halsec[20]=
 | 
|---|
| 71 |     {
 | 
|---|
| 72 |      0,
 | 
|---|
| 73 |      0,
 | 
|---|
| 74 |      .54119610014619698439972320536638942006107206337801,
 | 
|---|
| 75 |      .50979557910415916894193980398784391368261849190893,
 | 
|---|
| 76 |      .50241928618815570551167011928012092247859337193963,
 | 
|---|
| 77 |      .50060299823519630134550410676638239611758632599591,
 | 
|---|
| 78 |      .50015063602065098821477101271097658495974913010340,
 | 
|---|
| 79 |      .50003765191554772296778139077905492847503165398345,
 | 
|---|
| 80 |      .50000941253588775676512870469186533538523133757983,
 | 
|---|
| 81 |      .50000235310628608051401267171204408939326297376426,
 | 
|---|
| 82 |      .50000058827484117879868526730916804925780637276181,
 | 
|---|
| 83 |      .50000014706860214875463798283871198206179118093251,
 | 
|---|
| 84 |      .50000003676714377807315864400643020315103490883972,
 | 
|---|
| 85 |      .50000000919178552207366560348853455333939112569380,
 | 
|---|
| 86 |      .50000000229794635411562887767906868558991922348920,
 | 
|---|
| 87 |      .50000000057448658687873302235147272458812263401372
 | 
|---|
| 88 |     };
 | 
|---|
| 89 | static REAL costab[20]=
 | 
|---|
| 90 |     {
 | 
|---|
| 91 |      .00000000000000000000000000000000000000000000000000,
 | 
|---|
| 92 |      .70710678118654752440084436210484903928483593768847,
 | 
|---|
| 93 |      .92387953251128675612818318939678828682241662586364,
 | 
|---|
| 94 |      .98078528040323044912618223613423903697393373089333,
 | 
|---|
| 95 |      .99518472667219688624483695310947992157547486872985,
 | 
|---|
| 96 |      .99879545620517239271477160475910069444320361470461,
 | 
|---|
| 97 |      .99969881869620422011576564966617219685006108125772,
 | 
|---|
| 98 |      .99992470183914454092164649119638322435060646880221,
 | 
|---|
| 99 |      .99998117528260114265699043772856771617391725094433,
 | 
|---|
| 100 |      .99999529380957617151158012570011989955298763362218,
 | 
|---|
| 101 |      .99999882345170190992902571017152601904826792288976,
 | 
|---|
| 102 |      .99999970586288221916022821773876567711626389934930,
 | 
|---|
| 103 |      .99999992646571785114473148070738785694820115568892,
 | 
|---|
| 104 |      .99999998161642929380834691540290971450507605124278,
 | 
|---|
| 105 |      .99999999540410731289097193313960614895889430318945,
 | 
|---|
| 106 |      .99999999885102682756267330779455410840053741619428
 | 
|---|
| 107 |     };
 | 
|---|
| 108 | static REAL sintab[20]=
 | 
|---|
| 109 |     {
 | 
|---|
| 110 |      1.0000000000000000000000000000000000000000000000000,
 | 
|---|
| 111 |      .70710678118654752440084436210484903928483593768846,
 | 
|---|
| 112 |      .38268343236508977172845998403039886676134456248561,
 | 
|---|
| 113 |      .19509032201612826784828486847702224092769161775195,
 | 
|---|
| 114 |      .09801714032956060199419556388864184586113667316749,
 | 
|---|
| 115 |      .04906767432741801425495497694268265831474536302574,
 | 
|---|
| 116 |      .02454122852291228803173452945928292506546611923944,
 | 
|---|
| 117 |      .01227153828571992607940826195100321214037231959176,
 | 
|---|
| 118 |      .00613588464915447535964023459037258091705788631738,
 | 
|---|
| 119 |      .00306795676296597627014536549091984251894461021344,
 | 
|---|
| 120 |      .00153398018628476561230369715026407907995486457522,
 | 
|---|
| 121 |      .00076699031874270452693856835794857664314091945205,
 | 
|---|
| 122 |      .00038349518757139558907246168118138126339502603495,
 | 
|---|
| 123 |      .00019174759731070330743990956198900093346887403385,
 | 
|---|
| 124 |      .00009587379909597734587051721097647635118706561284,
 | 
|---|
| 125 |      .00004793689960306688454900399049465887274686668768
 | 
|---|
| 126 |     };
 | 
|---|
| 127 | static REAL coswrk[20]=
 | 
|---|
| 128 |     {
 | 
|---|
| 129 |      .00000000000000000000000000000000000000000000000000,
 | 
|---|
| 130 |      .70710678118654752440084436210484903928483593768847,
 | 
|---|
| 131 |      .92387953251128675612818318939678828682241662586364,
 | 
|---|
| 132 |      .98078528040323044912618223613423903697393373089333,
 | 
|---|
| 133 |      .99518472667219688624483695310947992157547486872985,
 | 
|---|
| 134 |      .99879545620517239271477160475910069444320361470461,
 | 
|---|
| 135 |      .99969881869620422011576564966617219685006108125772,
 | 
|---|
| 136 |      .99992470183914454092164649119638322435060646880221,
 | 
|---|
| 137 |      .99998117528260114265699043772856771617391725094433,
 | 
|---|
| 138 |      .99999529380957617151158012570011989955298763362218,
 | 
|---|
| 139 |      .99999882345170190992902571017152601904826792288976,
 | 
|---|
| 140 |      .99999970586288221916022821773876567711626389934930,
 | 
|---|
| 141 |      .99999992646571785114473148070738785694820115568892,
 | 
|---|
| 142 |      .99999998161642929380834691540290971450507605124278,
 | 
|---|
| 143 |      .99999999540410731289097193313960614895889430318945,
 | 
|---|
| 144 |      .99999999885102682756267330779455410840053741619428
 | 
|---|
| 145 |     };
 | 
|---|
| 146 | static REAL sinwrk[20]=
 | 
|---|
| 147 |     {
 | 
|---|
| 148 |      1.0000000000000000000000000000000000000000000000000,
 | 
|---|
| 149 |      .70710678118654752440084436210484903928483593768846,
 | 
|---|
| 150 |      .38268343236508977172845998403039886676134456248561,
 | 
|---|
| 151 |      .19509032201612826784828486847702224092769161775195,
 | 
|---|
| 152 |      .09801714032956060199419556388864184586113667316749,
 | 
|---|
| 153 |      .04906767432741801425495497694268265831474536302574,
 | 
|---|
| 154 |      .02454122852291228803173452945928292506546611923944,
 | 
|---|
| 155 |      .01227153828571992607940826195100321214037231959176,
 | 
|---|
| 156 |      .00613588464915447535964023459037258091705788631738,
 | 
|---|
| 157 |      .00306795676296597627014536549091984251894461021344,
 | 
|---|
| 158 |      .00153398018628476561230369715026407907995486457522,
 | 
|---|
| 159 |      .00076699031874270452693856835794857664314091945205,
 | 
|---|
| 160 |      .00038349518757139558907246168118138126339502603495,
 | 
|---|
| 161 |      .00019174759731070330743990956198900093346887403385,
 | 
|---|
| 162 |      .00009587379909597734587051721097647635118706561284,
 | 
|---|
| 163 |      .00004793689960306688454900399049465887274686668768
 | 
|---|
| 164 |     };
 | 
|---|