Generative Fast Fourier Transforms (GFFT)  0.3
gfftalg.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 __gfftalg_h
16 #define __gfftalg_h
17 
22 #include "gfftspec.h"
23 #include "gfftfactor.h"
24 #include "gfftswap.h"
25 
26 #include "metacomplex.h"
27 #include "metaroot.h"
28 
29 namespace GFFT {
30 
31 using namespace MF;
32 
33 
34 static const int_t StaticLoopLimit = 8;
35 
36 
37 template<int_t K, int_t M, typename T, int S, class W1, int NIter = 1, class W = W1>
38 class IterateInTime
39 {
40  typedef typename TempTypeTrait<T>::Result LocalVType;
41  typedef Compute<typename W::Re,2> WR;
42  typedef Compute<typename W::Im,2> WI;
43  static const int_t M2 = M*2;
44  static const int_t N = K*M;
45 
46  typedef typename GetNextRoot<NIter+1,N,W1,W,2>::Result Wnext;
47  IterateInTime<K,M,T,S,W1,NIter+1,Wnext> next;
48  DFTk_inp<K,M2,T,S> spec_inp;
49 
50 public:
51  void apply(T* data)
52  {
53  const LocalVType wr = WR::value();
54  const LocalVType wi = WI::value();
55 
56  spec_inp.apply(data + (NIter-1)*2, &wr, &wi);
57 
58  next.apply(data);
59  }
60 };
61 
62 // Last step of the loop
63 template<int_t K, int_t M, typename T, int S, class W1, class W>
64 class IterateInTime<K,M,T,S,W1,M,W>
65 {
66 // typedef typename RList::Head H;
67  typedef typename TempTypeTrait<T>::Result LocalVType;
68  typedef Compute<typename W::Re,2> WR;
69  typedef Compute<typename W::Im,2> WI;
70  static const int_t M2 = M*2;
71  static const int_t N = K*M;
72  DFTk_inp<K,M2,T,S> spec_inp;
73 public:
74  void apply(T* data)
75  {
76  const LocalVType wr = WR::value();
77  const LocalVType wi = WI::value();
78 
79  spec_inp.apply(data + (M-1)*2, &wr, &wi);
80  }
81 };
82 
83 // First step in the loop
84 template<int_t K, int_t M, typename T, int S, class W1, class W>
85 class IterateInTime<K,M,T,S,W1,1,W> {
86  static const int_t M2 = M*2;
87  DFTk_inp<K,M2,T,S> spec_inp;
88  IterateInTime<K,M,T,S,W1,2,W> next;
89 public:
90  void apply(T* data)
91  {
92  spec_inp.apply(data);
93  next.apply(data);
94  }
95 };
96 
97 
99 
112 template<int_t K, int_t M, typename T, int S, class W, bool doStaticLoop>
113 class DFTk_x_Im_T;
114 
115 // Rely on the static template loop
116 template<int_t K, int_t M, typename T, int S, class W>
117 class DFTk_x_Im_T<K,M,T,S,W,true> : public IterateInTime<K,M,T,S,W> {};
118 
119 // General implementation
120 template<int_t K, int_t M, typename T, int S, class W>
121 class DFTk_x_Im_T<K,M,T,S,W,false>
122 {
123  typedef typename TempTypeTrait<T>::Result LocalVType;
124  typedef Compute<typename W::Re,2> WR;
125  typedef Compute<typename W::Im,2> WI;
126  static const int_t N = K*M;
127  static const int_t M2 = M*2;
128  DFTk_inp<K,M2,T,S> spec_inp;
129 public:
130  void apply(T* data)
131  {
132  spec_inp.apply(data);
133 
134  LocalVType wr[K-1], wi[K-1], wpr[K-1], wpi[K-1], t;
135 
136  // W = (wpr[0], wpi[0])
137  wpr[0] = WR::value();
138  wpi[0] = WI::value();
139  //t = Sin<N,1,LocalVType>::value();
140 // wpr[0] = 1 - 2.0*t*t;
141 // wpi[0] = -S*Sin<N,2,LocalVType>::value();
142 
143  // W^i = (wpr[i], wpi[i])
144  for (int_t i=0; i<K-2; ++i) {
145  wpr[i+1] = wpr[i]*wpr[0] - wpi[i]*wpi[0];
146  wpi[i+1] = wpr[i]*wpi[0] + wpr[0]*wpi[i];
147  }
148 
149  for (int_t i=0; i<K-1; ++i) {
150  wr[i] = wpr[i];
151  wi[i] = wpi[i];
152  }
153 
154  for (int_t i=2; i<M2; i+=2) {
155  spec_inp.apply(data+i, wr, wi);
156 
157  for (int_t i=0; i<K-1; ++i) {
158  t = wr[i];
159  wr[i] = t*wpr[i] - wi[i]*wpi[i];
160  wi[i] = wi[i]*wpr[i] + t*wpi[i];
161  }
162  }
163  }
164 
165 };
166 
167 // Specialization for radix 3
168 template<int_t M, typename T, int S, class W>
169 class DFTk_x_Im_T<3,M,T,S,W,false> {
170  typedef typename TempTypeTrait<T>::Result LocalVType;
171  typedef Compute<typename W::Re,2> WR;
172  typedef Compute<typename W::Im,2> WI;
173  static const int_t N = 3*M;
174  static const int_t M2 = M*2;
175  DFTk_inp<3,M2,T,S> spec_inp;
176 public:
177  void apply(T* data)
178  {
179  spec_inp.apply(data);
180 
181  LocalVType wr[2],wi[2],t;
182 
183  // W = (wpr1, wpi1)
184 // t = Sin<N,1,LocalVType>::value();
185 // const LocalVType wpr1 = 1 - 2.0*t*t;
186 // const LocalVType wpi1 = -S*Sin<N,2,LocalVType>::value();
187  const LocalVType wpr1 = WR::value();
188  const LocalVType wpi1 = WI::value();
189 
190  // W^2 = (wpr2, wpi2)
191  const LocalVType wpr2 = wpr1*wpr1 - wpi1*wpi1;
192  const LocalVType wpi2 = 2*wpr1*wpi1;
193 
194  wr[0] = wpr1;
195  wi[0] = wpi1;
196  wr[1] = wpr2;
197  wi[1] = wpi2;
198  for (int_t i=2; i<M2; i+=2) {
199  spec_inp.apply(data+i, wr, wi);
200 
201  t = wr[0];
202  wr[0] = t*wpr1 - wi[0]*wpi1;
203  wi[0] = wi[0]*wpr1 + t*wpi1;
204  t = wr[1];
205  wr[1] = t*wpr2 - wi[1]*wpi2;
206  wi[1] = wi[1]*wpr2 + t*wpi2;
207  }
208  }
209 };
210 
211 // Specialization for radix 2
212 template<int_t M, typename T, int S, class W>
213 class DFTk_x_Im_T<2,M,T,S,W,false> {
214  typedef typename TempTypeTrait<T>::Result LocalVType;
215  typedef Compute<typename W::Re,2> WR;
216  typedef Compute<typename W::Im,2> WI;
217  static const int_t N = 2*M;
218  DFTk_inp<2,N,T,S> spec_inp;
219 public:
220  void apply(T* data)
221  {
222  spec_inp.apply(data);
223 
224  LocalVType wr,wi,t;
225 // t = Sin<N,1,LocalVType>::value();
226 // const LocalVType wpr = 1-2.0*t*t;
227 // const LocalVType wpi = -S*Sin<N,2,LocalVType>::value();
228  const LocalVType wpr = WR::value();
229  const LocalVType wpi = WI::value();
230  wr = wpr;
231  wi = wpi;
232  for (int_t i=2; i<N; i+=2) {
233  spec_inp.apply(data+i, &wr, &wi);
234 
235  t = wr;
236  wr = wr*wpr - wi*wpi;
237  wi = wi*wpr + t*wpi;
238  }
239  }
240 };
241 
243 
256 template<int_t N, typename NFact, typename T, int S, class W1, int_t LastK = 1>
257 class InTime;
258 
259 // template<int_t N, typename Head, typename Tail, typename T, int S, class W1, int_t LastK>
260 // class InTime<N, Loki::Typelist<Head,Tail>, T, S, W1, LastK>
261 // {
262 // // Not implemented, because not allowed
263 // // Transforms in-place are allowed for powers of primes only!!!
264 // };
265 
266 template<int_t N, typename Head, typename T, int S, class W1, int_t LastK>
267 class InTime<N, Loki::Typelist<Head,Loki::NullType>, T, S, W1, LastK>
268 {
269  typedef typename TempTypeTrait<T>::Result LocalVType;
270  static const int_t K = Head::first::value;
271  static const int_t M = N/K;
272  static const int_t M2 = M*2;
273  static const int_t N2 = N*2;
274 
275  typedef typename IPowBig<W1,K>::Result WK;
276  typedef Loki::Typelist<Pair<typename Head::first, SInt<Head::second::value-1> >, Loki::NullType> NFactNext;
277  InTime<M,NFactNext,T,S,WK,K*LastK> dft_str;
278  DFTk_x_Im_T<K,M,T,S,W1,(N<=StaticLoopLimit)> dft_scaled;
279 // DFTk_x_Im_T<K,M,T,S,W1,false> dft_scaled;
280 public:
281  void apply(T* data)
282  {
283  // run strided DFT recursively K times
284  for (int_t m=0; m < N2; m+=M2)
285  dft_str.apply(data + m);
286 
287  dft_scaled.apply(data);
288  }
289 };
290 
291 // Take the next factor from the list
292 template<int_t N, int_t K, typename Tail, typename T, int S, class W1, int_t LastK>
293 class InTime<N, Loki::Typelist<Pair<SInt<K>, SInt<0> >,Tail>, T, S, W1, LastK>
294 : public InTime<N, Tail, T, S, W1, LastK> {};
295 
296 
297 // Specialization for a prime N
298 template<int_t N, typename T, int S, class W1, int_t LastK>
299 class InTime<N,Loki::Typelist<Pair<SInt<N>, SInt<1> >, Loki::NullType>,T,S,W1,LastK> {
300  DFTk_inp<N, 2, T, S> spec_inp;
301 public:
302  void apply(T* data)
303  {
304  spec_inp.apply(data);
305  }
306 };
307 
308 
310 
323 template<int_t N, typename NFact, typename T, int S, class W1, int_t LastK = 1>
324 class InTimeOOP;
325 
326 template<int_t N, typename Head, typename Tail, typename T, int S, class W1, int_t LastK>
327 class InTimeOOP<N, Loki::Typelist<Head,Tail>, T, S, W1, LastK>
328 {
329  typedef typename TempTypeTrait<T>::Result LocalVType;
330  static const int_t K = Head::first::value;
331  static const int_t M = N/K;
332  static const int_t M2 = M*2;
333  static const int_t N2 = N*2;
334  static const int_t LastK2 = LastK*2;
335 
336  typedef typename IPowBig<W1,K>::Result WK;
337  typedef Loki::Typelist<Pair<typename Head::first, SInt<Head::second::value-1> >, Tail> NFactNext;
338  InTimeOOP<M,NFactNext,T,S,WK,K*LastK> dft_str;
339  DFTk_x_Im_T<K,M,T,S,W1,(N<=StaticLoopLimit)> dft_scaled;
340 // DFTk_x_Im_T<K,M,T,S,W1,false> dft_scaled;
341 public:
342 
343  void apply(const T* src, T* dst)
344  {
345  // run strided DFT recursively K times
346  int_t lk = 0;
347  for (int_t m = 0; m < N2; m+=M2, lk+=LastK2)
348  dft_str.apply(src + lk, dst + m);
349 
350  dft_scaled.apply(dst);
351  }
352 };
353 
354 // Take the next factor from the list
355 template<int_t N, int_t K, typename Tail, typename T, int S, class W1, int_t LastK>
356 class InTimeOOP<N, Loki::Typelist<Pair<SInt<K>, SInt<0> >,Tail>, T, S, W1, LastK>
357 : public InTimeOOP<N, Tail, T, S, W1, LastK> {};
358 
359 
360 // Specialization for prime N
361 template<int_t N, typename T, int S, class W1, int_t LastK>
362 class InTimeOOP<N,Loki::Typelist<Pair<SInt<N>, SInt<1> >, Loki::NullType>,T,S,W1,LastK>
363 : public DFTk<N, LastK*2, 2, T, S> {};
364 
365 
366 } //namespace DFT
367 
368 #endif /*__gfftalg_h*/
Generated on Mon Feb 10 2014 for Generative Fast Fourier Transforms (GFFT) by DoxyGen 1.8.3.1