Generative Fast Fourier Transforms (GFFT)  0.3
gfftalgfreq.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 __gfftalgfreq_h
16 #define __gfftalgfreq_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 template<int_t K, int_t M, typename T, int S, class W1, int NIter = 1, class W = W1>
35 class IterateInFreq
36 {
37  typedef typename TempTypeTrait<T>::Result LocalVType;
38  typedef Compute<typename W::Re,2> WR;
39  typedef Compute<typename W::Im,2> WI;
40  static const int_t M2 = M*2;
41  static const int_t N = K*M;
42 
43  typedef typename GetNextRoot<NIter+1,N,W1,W,2>::Result Wnext;
44  IterateInFreq<K,M,T,S,W1,NIter+1,Wnext> next;
45  DFTk_inp<K,M2,T,S> spec_inp;
46 
47 public:
48  void apply(T* data)
49  {
50  const LocalVType wr = WR::value();
51  const LocalVType wi = WI::value();
52 
53  spec_inp.apply(&wr, &wi, data + (NIter-1)*2);
54 
55  next.apply(data);
56  }
57 };
58 
59 // Last step of the loop
60 template<int_t K, int_t M, typename T, int S, class W1, class W>
61 class IterateInFreq<K,M,T,S,W1,M,W>
62 {
63 // typedef typename RList::Head H;
64  typedef typename TempTypeTrait<T>::Result LocalVType;
65  typedef Compute<typename W::Re,2> WR;
66  typedef Compute<typename W::Im,2> WI;
67  static const int_t M2 = M*2;
68  static const int_t N = K*M;
69  DFTk_inp<K,M2,T,S> spec_inp;
70 public:
71  void apply(T* data)
72  {
73 // const LocalVType t = Sin<N,M-1,LocalVType>::value();
74 // const LocalVType wr = 1 - 2.0*t*t;
75 // const LocalVType wi = -S*Sin<N,2*(M-1),LocalVType>::value();
76  const LocalVType wr = WR::value();
77  const LocalVType wi = WI::value();
78 
79  spec_inp.apply(&wr, &wi, data + (M-1)*2);
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 IterateInFreq<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  IterateInFreq<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 
98 template<int_t K, int_t M, typename T, int S, class W, bool doStaticLoop>
99 class T_DFTk_x_Im;
100 
101 template<int_t K, int_t M, typename T, int S, class W>
102 class T_DFTk_x_Im<K,M,T,S,W,true>
103 {
104  IterateInFreq<K,M,T,S,W> iterate;
105 public:
106  void apply(T* data)
107  {
108  iterate.apply(data);
109  }
110 };
111 
112 template<int_t K, int_t M, typename T, int S, class W>
113 class T_DFTk_x_Im<K,M,T,S,W,false>
114 {
115  typedef typename TempTypeTrait<T>::Result LocalVType;
116  static const int_t N = K*M;
117  static const int_t M2 = M*2;
118  DFTk_inp<K,M2,T,S> spec_inp;
119 public:
120  void apply(T* data)
121  {
122  spec_inp.apply(data);
123 
124  LocalVType wr[K-1], wi[K-1], wpr[K-1], wpi[K-1], t;
125  t = Sin<N,1,LocalVType>::value();
126 
127  // W = (wpr[0], wpi[0])
128  wpr[0] = 1 - 2.0*t*t;
129  wpi[0] = -S*Sin<N,2,LocalVType>::value();
130 
131  // W^i = (wpr2, wpi2)
132  for (int_t i=0; i<K-2; ++i) {
133  wpr[i+1] = wpr[i]*wpr[0] - wpi[i]*wpi[0];
134  wpi[i+1] = wpr[i]*wpi[0] + wpr[0]*wpi[i];
135  }
136 
137  for (int_t i=0; i<K-1; ++i) {
138  wr[i] = wpr[i];
139  wi[i] = wpi[i];
140  }
141 
142  for (int_t i=2; i<M2; i+=2) {
143  spec_inp.apply(wr, wi, data+i);
144 
145  for (int_t i=0; i<K-1; ++i) {
146  t = wr[i];
147  wr[i] = t*wpr[i] - wi[i]*wpi[i];
148  wi[i] = wi[i]*wpr[i] + t*wpi[i];
149  }
150  }
151  }
152 };
153 
154 template<int_t M, typename T, int S, class W>
155 class T_DFTk_x_Im<3,M,T,S,W,false> {
156  typedef typename TempTypeTrait<T>::Result LocalVType;
157  static const int_t N = 3*M;
158  static const int_t M2 = M*2;
159  DFTk_inp<3,M2,T,S> spec_inp;
160 public:
161  void apply(T* data)
162  {
163  spec_inp.apply(data);
164 
165  LocalVType wr[2],wi[2],t;
166  t = Sin<N,1,LocalVType>::value();
167 
168  // W = (wpr1, wpi1)
169  const LocalVType wpr1 = 1 - 2.0*t*t;
170  const LocalVType wpi1 = -S*Sin<N,2,LocalVType>::value();
171 
172  // W^2 = (wpr2, wpi2)
173  const LocalVType wpr2 = wpr1*wpr1 - wpi1*wpi1;
174  const LocalVType wpi2 = 2*wpr1*wpi1;
175 
176  wr[0] = wpr1;
177  wi[0] = wpi1;
178  wr[1] = wpr2;
179  wi[1] = wpi2;
180  for (int_t i=2; i<M2; i+=2) {
181  spec_inp.apply(wr, wi, data+i);
182 
183  t = wr[0];
184  wr[0] = t*wpr1 - wi[0]*wpi1;
185  wi[0] = wi[0]*wpr1 + t*wpi1;
186  t = wr[1];
187  wr[1] = t*wpr2 - wi[1]*wpi2;
188  wi[1] = wi[1]*wpr2 + t*wpi2;
189  }
190  }
191 };
192 
193 template<int_t M, typename T, int S, class W>
194 class T_DFTk_x_Im<2,M,T,S,W,false>
195 {
196  typedef typename TempTypeTrait<T>::Result LocalVType;
197  static const int_t N = 2*M;
198  DFTk_inp<2,N,T,S> spec_inp;
199 
200 // IterateInFreq<2,M,T,S,W> iterate;
201 public:
202  void apply(T* data)
203  {
204 // iterate.apply(data);
205  spec_inp.apply(data);
206 
207  LocalVType t,wr,wi;
208  t = Sin<N,1,LocalVType>::value();
209  const LocalVType wpr = -2.0*t*t;
210  const LocalVType wpi = -S*Sin<N,2,LocalVType>::value();
211  wr = 1+wpr;
212  wi = wpi;
213  for (int_t i=2; i<N; i+=2) {
214 //std::cout << i << "/" << N << ": " << t << " --- (" << wr << ", " << wi << ")" << std::endl;
215  spec_inp.apply(&wr, &wi, data+i);
216 
217  t = wr;
218  wr += t*wpr - wi*wpi;
219  wi += wi*wpr + t*wpi;
220  }
221  }
222 };
223 
224 
225 /*
226 template<int_t K, int_t M, typename T, int S, class W, bool doStaticLoop>
227 class T_DFTk_x_Im_OOP;
228 
229 template<int_t K, int_t M, typename T, int S, class W>
230 class T_DFTk_x_Im_OOP<K,M,T,S,W,true>
231 {
232  IterateInFreqOOP<K,M,T,S,W> iterate;
233 public:
234  void apply(const T* src, T* dst)
235  {
236  iterate.apply(src,dst);
237  }
238 };
239 
240 template<int_t K, int_t M, typename T, int S, class W>
241 class T_DFTk_x_Im_OOP<K,M,T,S,W,false>
242 {
243  typedef typename TempTypeTrait<T>::Result LocalVType;
244  static const int_t N = K*M;
245  static const int_t M2 = M*2;
246  DFTk<K,M2,M2,T,S> spec;
247 public:
248  void apply(const T* src, T* dst)
249  {
250  spec.apply(src, dst);
251 
252  LocalVType wr[K-1], wi[K-1], wpr[K-1], wpi[K-1], t;
253  t = Sin<N,1,LocalVType>::value();
254 
255  // W = (wpr[0], wpi[0])
256  wpr[0] = 1 - 2.0*t*t;
257  wpi[0] = -S*Sin<N,2,LocalVType>::value();
258 
259  // W^i = (wpr2, wpi2)
260  for (int_t i=0; i<K-2; ++i) {
261  wpr[i+1] = wpr[i]*wpr[0] - wpi[i]*wpi[0];
262  wpi[i+1] = wpr[i]*wpi[0] + wpr[0]*wpi[i];
263  }
264 
265  for (int_t i=0; i<K-1; ++i) {
266  wr[i] = wpr[i];
267  wi[i] = wpi[i];
268  }
269  for (int_t i=2; i<M2; i+=2) {
270  spec.apply(wr, wi, src+i, dst+i);
271 
272  for (int_t i=0; i<K-1; ++i) {
273  t = wr[i];
274  wr[i] = t*wpr[i] - wi[i]*wpi[i];
275  wi[i] = wi[i]*wpr[i] + t*wpi[i];
276  }
277  }
278  }
279 };
280 
281 template<int_t M, typename T, int S, class W>
282 class T_DFTk_x_Im_OOP<3,M,T,S,W,false> {
283  typedef typename TempTypeTrait<T>::Result LocalVType;
284  static const int_t N = 3*M;
285  static const int_t M2 = M*2;
286  DFTk<3,M2,M2,T,S> spec;
287 public:
288  void apply(const T* src, T* dst)
289  {
290  spec.apply(src, dst);
291 
292  LocalVType wr[2],wi[2],t;
293  t = Sin<N,1,LocalVType>::value();
294 
295  // W = (wpr1, wpi1)
296  const LocalVType wpr1 = 1 - 2.0*t*t;
297  const LocalVType wpi1 = -S*Sin<N,2,LocalVType>::value();
298 
299  // W^2 = (wpr2, wpi2)
300  const LocalVType wpr2 = wpr1*wpr1 - wpi1*wpi1;
301  const LocalVType wpi2 = 2*wpr1*wpi1;
302 
303  wr[0] = wpr1;
304  wi[0] = wpi1;
305  wr[1] = wpr2;
306  wi[1] = wpi2;
307  for (int_t i=2; i<M2; i+=2) {
308  spec.apply(wr, wi, src+i, dst+i);
309 
310  t = wr[0];
311  wr[0] = t*wpr1 - wi[0]*wpi1;
312  wi[0] = wi[0]*wpr1 + t*wpi1;
313  t = wr[1];
314  wr[1] = t*wpr2 - wi[1]*wpi2;
315  wi[1] = wi[1]*wpr2 + t*wpi2;
316  }
317  }
318 };
319 
320 template<int_t M, typename T, int S, class W>
321 class T_DFTk_x_Im_OOP<2,M,T,S,W,false>
322 {
323  typedef typename TempTypeTrait<T>::Result LocalVType;
324  static const int_t N = 2*M;
325  DFTk<2,N,N,T,S> spec;
326 
327 // IterateInFreq<2,M,T,S,W> iterate;
328 public:
329  void apply(const T* src, T* dst)
330  {
331  spec.apply(src, dst);
332 
333  LocalVType t,wr,wi;
334  t = Sin<N,1,LocalVType>::value();
335  const LocalVType wpr = -2.0*t*t;
336  const LocalVType wpi = -S*Sin<N,2,LocalVType>::value();
337  wr = 1+wpr;
338  wi = wpi;
339  for (int_t i=2; i<N; i+=2) {
340  spec.apply(&wr, &wi, src+i, dst+i);
341 
342  t = wr;
343  wr += t*wpr - wi*wpi;
344  wi += wi*wpr + t*wpi;
345  }
346  }
347 };
348 */
350 
362 template<int_t N, typename NFact, typename T, int S, class W1, int_t LastK = 1>
363 class InFreq;
364 
365 // template<int_t N, typename Head, typename Tail, typename T, int S, class W1, int_t LastK>
366 // class InFreq<N, Loki::Typelist<Head,Tail>, T, S, W1, LastK>
367 // {
368 // // Not implemented, because not allowed
369 // // Transforms in-place are allowed for powers of primes only!!!
370 // };
371 
372 template<int_t N, typename Head, typename T, int S, class W1, int_t LastK>
373 class InFreq<N, Loki::Typelist<Head,Loki::NullType>, T, S, W1, LastK>
374 {
375  typedef typename TempTypeTrait<T>::Result LocalVType;
376  static const int_t K = Head::first::value;
377  static const int_t M = N/K;
378  static const int_t M2 = M*2;
379  static const int_t N2 = N*2;
380  static const int_t LastK2 = LastK*2;
381 
382 // typedef typename Loki::TL::Next<RList,K-1>::Result RListK;
383  typedef typename IPowBig<W1,K>::Result WK;
384  typedef Loki::Typelist<Pair<typename Head::first, SInt<Head::second::value-1> >, Loki::NullType> NFactNext;
385  InFreq<M,NFactNext,T,S,WK,K*LastK> dft_str;
386  T_DFTk_x_Im<K,M,T,S,W1,(N<=StaticLoopLimit)> dft_scaled;
387 
388 public:
389  void apply(T* data)
390  {
391  dft_scaled.apply(data);
392 
393  // K times call to dft_str.apply()
394  for (int_t m = 0; m < N2; m+=M2)
395  dft_str.apply(data + m);
396  }
397 };
398 
399 
400 // Take the next factor from the list
401 template<int_t N, int_t K, typename Tail, typename T, int S, class W1, int_t LastK>
402 class InFreq<N, Loki::Typelist<Pair<SInt<K>, SInt<0> >,Tail>, T, S, W1, LastK>
403 : public InFreq<N, Tail, T, S, W1, LastK> {};
404 
405 
406 // Specialization for prime N
407 template<int_t N, typename T, int S, class W1, int_t LastK>
408 class InFreq<N, Loki::Typelist<Pair<SInt<N>, SInt<1> >, Loki::NullType>,T,S,W1,LastK> {
409  DFTk_inp<N, 2, T, S> spec_inp;
410 public:
411  void apply(T* data)
412  {
413  spec_inp.apply(data);
414  }
415 };
416 
417 
418 /*
419 template<int_t N, typename NFact, typename T, int S, class W1, int_t LastK = 1>
420 class InFreqOOP;
421 
422 template<int_t N, typename Head, typename Tail, typename T, int S, class W1, int_t LastK>
423 class InFreqOOP<N, Loki::Typelist<Head,Tail>, T, S, W1, LastK>
424 {
425  typedef typename TempTypeTrait<T>::Result LocalVType;
426  static const int_t K = Head::first::value;
427  static const int_t M = N/K;
428  static const int_t M2 = M*2;
429  static const int_t N2 = N*2;
430  static const int_t LastK2 = LastK*2;
431 
432 // typedef typename Loki::TL::Next<RList,K-1>::Result RListK;
433  typedef typename IPowBig<W1,K>::Result WK;
434  typedef Loki::Typelist<Pair<typename Head::first, SInt<Head::second::value-1> >, Tail> NFactNext;
435  InFreqOOP<M,NFactNext,T,S,WK,K*LastK> dft_str;
436  T_DFTk_x_Im_OOP<K,M,T,S,W1,(N<=StaticLoopLimit)> dft_scaled;
437 
438 public:
439  void apply(const T* src, T* dst, T* buf)
440  {
441  dft_scaled.apply(src, buf);
442 
443  int_t lk = 0;
444  for (int_t m = 0; m < N2; m+=M2, lk+=LastK2)
445  dft_str.apply(buf + m, dst + lk, buf + m);
446  }
447 };
448 
449 
450 // Take the next factor from the list
451 template<int_t N, int_t K, typename Tail, typename T, int S, class W1, int_t LastK>
452 class InFreqOOP<N, Loki::Typelist<Pair<SInt<K>, SInt<0> >,Tail>, T, S, W1, LastK>
453 : public InFreqOOP<N, Tail, T, S, W1, LastK> {};
454 
455 
456 // Specialization for prime N
457 template<int_t N, typename T, int S, class W1, int_t LastK>
458 class InFreqOOP<N, Loki::Typelist<Pair<SInt<N>, SInt<1> >, Loki::NullType>,T,S,W1,LastK> {
459  DFTk<N, 2, LastK*2, T, S> spec;
460 public:
461  void apply(const T* src, T* dst, T*)
462  {
463  spec.apply(src, dst);
464  }
465 };
466 */
467 
468 } //namespace DFT
469 
470 #endif /*__gfftalgfreq_h*/
Generated on Mon Feb 10 2014 for Generative Fast Fourier Transforms (GFFT) by DoxyGen 1.8.3.1