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-2001 Kevin Rosenberg
12 ** $Id: fourier.cpp,v 1.4 2001/01/28 19:10:18 kevin Exp $
14 ** This program is free software; you can redistribute it and/or modify
15 ** it under the terms of the GNU General Public License (version 2) as
16 ** published by the Free Software Foundation.
18 ** This program is distributed in the hope that it will be useful,
19 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
20 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21 ** GNU General Public License for more details.
23 ** You should have received a copy of the GNU General Public License
24 ** along with this program; if not, write to the Free Software
25 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
26 ******************************************************************************/
33 Fourier::shuffleFourierToNaturalOrder (ImageFile& im)
35 ImageFileArray vReal = im.getArray();
36 ImageFileArray vImag = im.getImaginaryArray();
38 unsigned int nx = im.nx();
39 unsigned int ny = im.ny();
41 // shuffle each column
42 for (ix = 0; ix < nx; ix++) {
43 Fourier::shuffleFourierToNaturalOrder (vReal[ix], ny);
45 Fourier::shuffleFourierToNaturalOrder (vImag[ix], ny);
49 float* pRow = new float [nx];
50 for (iy = 0; iy < ny; iy++) {
51 for (ix = 0; ix < nx; ix++)
52 pRow[ix] = vReal[ix][iy];
53 Fourier::shuffleFourierToNaturalOrder (pRow, nx);
54 for (ix = 0; ix < nx; ix++)
55 vReal[ix][iy] = pRow[ix];
57 for (ix = 0; ix < nx; ix++)
58 pRow[ix] = vImag[ix][iy];
59 Fourier::shuffleFourierToNaturalOrder (pRow, nx);
60 for (ix = 0; ix < nx; ix++)
61 vImag[ix][iy] = pRow[ix];
68 Fourier::shuffleNaturalToFourierOrder (ImageFile& im)
70 ImageFileArray vReal = im.getArray();
71 ImageFileArray vImag = im.getImaginaryArray();
73 unsigned int nx = im.nx();
74 unsigned int ny = im.ny();
76 // shuffle each x column
77 for (ix = 0; ix < nx; ix++) {
78 Fourier::shuffleNaturalToFourierOrder (vReal[ix], ny);
80 Fourier::shuffleNaturalToFourierOrder (vImag[ix], ny);
84 float* pRow = new float [nx];
85 for (iy = 0; iy < ny; iy++) {
86 for (ix = 0; ix < nx; ix++)
87 pRow[ix] = vReal[ix][iy];
88 Fourier::shuffleNaturalToFourierOrder (pRow, nx);
89 for (ix = 0; ix < nx; ix++)
90 vReal[ix][iy] = pRow[ix];
92 for (ix = 0; ix < nx; ix++)
93 pRow[ix] = vImag[ix][iy];
94 Fourier::shuffleNaturalToFourierOrder (pRow, nx);
95 for (ix = 0; ix < nx; ix++)
96 vImag[ix][iy] = pRow[ix];
103 // Odd Number of Points
104 // Natural Frequency Order: -(n-1)/2...-1,0,1...(n-1)/2
105 // Fourier Frequency Order: 0, 1..(n-1)/2,-(n-1)/2...-1
106 // Even Number of Points
107 // Natural Frequency Order: -n/2...-1,0,1...((n/2)-1)
108 // Fourier Frequency Order: 0,1...((n/2)-1),-n/2...-1
111 Fourier::shuffleNaturalToFourierOrder (double* pdVector, const int n)
113 double* pdTemp = new double [n];
116 int iHalfN = (n - 1) / 2;
118 pdTemp[0] = pdVector[iHalfN];
119 for (i = 0; i < iHalfN; i++)
120 pdTemp[i + 1] = pdVector[i + 1 + iHalfN];
121 for (i = 0; i < iHalfN; i++)
122 pdTemp[i + iHalfN + 1] = pdVector[i];
125 pdTemp[0] = pdVector[iHalfN];
126 for (i = 0; i < iHalfN - 1; i++)
127 pdTemp[i + 1] = pdVector[i + iHalfN + 1];
128 for (i = 0; i < iHalfN; i++)
129 pdTemp[i + iHalfN] = pdVector[i];
132 for (i = 0; i < n; i++)
133 pdVector[i] = pdTemp[i];
138 Fourier::shuffleNaturalToFourierOrder (std::complex<double>* pdVector, const int n)
140 std::complex<double>* pdTemp = new std::complex<double> [n];
143 int iHalfN = (n - 1) / 2;
145 pdTemp[0] = pdVector[iHalfN];
146 for (i = 0; i < iHalfN; i++)
147 pdTemp[i + 1] = pdVector[i + 1 + iHalfN];
148 for (i = 0; i < iHalfN; i++)
149 pdTemp[i + iHalfN + 1] = pdVector[i];
152 pdTemp[0] = pdVector[iHalfN];
153 for (i = 0; i < iHalfN - 1; i++)
154 pdTemp[i + 1] = pdVector[i + iHalfN + 1];
155 for (i = 0; i < iHalfN; i++)
156 pdTemp[i + iHalfN] = pdVector[i];
159 for (i = 0; i < n; i++)
160 pdVector[i] = pdTemp[i];
166 Fourier::shuffleNaturalToFourierOrder (float* pdVector, const int n)
168 float* pdTemp = new float [n];
171 int iHalfN = (n - 1) / 2;
173 pdTemp[0] = pdVector[iHalfN];
174 for (i = 0; i < iHalfN; i++)
175 pdTemp[i + 1] = pdVector[i + 1 + iHalfN];
176 for (i = 0; i < iHalfN; i++)
177 pdTemp[i + iHalfN + 1] = pdVector[i];
180 pdTemp[0] = pdVector[iHalfN];
181 for (i = 0; i < iHalfN - 1; i++)
182 pdTemp[i + 1] = pdVector[i + iHalfN + 1];
183 for (i = 0; i < iHalfN; i++)
184 pdTemp[i + iHalfN] = pdVector[i];
187 for (i = 0; i < n; i++)
188 pdVector[i] = pdTemp[i];
195 Fourier::shuffleFourierToNaturalOrder (double* pdVector, const int n)
197 double* pdTemp = new double [n];
200 int iHalfN = (n - 1) / 2;
202 pdTemp[iHalfN] = pdVector[0];
203 for (i = 0; i < iHalfN; i++)
204 pdTemp[i + 1 + iHalfN] = pdVector[i + 1];
205 for (i = 0; i < iHalfN; i++)
206 pdTemp[i] = pdVector[i + iHalfN + 1];
209 pdTemp[iHalfN] = pdVector[0];
210 for (i = 0; i < iHalfN; i++)
211 pdTemp[i] = pdVector[i + iHalfN];
212 for (i = 0; i < iHalfN - 1; i++)
213 pdTemp[i + iHalfN + 1] = pdVector[i+1];
216 for (i = 0; i < n; i++)
217 pdVector[i] = pdTemp[i];
223 Fourier::shuffleFourierToNaturalOrder (std::complex<double>* pdVector, const int n)
225 std::complex<double>* pdTemp = new std::complex<double> [n];
228 int iHalfN = (n - 1) / 2;
230 pdTemp[iHalfN] = pdVector[0];
231 for (i = 0; i < iHalfN; i++)
232 pdTemp[i + 1 + iHalfN] = pdVector[i + 1];
233 for (i = 0; i < iHalfN; i++)
234 pdTemp[i] = pdVector[i + iHalfN + 1];
237 pdTemp[iHalfN] = pdVector[0];
238 for (i = 0; i < iHalfN; i++)
239 pdTemp[i] = pdVector[i + iHalfN];
240 for (i = 0; i < iHalfN - 1; i++)
241 pdTemp[i + iHalfN + 1] = pdVector[i+1];
244 for (i = 0; i < n; i++)
245 pdVector[i] = pdTemp[i];
253 Fourier::shuffleFourierToNaturalOrder (float* pVector, const int n)
255 float* pTemp = new float [n];
258 int iHalfN = (n - 1) / 2;
260 pTemp[iHalfN] = pVector[0];
261 for (i = 0; i < iHalfN; i++)
262 pTemp[i + 1 + iHalfN] = pVector[i + 1];
263 for (i = 0; i < iHalfN; i++)
264 pTemp[i] = pVector[i + iHalfN + 1];
267 pTemp[iHalfN] = pVector[0];
268 for (i = 0; i < iHalfN; i++)
269 pTemp[i] = pVector[i + iHalfN];
270 for (i = 0; i < iHalfN - 1; i++)
271 pTemp[i + iHalfN + 1] = pVector[i+1];
274 for (i = 0; i < n; i++)
275 pVector[i] = pTemp[i];