1 /*****************************************************************************
5 ** Purpose: Fourier transform functions
6 ** Programmer: Kevin Rosenberg
7 ** Date Started: Dec 2000
9 ** This is part of the CTSim program
10 ** Copyright (c) 1983-2009 Kevin Rosenberg
12 ** This program is free software; you can redistribute it and/or modify
13 ** it under the terms of the GNU General Public License (version 2) as
14 ** published by the Free Software Foundation.
16 ** This program is distributed in the hope that it will be useful,
17 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
18 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 ** GNU General Public License for more details.
21 ** You should have received a copy of the GNU General Public License
22 ** along with this program; if not, write to the Free Software
23 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
24 ******************************************************************************/
31 Fourier::shuffleFourierToNaturalOrder (ImageFile& im)
33 ImageFileArray vReal = im.getArray();
34 ImageFileArray vImag = im.getImaginaryArray();
36 unsigned int nx = im.nx();
37 unsigned int ny = im.ny();
39 // shuffle each column
40 for (ix = 0; ix < nx; ix++) {
41 Fourier::shuffleFourierToNaturalOrder (vReal[ix], ny);
43 Fourier::shuffleFourierToNaturalOrder (vImag[ix], ny);
47 float* pRow = new float [nx];
48 for (iy = 0; iy < ny; iy++) {
49 for (ix = 0; ix < nx; ix++)
50 pRow[ix] = vReal[ix][iy];
51 Fourier::shuffleFourierToNaturalOrder (pRow, nx);
52 for (ix = 0; ix < nx; ix++)
53 vReal[ix][iy] = pRow[ix];
55 for (ix = 0; ix < nx; ix++)
56 pRow[ix] = vImag[ix][iy];
57 Fourier::shuffleFourierToNaturalOrder (pRow, nx);
58 for (ix = 0; ix < nx; ix++)
59 vImag[ix][iy] = pRow[ix];
66 Fourier::shuffleNaturalToFourierOrder (ImageFile& im)
68 ImageFileArray vReal = im.getArray();
69 ImageFileArray vImag = im.getImaginaryArray();
71 unsigned int nx = im.nx();
72 unsigned int ny = im.ny();
74 // shuffle each x column
75 for (ix = 0; ix < nx; ix++) {
76 Fourier::shuffleNaturalToFourierOrder (vReal[ix], ny);
78 Fourier::shuffleNaturalToFourierOrder (vImag[ix], ny);
82 float* pRow = new float [nx];
83 for (iy = 0; iy < ny; iy++) {
84 for (ix = 0; ix < nx; ix++)
85 pRow[ix] = vReal[ix][iy];
86 Fourier::shuffleNaturalToFourierOrder (pRow, nx);
87 for (ix = 0; ix < nx; ix++)
88 vReal[ix][iy] = pRow[ix];
90 for (ix = 0; ix < nx; ix++)
91 pRow[ix] = vImag[ix][iy];
92 Fourier::shuffleNaturalToFourierOrder (pRow, nx);
93 for (ix = 0; ix < nx; ix++)
94 vImag[ix][iy] = pRow[ix];
101 void Fourier::shuffleNaturalToFourierOrder (fftw_complex* pVector, const int n)
103 fftw_complex* pTemp = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * n));
106 if (isOdd(n)) { // Odd
107 int iHalfN = (n - 1) / 2;
109 pTemp[0][0] = pVector[iHalfN][0];
110 pTemp[0][1] = pVector[iHalfN][1];
111 for (i = 0; i < iHalfN; i++) {
112 pTemp[i + 1][0] = pVector[i + 1 + iHalfN][0];
113 pTemp[i + 1][1] = pVector[i + 1 + iHalfN][1];
115 for (i = 0; i < iHalfN; i++) {
116 pTemp[i + iHalfN + 1][0] = pVector[i][0];
117 pTemp[i + iHalfN + 1][1] = pVector[i][1];
121 pTemp[0][0] = pVector[iHalfN][0];
122 pTemp[0][1] = pVector[iHalfN][1];
123 for (i = 0; i < iHalfN - 1; i++) {
124 pTemp[i + 1][0] = pVector[i + iHalfN + 1][0];
125 pTemp[i + 1][1] = pVector[i + iHalfN + 1][1];
127 for (i = 0; i < iHalfN; i++) {
128 pTemp[i + iHalfN][0] = pVector[i][0];
129 pTemp[i + iHalfN][1] = pVector[i][1];
133 for (i = 0; i < n; i++) {
134 pVector[i][0] = pTemp[i][0];
135 pVector[i][1] = pTemp[i][1];
140 void Fourier::shuffleFourierToNaturalOrder (fftw_complex* pVector, const int n)
142 fftw_complex* pTemp = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * n));
144 if (isOdd(n)) { // Odd
145 int iHalfN = (n - 1) / 2;
147 pTemp[iHalfN][0] = pVector[0][0];
148 pTemp[iHalfN][1] = pVector[0][1];
149 for (i = 0; i < iHalfN; i++) {
150 pTemp[i + 1 + iHalfN][0] = pVector[i + 1][0];
151 pTemp[i + 1 + iHalfN][1] = pVector[i + 1][1];
153 for (i = 0; i < iHalfN; i++) {
154 pTemp[i][0] = pVector[i + iHalfN + 1][0];
155 pTemp[i][1] = pVector[i + iHalfN + 1][1];
159 pTemp[iHalfN][0] = pVector[0][0];
160 pTemp[iHalfN][1] = pVector[0][1];
161 for (i = 0; i < iHalfN; i++) {
162 pTemp[i][0] = pVector[i + iHalfN][0];
163 pTemp[i][1] = pVector[i + iHalfN][1];
165 for (i = 0; i < iHalfN - 1; i++) {
166 pTemp[i + iHalfN + 1][0] = pVector[i+1][0];
167 pTemp[i + iHalfN + 1][1] = pVector[i+1][1];
171 for (i = 0; i < n; i++) {
172 pVector[i][0] = pTemp[i][0];
173 pVector[i][1] = pTemp[i][1];