summaryrefslogtreecommitdiffstats
path: root/include
diff options
context:
space:
mode:
authorDaniel M. Pelt <D.M.Pelt@cwi.nl>2016-04-22 14:22:53 +0200
committerDaniel M. Pelt <D.M.Pelt@cwi.nl>2016-04-22 14:27:15 +0200
commit6ccde536191676f9b504055b16c68786858b693d (patch)
tree581e747a80ae5431dd6ad3d54d55dbbe0df84e75 /include
parent547def0ea6e3eab07b7e4c48cee6d6a81f6155e1 (diff)
downloadastra-6ccde536191676f9b504055b16c68786858b693d.tar.gz
astra-6ccde536191676f9b504055b16c68786858b693d.tar.bz2
astra-6ccde536191676f9b504055b16c68786858b693d.tar.xz
astra-6ccde536191676f9b504055b16c68786858b693d.zip
Change CPU FFT implementation
Diffstat (limited to 'include')
-rw-r--r--include/astra/Fourier.h132
1 files changed, 44 insertions, 88 deletions
diff --git a/include/astra/Fourier.h b/include/astra/Fourier.h
index b515dc6..ff26f82 100644
--- a/include/astra/Fourier.h
+++ b/include/astra/Fourier.h
@@ -33,94 +33,50 @@ $Id$
namespace astra {
-
-/**
- * Perform a 1D DFT or inverse DFT.
- *
- * @param iLength number of elements
- * @param pfRealIn real part of input
- * @param pfImaginaryIn imaginary part of input
- * @param pfRealOut real part of output
- * @param pfImaginaryOut imaginary part of output
- * @param iStrideIn distance between elements in pf*In
- * @param iStrideOut distance between elements in pf*Out
- * @param bInverse if true, perform an inverse DFT
- */
-
-void _AstraExport discreteFourierTransform1D(unsigned int iLength,
- const float32* pfRealIn,
- const float32* pfImaginaryIn,
- float32* pfRealOut,
- float32* pfImaginaryOut,
- unsigned int iStrideIn,
- unsigned int iStrideOut,
- bool bInverse);
-
-/**
- * Perform a 2D DFT or inverse DFT.
- *
- * @param iHeight number of rows
- * @param iWidth number of columns
- * @param pfRealIn real part of input
- * @param pfImaginaryIn imaginary part of input
- * @param pfRealOut real part of output
- * @param pfImaginaryOut imaginary part of output
- * @param bInverse if true, perform an inverse DFT
- */
-
-void _AstraExport discreteFourierTransform2D(unsigned int iHeight, unsigned int iWidth,
- const float32* pfRealIn,
- const float32* pfImaginaryIn,
- float32* pfRealOut,
- float32* pfImaginaryOut,
- bool bInverse);
-
-/**
- * Perform a 1D FFT or inverse FFT. The size must be a power of two.
- * This transform can be done in-place, so the input and output pointers
- * may point to the same data.
- *
- * @param iLength number of elements, must be a power of two
- * @param pfRealIn real part of input
- * @param pfImaginaryIn imaginary part of input
- * @param pfRealOut real part of output
- * @param pfImaginaryOut imaginary part of output
- * @param iStrideIn distance between elements in pf*In
- * @param iStrideOut distance between elements in pf*Out
- * @param bInverse if true, perform an inverse DFT
- */
-
-void _AstraExport fastTwoPowerFourierTransform1D(unsigned int iLength,
- const float32* pfRealIn,
- const float32* pfImaginaryIn,
- float32* pfRealOut,
- float32* pfImaginaryOut,
- unsigned int iStrideIn,
- unsigned int iStrideOut,
- bool bInverse);
-
-/**
- * Perform a 2D FFT or inverse FFT. The size must be a power of two.
- * This transform can be done in-place, so the input and output pointers
- * may point to the same data.
- *
- * @param iHeight number of rows, must be a power of two
- * @param iWidth number of columns, must be a power of two
- * @param pfRealIn real part of input
- * @param pfImaginaryIn imaginary part of input
- * @param pfRealOut real part of output
- * @param pfImaginaryOut imaginary part of output
- * @param bInverse if true, perform an inverse DFT
- */
-
-void _AstraExport fastTwoPowerFourierTransform2D(unsigned int iHeight,
- unsigned int iWidth,
- const float32* pfRealIn,
- const float32* pfImaginaryIn,
- float32* pfRealOut,
- float32* pfImaginaryOut,
- bool bInverse);
-
+/*
+-------- Complex DFT (Discrete Fourier Transform) --------
+ [definition]
+ <case1>
+ X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
+ <case2>
+ X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
+ (notes: sum_j=0^n-1 is a summation from j=0 to n-1)
+ [usage]
+ <case1>
+ ip[0] = 0; // first time only
+ cdft(2*n, 1, a, ip, w);
+ <case2>
+ ip[0] = 0; // first time only
+ cdft(2*n, -1, a, ip, w);
+ [parameters]
+ 2*n :data length (int)
+ n >= 1, n = power of 2
+ a[0...2*n-1] :input/output data (float32 *)
+ input data
+ a[2*j] = Re(x[j]),
+ a[2*j+1] = Im(x[j]), 0<=j<n
+ output data
+ a[2*k] = Re(X[k]),
+ a[2*k+1] = Im(X[k]), 0<=k<n
+ ip[0...*] :work area for bit reversal (int *)
+ length of ip >= 2+sqrt(n)
+ strictly,
+ length of ip >=
+ 2+(1<<(int)(log(n+0.5)/log(2))/2).
+ ip[0],ip[1] are pointers of the cos/sin table.
+ w[0...n/2-1] :cos/sin table (float32 *)
+ w[],ip[] are initialized if ip[0] == 0.
+ [remark]
+ Inverse of
+ cdft(2*n, -1, a, ip, w);
+ is
+ cdft(2*n, 1, a, ip, w);
+ for (j = 0; j <= 2 * n - 1; j++) {
+ a[j] *= 1.0 / n;
+ }
+ .
+*/
+void cdft(int n, int isgn, float32 *a, int *ip, float32 *w);
}