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