Isis 3 Programmer Reference
FourierTransform.cpp
Go to the documentation of this file.
1 
23 #include "FourierTransform.h"
24 
25 using namespace std;
26 
27 namespace Isis {
29  FourierTransform::FourierTransform() {};
30 
32  FourierTransform::~FourierTransform() {};
33 
42  std::vector< std::complex<double> >
43  FourierTransform::Transform(std::vector< std::complex<double> > input) {
44  // data length must be a power of two
45  // any extra space is filled with zeroes
46  int n = NextPowerOfTwo(input.size());
47  vector< std::complex<double> > output(n);
48  input.resize(n);
49 
50  // rearrange the data to fit the iterative algorithm
51  // which will apply the transform from the bottom up
52  for(int i = 0; i < n; i++) {
53  if(output[i] == 0.0) {
54  int j = BitReverse(n, i);
55  output[i] = input[j];
56  output[j] = input[i];
57  }
58  }
59 
60  // do the iterative fft calculation by first combining
61  // subarrays of length 2, then 4, 8, etc.
62  for(int m = 1; m < n; m *= 2) {
63  complex<double> Wm(polar(1.0, -1.0 * PI / m)); // Wm = e^(-PI/m *i)
64  for(int k = 0; k < n; k += 2 * m) {
65  // W = Wm^j, the roots of unity for x^m=1
66  complex<double> W(1.0);
67  for(int j = 0; j < m; j++) {
68  complex<double> t = W * output[k+j+m]; // the "twiddle" factor
69  complex<double> u = output[k+j];
70  output[k+j] = u + t; // a[k+j]+Wm^j*a[k+j+m]
71  output[k+j+m] = u - t; // a[k+j]+Wm^(j+m)*[k+j+m] = a[k+j]-Wm^j*[k+j+m]
72  W = W * Wm;
73  }
74  }
75  }
76 
77  return output;
78  }
79 
80 
89  std::vector< std::complex<double> >
90  FourierTransform::Inverse(std::vector< std::complex<double> > input) {
91  // Inverse(input) = 1/n*conj(Transform(conj(input)))
92  int n = input.size();
93  std::vector< std::complex<double> > temp(n);
94  for(int i = 0; i < n; i++) {
95  temp[i] = conj(input[i]);
96  }
97 
98  std::vector< std::complex<double> > output = Transform(temp);
99 
100  for(int i = 0; i < n; i++) {
101  output[i] = conj(output[i]) / ((double)n);
102  }
103 
104  return output;
105  }
106 
115  bool FourierTransform:: IsPowerOfTwo(int n) {
116  while(n > 1) {
117  if(n % 2 == 1) return false;
118  n /= 2;
119  }
120 
121  return true;
122  }
123 
131  int FourierTransform:: lg(int n) {
132  int k = 0;
133  while(n > 1) {
134  n = n / 2;
135  k++;
136  }
137  return k;
138  }
139 
152  int FourierTransform::BitReverse(int n, int x) {
153  double ans = 0.0;
154  double a = 1.0;
155  while(x > 0) {
156  if(x % 2 == 1) {
157  ans += a;
158  x -= 1;
159  }
160  x = x / 2;
161  a = a / 2.0;
162  }
163 
164  return(int)(ans * n / 2);
165  }
166 
175  int FourierTransform::NextPowerOfTwo(int n) {
176  if(IsPowerOfTwo(n)) return n;
177  return(int)pow(2.0, lg(n) + 1);
178  }
179 }
const double PI
The mathematical constant PI.
Definition: Constants.h:56
Namespace for the standard library.
Namespace for ISIS/Bullet specific routines.
Definition: Apollo.h:31
Pixel transformation.
Definition: Transform.h:89