** Date Started: Dec 2000
**
** This is part of the CTSim program
-** Copyright (C) 1983-2001 Kevin Rosenberg
-**
-** $Id: fourier.cpp,v 1.2 2001/01/02 16:02:13 kevin Exp $
+** Copyright (c) 1983-2009 Kevin Rosenberg
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License (version 2) as
if (im.isComplex()) {
for (ix = 0; ix < nx; ix++)
pRow[ix] = vImag[ix][iy];
- Fourier::shuffleFourierToNaturalOrder (pRow, nx);;
+ Fourier::shuffleFourierToNaturalOrder (pRow, nx);
for (ix = 0; ix < nx; ix++)
vImag[ix][iy] = pRow[ix];
}
}
delete pRow;
}
-
+
void
Fourier::shuffleNaturalToFourierOrder (ImageFile& im)
{
delete [] pRow;
}
-
-// Odd Number of Points
-// Natural Frequency Order: -(n-1)/2...-1,0,1...(n-1)/2
-// Fourier Frequency Order: 0, 1..(n-1)/2,-(n-1)/2...-1
-// Even Number of Points
-// Natural Frequency Order: -n/2...-1,0,1...((n/2)-1)
-// Fourier Frequency Order: 0,1...((n/2)-1),-n/2...-1
-
-void
-Fourier::shuffleNaturalToFourierOrder (double* pdVector, const int n)
+#ifdef HAVE_FFTW
+void Fourier::shuffleNaturalToFourierOrder (fftw_complex* pVector, const int n)
{
- double* pdTemp = new double [n];
+ fftw_complex* pTemp = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * n));
int i;
- if (n % 2) { // Odd
- int iHalfN = (n - 1) / 2;
-
- pdTemp[0] = pdVector[iHalfN];
- for (i = 0; i < iHalfN; i++)
- pdTemp[i + 1] = pdVector[i + 1 + iHalfN];
- for (i = 0; i < iHalfN; i++)
- pdTemp[i + iHalfN + 1] = pdVector[i];
- } else { // Even
- int iHalfN = n / 2;
- pdTemp[0] = pdVector[iHalfN];
- for (i = 0; i < iHalfN - 1; i++)
- pdTemp[i + 1] = pdVector[i + iHalfN + 1];
- for (i = 0; i < iHalfN; i++)
- pdTemp[i + iHalfN] = pdVector[i];
- }
-
- for (i = 0; i < n; i++)
- pdVector[i] = pdTemp[i];
- delete pdTemp;
-}
-void
-Fourier::shuffleNaturalToFourierOrder (std::complex<double>* pdVector, const int n)
-{
- std::complex<double>* pdTemp = new std::complex<double> [n];
- int i;
- if (n % 2) { // Odd
+ if (isOdd(n)) { // Odd
int iHalfN = (n - 1) / 2;
-
- pdTemp[0] = pdVector[iHalfN];
- for (i = 0; i < iHalfN; i++)
- pdTemp[i + 1] = pdVector[i + 1 + iHalfN];
- for (i = 0; i < iHalfN; i++)
- pdTemp[i + iHalfN + 1] = pdVector[i];
- } else { // Even
- int iHalfN = n / 2;
- pdTemp[0] = pdVector[iHalfN];
- for (i = 0; i < iHalfN - 1; i++)
- pdTemp[i + 1] = pdVector[i + iHalfN + 1];
- for (i = 0; i < iHalfN; i++)
- pdTemp[i + iHalfN] = pdVector[i];
- }
-
- for (i = 0; i < n; i++)
- pdVector[i] = pdTemp[i];
- delete [] pdTemp;
-}
-
-void
-Fourier::shuffleNaturalToFourierOrder (float* pdVector, const int n)
-{
- float* pdTemp = new float [n];
- int i;
- if (n % 2) { // Odd
- int iHalfN = (n - 1) / 2;
-
- pdTemp[0] = pdVector[iHalfN];
- for (i = 0; i < iHalfN; i++)
- pdTemp[i + 1] = pdVector[i + 1 + iHalfN];
- for (i = 0; i < iHalfN; i++)
- pdTemp[i + iHalfN + 1] = pdVector[i];
+ pTemp[0][0] = pVector[iHalfN][0];
+ pTemp[0][1] = pVector[iHalfN][1];
+ for (i = 0; i < iHalfN; i++) {
+ pTemp[i + 1][0] = pVector[i + 1 + iHalfN][0];
+ pTemp[i + 1][1] = pVector[i + 1 + iHalfN][1];
+ }
+ for (i = 0; i < iHalfN; i++) {
+ pTemp[i + iHalfN + 1][0] = pVector[i][0];
+ pTemp[i + iHalfN + 1][1] = pVector[i][1];
+ }
} else { // Even
int iHalfN = n / 2;
- pdTemp[0] = pdVector[iHalfN];
- for (i = 0; i < iHalfN - 1; i++)
- pdTemp[i + 1] = pdVector[i + iHalfN + 1];
- for (i = 0; i < iHalfN; i++)
- pdTemp[i + iHalfN] = pdVector[i];
+ pTemp[0][0] = pVector[iHalfN][0];
+ pTemp[0][1] = pVector[iHalfN][1];
+ for (i = 0; i < iHalfN - 1; i++) {
+ pTemp[i + 1][0] = pVector[i + iHalfN + 1][0];
+ pTemp[i + 1][1] = pVector[i + iHalfN + 1][1];
+ }
+ for (i = 0; i < iHalfN; i++) {
+ pTemp[i + iHalfN][0] = pVector[i][0];
+ pTemp[i + iHalfN][1] = pVector[i][1];
+ }
}
-
- for (i = 0; i < n; i++)
- pdVector[i] = pdTemp[i];
- delete pdTemp;
-}
-
-
-void
-Fourier::shuffleFourierToNaturalOrder (double* pdVector, const int n)
-{
- double* pdTemp = new double [n];
- int i;
- if (n % 2) { // Odd
- int iHalfN = (n - 1) / 2;
-
- pdTemp[iHalfN] = pdVector[0];
- for (i = 0; i < iHalfN; i++)
- pdTemp[i + 1 + iHalfN] = pdVector[i + 1];
- for (i = 0; i < iHalfN; i++)
- pdTemp[i] = pdVector[i + iHalfN + 1];
- } else { // Even
- int iHalfN = n / 2;
- pdTemp[iHalfN] = pdVector[0];
- for (i = 0; i < iHalfN; i++)
- pdTemp[i] = pdVector[i + iHalfN];
- for (i = 0; i < iHalfN - 1; i++)
- pdTemp[i + iHalfN + 1] = pdVector[i+1];
+ for (i = 0; i < n; i++) {
+ pVector[i][0] = pTemp[i][0];
+ pVector[i][1] = pTemp[i][1];
}
-
- for (i = 0; i < n; i++)
- pdVector[i] = pdTemp[i];
- delete pdTemp;
+ fftw_free(pTemp);
}
-
-void
-Fourier::shuffleFourierToNaturalOrder (std::complex<double>* pdVector, const int n)
+void Fourier::shuffleFourierToNaturalOrder (fftw_complex* pVector, const int n)
{
- std::complex<double>* pdTemp = new std::complex<double> [n];
+ fftw_complex* pTemp = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * n));
int i;
- if (n % 2) { // Odd
+ if (isOdd(n)) { // Odd
int iHalfN = (n - 1) / 2;
-
- pdTemp[iHalfN] = pdVector[0];
- for (i = 0; i < iHalfN; i++)
- pdTemp[i + 1 + iHalfN] = pdVector[i + 1];
- for (i = 0; i < iHalfN; i++)
- pdTemp[i] = pdVector[i + iHalfN + 1];
+
+ pTemp[iHalfN][0] = pVector[0][0];
+ pTemp[iHalfN][1] = pVector[0][1];
+ for (i = 0; i < iHalfN; i++) {
+ pTemp[i + 1 + iHalfN][0] = pVector[i + 1][0];
+ pTemp[i + 1 + iHalfN][1] = pVector[i + 1][1];
+ }
+ for (i = 0; i < iHalfN; i++) {
+ pTemp[i][0] = pVector[i + iHalfN + 1][0];
+ pTemp[i][1] = pVector[i + iHalfN + 1][1];
+ }
} else { // Even
int iHalfN = n / 2;
- pdTemp[iHalfN] = pdVector[0];
- for (i = 0; i < iHalfN; i++)
- pdTemp[i] = pdVector[i + iHalfN];
- for (i = 0; i < iHalfN - 1; i++)
- pdTemp[i + iHalfN + 1] = pdVector[i+1];
+ pTemp[iHalfN][0] = pVector[0][0];
+ pTemp[iHalfN][1] = pVector[0][1];
+ for (i = 0; i < iHalfN; i++) {
+ pTemp[i][0] = pVector[i + iHalfN][0];
+ pTemp[i][1] = pVector[i + iHalfN][1];
+ }
+ for (i = 0; i < iHalfN - 1; i++) {
+ pTemp[i + iHalfN + 1][0] = pVector[i+1][0];
+ pTemp[i + iHalfN + 1][1] = pVector[i+1][1];
+ }
}
-
- for (i = 0; i < n; i++)
- pdVector[i] = pdTemp[i];
- delete [] pdTemp;
-}
-
-
-
-void
-Fourier::shuffleFourierToNaturalOrder (float* pVector, const int n)
-{
- float* pTemp = new float [n];
- int i;
- if (n % 2) { // Odd
- int iHalfN = (n - 1) / 2;
-
- pTemp[iHalfN] = pVector[0];
- for (i = 0; i < iHalfN; i++)
- pTemp[i + 1 + iHalfN] = pVector[i + 1];
- for (i = 0; i < iHalfN; i++)
- pTemp[i] = pVector[i + iHalfN + 1];
- } else { // Even
- int iHalfN = n / 2;
- pTemp[iHalfN] = pVector[0];
- for (i = 0; i < iHalfN; i++)
- pTemp[i] = pVector[i + iHalfN];
- for (i = 0; i < iHalfN - 1; i++)
- pTemp[i + iHalfN + 1] = pVector[i+1];
+ for (i = 0; i < n; i++) {
+ pVector[i][0] = pTemp[i][0];
+ pVector[i][1] = pTemp[i][1];
}
-
- for (i = 0; i < n; i++)
- pVector[i] = pTemp[i];
- delete [] pTemp;
-}
+ fftw_free(pTemp);
+}
+#endif