| 39 | // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // |
| 40 | |
| 41 | void fft::transform |
| 42 | ( |
| 43 | complexField& field, |
| 44 | const labelList& nn, |
| 45 | transformDirection isign |
| 46 | ) |
| 47 | { |
| 48 | forAll(nn, idim) |
| 49 | { |
| 50 | // Check for power of two |
| 51 | unsigned int dimCount = nn[idim]; |
| 52 | if (!dimCount || (dimCount & (dimCount - 1))) |
| 53 | { |
| 54 | FatalErrorInFunction |
| 55 | << "number of elements in direction " << idim |
| 56 | << " is not a power of 2" << endl |
| 57 | << " Number of elements in each direction = " << nn |
| 58 | << abort(FatalError); |
| 59 | } |
| 60 | } |
| 61 | |
| 62 | const label ndim = nn.size(); |
| 63 | |
| 64 | label i1, i2, i3, i2rev, i3rev, ip1, ip2, ip3, ifp1, ifp2; |
| 65 | label ibit, k1, k2, n, nprev, nrem, idim; |
| 66 | scalar tempi, tempr; |
| 67 | scalar theta, wi, wpi, wpr, wr, wtemp; |
| 68 | scalar* data = reinterpret_cast<scalar*>(field.begin()) - 1; |
| 69 | |
| 70 | |
| 71 | // if inverse transform : renumber before transform |
| 72 | |
| 73 | if (isign == REVERSE_TRANSFORM) |
| 74 | { |
| 75 | fftRenumber(field, nn); |
| 76 | } |
| 77 | |
| 78 | |
| 79 | label ntot = 1; |
| 80 | forAll(nn, idim) |
| 81 | { |
| 82 | ntot *= nn[idim]; |
| 83 | } |
| 84 | |
| 85 | |
| 86 | nprev = 1; |
| 87 | |
| 88 | for (idim=ndim; idim>=1; idim--) |
| 89 | { |
| 90 | n = nn[idim-1]; |
| 91 | nrem = ntot/(n*nprev); |
| 92 | ip1 = nprev << 1; |
| 93 | ip2 = ip1*n; |
| 94 | ip3 = ip2*nrem; |
| 95 | i2rev = 1; |
| 96 | |
| 97 | for (i2=1; i2<=ip2; i2+=ip1) |
| 98 | { |