Generative Fast Fourier Transforms (GFFT)  0.3
gfftswap.h
Go to the documentation of this file.
1 /***************************************************************************
2  * Copyright (C) 2006-2014 by Vladimir Mirnyy *
3  * *
4  * This program is free software; you can redistribute it and/or modify *
5  * it under the terms of the GNU General Public License as published by *
6  * the Free Software Foundation; either version 2 of the License, or *
7  * (at your option) any later version. *
8  * *
9  * This program is distributed in the hope that it will be useful, *
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
12  * GNU General Public License for more details. *
13  ***************************************************************************/
14 
15 #ifndef __gffswap_h
16 #define __gffswap_h
17 
23 namespace GFFT {
24 
25 using namespace MF;
26 
27 
28 
30 
39 template<uint_t M, uint_t P, typename T>
40 class SwapNR {
41  static const uint_t N = IPow<M,P>::value;
42 public:
43  void apply(T* data) {
44  int_t m, j = 0;
45  for (int_t i=0; i<2*N-1; i+=2) {
46  if (j>i) {
47  std::swap(data[j], data[i]);
48  std::swap(data[j+1], data[i+1]);
49  }
50  m = N;
51  while (m>=2 && j>=m) {
52  j -= m;
53  m >>= 1;
54  }
55  j += m;
56  }
57  }
58 };
59 
60 
62 
74 template<uint_t M, uint_t P, typename T, uint_t I=0>
75 class GFFTswap2 {
76  static const int_t BN = IPow<M,I>::value;
77  static const int_t BR = IPow<M,P-I-1>::value;
78  GFFTswap2<M,P,T,I+1> next;
79 public:
80  void apply(T* data, int_t n=0, int_t r=0) {
81  const int_t qn = n/BN;
82  const int_t rn = n%BN;
83  const int_t qr = r/BR;
84  const int_t rr = r%BR;
85  for (uint_t i = 0; i < M; ++i) {
86  next.apply(data,(qn+i)*BN+rn,(qr+i)*BR+rr);
87  }
88  }
89 };
90 
91 template<uint_t M, uint_t P, typename T>
92 class GFFTswap2<M,P,T,P> {
93 public:
94  void apply(T* data, int_t n=0, int_t r=0) {
95  if (n>r) {
96  const int_t n2 = 2*n;
97  const int_t r2 = 2*r;
98  std::swap(data[n2],data[r2]);
99  std::swap(data[n2+1],data[r2+1]);
100  }
101  }
102 };
103 
104 template<uint_t P, typename T, uint_t I>
105 class GFFTswap2<2,P,T,I> {
106  static const int_t BN = 1<<(I+1);
107  static const int_t BR = 1<<(P-I);
108  GFFTswap2<2,P,T,I+1> next;
109 public:
110  void apply(T* data, int_t n=0, int_t r=0) {
111  next.apply(data,n,r);
112  next.apply(data,n|BN,r|BR);
113  }
114 };
115 
116 template<uint_t P, typename T>
117 class GFFTswap2<2,P,T,P> {
118 public:
119  void apply(T* data, int_t n=0, int_t r=0) {
120  if (n>r) {
121  std::swap(data[n],data[r]);
122  std::swap(data[n+1],data[r+1]);
123  }
124  }
125 };
126 
128 
132 template<uint_t M, uint_t P, typename T,
133 template<typename> class Complex>
134 class SwapNR<M,P,Complex<T> >
135 {
136  static const uint_t N = IPow<M,P>::value;
137 public:
138  void apply(Complex<T>* data) {
139  int_t m,j=0;
140  for (int_t i=0; i<N; ++i) {
141  if (j>i) {
142  std::swap(data[j], data[i]);
143  }
144  m = N/2;
145  while (m>=1 && j>=m) {
146  j -= m;
147  m >>= 1;
148  }
149  j += m;
150  }
151  }
152 };
153 
154 template<uint_t M, uint_t P, typename T, uint_t I,
155 template<typename> class Complex>
156 class GFFTswap2<M,P,Complex<T>,I> {
157  static const int_t BN = IPow<M,I>::value;
158  static const int_t BR = IPow<M,P-I-1>::value;
159  GFFTswap2<M,P,Complex<T>,I+1> next;
160 public:
161  void apply(Complex<T>* data, int_t n=0, int_t r=0) {
162  const int_t qn = n/BN;
163  const int_t rn = n%BN;
164  const int_t qr = r/BR;
165  const int_t rr = r%BR;
166  for (uint_t i = 0; i < M; ++i) {
167  next.apply(data,(qn+i)*BN+rn,(qr+i)*BR+rr);
168  }
169  }
170 };
171 
172 template<uint_t M, uint_t P, typename T,
173 template<typename> class Complex>
174 class GFFTswap2<M,P,Complex<T>,P> {
175 public:
176  void apply(Complex<T>* data, int_t n=0, int_t r=0) {
177  if (n>r)
178  std::swap(data[n],data[r]);
179  }
180 };
181 
182 template<uint_t P, typename T, uint_t I,
183 template<typename> class Complex>
184 class GFFTswap2<2,P,Complex<T>,I> {
185  static const int_t BN = 1<<I;
186  static const int_t BR = 1<<(P-I-1);
187  GFFTswap2<2,P,Complex<T>,I+1> next;
188 public:
189  void apply(Complex<T>* data, int_t n=0, int_t r=0) {
190  next.apply(data,n,r);
191  next.apply(data,n|BN,r|BR);
192  }
193 };
194 
195 template<uint_t P, typename T,
196 template<typename> class Complex>
197 class GFFTswap2<2,P,Complex<T>,P> {
198 public:
199  void apply(Complex<T>* data, int_t n=0, int_t r=0) {
200  if (n>r)
201  swap(data[n],data[r]);
202  }
203 };
204 
205 
207 
212 template<int_t N, typename T, int S>
213 class Separate {
214  typedef typename TempTypeTrait<T>::Result LocalVType;
215  static const int M = (S==1) ? 2 : 1;
216 public:
217  void apply(T* data) {
218  int_t i,i1,i2,i3,i4;
219  LocalVType wtemp,wr,wi,wpr,wpi;
220  LocalVType h1r,h1i,h2r,h2i,h3r,h3i;
221  wtemp = Sin<2*N,1,LocalVType>::value();
222  wpr = -2.*wtemp*wtemp;
223  wpi = -S*Sin<N,1,LocalVType>::value();
224  wr = 1.+wpr;
225  wi = wpi;
226  for (i=1; i<N/2; ++i) {
227  i1 = i+i;
228  i2 = i1+1;
229  i3 = 2*N-i1;
230  i4 = i3+1;
231  h1r = 0.5*(data[i1]+data[i3]);
232  h1i = 0.5*(data[i2]-data[i4]);
233  h2r = S*0.5*(data[i2]+data[i4]);
234  h2i =-S*0.5*(data[i1]-data[i3]);
235  h3r = wr*h2r - wi*h2i;
236  h3i = wr*h2i + wi*h2r;
237  data[i1] = h1r + h3r;
238  data[i2] = h1i + h3i;
239  data[i3] = h1r - h3r;
240  data[i4] =-h1i + h3i;
241 
242  wtemp = wr;
243  wr += wr*wpr - wi*wpi;
244  wi += wi*wpr + wtemp*wpi;
245  }
246  h1r = data[0];
247  data[0] = M*0.5*(h1r + data[1]);
248  data[1] = M*0.5*(h1r - data[1]);
249 
250  if (N>1) data[N+1] = -data[N+1];
251  }
252 
253  void apply(const T*, T*)
254  {
255  std::cout << "Not implemented!" << std::endl;
256  exit(1);
257  }
258 };
259 
260 
262 
267 template<int_t N, typename T, int S,
268 template<typename> class Complex>
269 class Separate<N,Complex<T>,S> {
270  typedef typename TempTypeTrait<T>::Result LocalVType;
271  typedef Complex<LocalVType> LocalComplex;
272  static const int M = (S==1) ? 2 : 1;
273 public:
274  void apply(Complex<T>* data) {
275  int_t i,i1;
276  LocalComplex h1,h2,h3;
277  LocalVType wtemp = Sin<2*N,1,LocalVType>::value();
278  LocalComplex wp(-2.*wtemp*wtemp,-S*Sin<N,1,LocalVType>::value());
279  LocalComplex w(1.+wp.real(),wp.imag());
280 
281  for (i=1; i<N/2; ++i) {
282  i1 = N-i;
283  h1 = Complex<LocalVType>(static_cast<LocalVType>(0.5*(data[i].real()+data[i1].real())),
284  static_cast<LocalVType>(0.5*(data[i].imag()-data[i1].imag())));
285  h2 = Complex<LocalVType>(static_cast<LocalVType>( S*0.5*(data[i].imag()+data[i1].imag())),
286  static_cast<LocalVType>(-S*0.5*(data[i].real()-data[i1].real())));
287  h3 = w*h2;
288  data[i] = h1 + h3;
289  data[i1]= h1 - h3;
290  data[i1] = Complex<T>(data[i1].real(), -data[i1].imag());
291 
292  w += w*wp;
293  }
294  wtemp = data[0].real();
295  data[0] = Complex<T>(M*0.5*(wtemp + data[0].imag()), M*0.5*(wtemp - data[0].imag()));
296 
297  data[N/2] = Complex<T>(data[N/2].real(), -data[N/2].imag());
298  }
299 
300  void apply(const Complex<T>*, Complex<T>*)
301  {
302  std::cout << "Not implemented!" << std::endl;
303  exit(1);
304  }
305 };
306 
307 // Policy for a definition of forward FFT
308 template<int_t N, typename T>
309 struct Forward {
310  enum { Sign = 1 };
311  void apply(T*) { }
312  void apply(const T*, T*) { }
313  void apply(const T*, T*, T*) { }
314 };
315 
316 template<int_t N, typename T,
317 template<typename> class Complex>
318 struct Forward<N,Complex<T> > {
319  enum { Sign = 1 };
320  void apply(Complex<T>*) { }
321  void apply(const Complex<T>*, Complex<T>*) { }
322  void apply(const Complex<T>*, Complex<T>*, Complex<T>*) { }
323 };
324 
325 // Policy for a definition of backward FFT
326 template<int_t N, typename T>
327 struct Backward {
328  enum { Sign = -1 };
329  void apply(T* data) {
330  for (T* i=data; i<data+2*N; ++i) *i/=N;
331  }
332  void apply(const T*, T* dst) { apply(dst); }
333  void apply(const T*, T* dst, T*) { apply(dst); }
334 };
335 
336 template<int_t N, typename T,
337 template<typename> class Complex>
338 struct Backward<N,Complex<T> > {
339  enum { Sign = -1 };
340  void apply(Complex<T>* data) {
341  for (int_t i=0; i<N; ++i) {
342  data[i]/=N;
343  }
344  }
345  void apply(const Complex<T>*, Complex<T>* dst) { apply(dst); }
346  void apply(const Complex<T>*, Complex<T>* dst, Complex<T>*) { apply(dst); }
347 };
348 
349 
350 } //namespace DFT
351 
352 #endif /*__gfftalg_h*/
Generated on Mon Feb 10 2014 for Generative Fast Fourier Transforms (GFFT) by DoxyGen 1.8.3.1