Performs an optimized 1D FHT of an array or part of an array. @param x Input array; will be overwritten by the output in the range given by base and maxN. @param base First index from where data of the input array should be read. @param inverse True for inverse transform. @param maxN
(float[] x, int base, boolean inverse, int maxN)
| 278 | * Note that all amplitudes in the output 'x' are multiplied by maxN. |
| 279 | */ |
| 280 | public void dfht3(float[] x, int base, boolean inverse, int maxN) { |
| 281 | int i, stage, gpNum, gpIndex, gpSize, numGps, Nlog2; |
| 282 | int bfNum, numBfs; |
| 283 | int Ad0, Ad1, Ad2, Ad3, Ad4, CSAd; |
| 284 | float rt1, rt2, rt3, rt4; |
| 285 | |
| 286 | if (S==null) initializeTables(maxN); |
| 287 | Nlog2 = log2(maxN); |
| 288 | BitRevRArr(x, base, Nlog2, maxN); //bitReverse the input array |
| 289 | gpSize = 2; //first & second stages - do radix 4 butterflies once thru |
| 290 | numGps = maxN / 4; |
| 291 | for (gpNum=0; gpNum<numGps; gpNum++) { |
| 292 | Ad1 = gpNum * 4; |
| 293 | Ad2 = Ad1 + 1; |
| 294 | Ad3 = Ad1 + gpSize; |
| 295 | Ad4 = Ad2 + gpSize; |
| 296 | rt1 = x[base+Ad1] + x[base+Ad2]; // a + b |
| 297 | rt2 = x[base+Ad1] - x[base+Ad2]; // a - b |
| 298 | rt3 = x[base+Ad3] + x[base+Ad4]; // c + d |
| 299 | rt4 = x[base+Ad3] - x[base+Ad4]; // c - d |
| 300 | x[base+Ad1] = rt1 + rt3; // a + b + (c + d) |
| 301 | x[base+Ad2] = rt2 + rt4; // a - b + (c - d) |
| 302 | x[base+Ad3] = rt1 - rt3; // a + b - (c + d) |
| 303 | x[base+Ad4] = rt2 - rt4; // a - b - (c - d) |
| 304 | } |
| 305 | |
| 306 | if (Nlog2 > 2) { |
| 307 | // third + stages computed here |
| 308 | gpSize = 4; |
| 309 | numBfs = 2; |
| 310 | numGps = numGps / 2; |
| 311 | for (stage=2; stage<Nlog2; stage++) { |
| 312 | for (gpNum=0; gpNum<numGps; gpNum++) { |
| 313 | Ad0 = gpNum * gpSize * 2; |
| 314 | Ad1 = Ad0; // 1st butterfly is different from others - no mults needed |
| 315 | Ad2 = Ad1 + gpSize; |
| 316 | Ad3 = Ad1 + gpSize / 2; |
| 317 | Ad4 = Ad3 + gpSize; |
| 318 | rt1 = x[base+Ad1]; |
| 319 | x[base+Ad1] = x[base+Ad1] + x[base+Ad2]; |
| 320 | x[base+Ad2] = rt1 - x[base+Ad2]; |
| 321 | rt1 = x[base+Ad3]; |
| 322 | x[base+Ad3] = x[base+Ad3] + x[base+Ad4]; |
| 323 | x[base+Ad4] = rt1 - x[base+Ad4]; |
| 324 | for (bfNum=1; bfNum<numBfs; bfNum++) { |
| 325 | // subsequent BF's dealt with together |
| 326 | Ad1 = bfNum + Ad0; |
| 327 | Ad2 = Ad1 + gpSize; |
| 328 | Ad3 = gpSize - bfNum + Ad0; |
| 329 | Ad4 = Ad3 + gpSize; |
| 330 | |
| 331 | CSAd = bfNum * numGps; |
| 332 | rt1 = x[base+Ad2] * C[CSAd] + x[base+Ad4] * S[CSAd]; |
| 333 | rt2 = x[base+Ad4] * C[CSAd] - x[base+Ad2] * S[CSAd]; |
| 334 | |
| 335 | x[base+Ad2] = x[base+Ad1] - rt1; |
| 336 | x[base+Ad1] = x[base+Ad1] + rt1; |
| 337 | x[base+Ad4] = x[base+Ad3] + rt2; |
no test coverage detected