Generative Fast Fourier Transforms (GFFT)  0.3
gfftstdspec.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 __gfftstdspec_h
16 #define __gfftstdspec_h
17 
22 #include "metafunc.h"
23 
24 namespace GFFT {
25 
26 using namespace MF;
27 
28 
30 
39 template<int_t N, int_t M, typename T, int S,
40 template<typename> class Complex>
41 class DFTk_inp<N,M,Complex<T>,S>
42 {
43  //GFFT_STATIC_ASSERT((N%2 == 1)) // N is assumed odd, otherwise compiler would not come here
44 
45  //typedef typename TempTypeTrait<T>::Result LocalVType;
46  static const int_t K = (N-1)/2;
47  static const int_t NM = N*M;
48 
49  T m_c[K], m_s[K];
50 
51 public:
52  DFTk_inp()
53  {
54  ComputeTwiddles<T, N, S, K>::apply(m_c, m_s);
55  }
56 
57  void apply(Complex<T>* data)
58  {
59  Complex<T> s[K], d[K];
60  for (int_t i=0; i<K; ++i) {
61  const int_t k = (i+1)*M;
62  s[i] = data[k] + data[NM-k];
63  d[i] = data[k] - data[NM-k];
64  }
65 
66  for (int_t i=1; i<K+1; ++i) {
67  Complex<T> t1(0,0), t2(0,0);
68  for (int_t j=0; j<K; ++j) {
69  const bool sign_change = (i*(j+1) % N) > K;
70  const int_t kk = (i+j*i)%N;
71  const int_t k = (kk>K) ? N-kk-1 : kk-1;
72  const T s1 = m_s[k]*d[j].imag();
73  const T s2 = m_s[k]*d[j].real();
74  t1 += m_c[k]*s[j];
75  Complex<T> tt(sign_change ? -s1 : s1, sign_change ? s2 : -s2);
76  t2 += tt;
77  }
78  const int_t k = i*M;
79  data[k] = data[0] + t1 + t2;
80  data[NM-k] = data[0] + t1 - t2;
81  }
82 
83  for (int_t i=0; i<K; ++i)
84  data[0] += s[i];
85  }
86 
87  void apply(Complex<T>* data, const Complex<T>* w)
88  {
89  Complex<T> s[K], d[K];
90  for (int_t i=0; i<K; ++i) {
91  const int_t k = (i+1)*M;
92  Complex<T> t1(data[k]*w[i]);
93  Complex<T> t2(data[NM-k]*w[N-i-2]);
94  s[i] = t1 + t2;
95  d[i] = t1 - t2;
96  }
97 
98  for (int_t i=1; i<K+1; ++i) {
99  Complex<T> t1(0,0), t2(0,0);
100  for (int_t j=0; j<K; ++j) {
101  const bool sign_change = (i*(j+1) % N) > K;
102  const int_t kk = (i+j*i)%N;
103  const int_t k = (kk>K) ? N-kk-1 : kk-1;
104  const T s1 = m_s[k]*d[j].imag();
105  const T s2 = m_s[k]*d[j].real();
106  t1 += m_c[k]*s[j];
107  Complex<T> tt(sign_change ? -s1 : s1, sign_change ? s2 : -s2);
108  t2 += tt;
109  }
110  const int_t k = i*M;
111  data[k] = data[0] + t1 + t2;
112  data[NM-k] = data[0] + t1 - t2;
113  }
114 
115  for (int_t i=0; i<K; ++i)
116  data[0] += s[i];
117  }
118 };
119 
120 template<int_t M, typename T, int S,
121 template<typename> class Complex>
122 class DFTk_inp<3,M,Complex<T>,S>
123 {
124  static const int_t I10 = M;
125  static const int_t I20 = M+M;
126 
127  T m_coef;
128 
129 public:
130  DFTk_inp() : m_coef(S * Sqrt<3, T>::value() * 0.5) { }
131 
132  void apply(Complex<T>* data)
133  {
134  Complex<T> sum(data[I10] + data[I20]);
135  Complex<T> dif(m_coef * (data[I10] - data[I20]));
136  Complex<T> t(data[0] - 0.5*sum);
137  data[0] += sum;
138  data[I10] = Complex<T>(t.real() + dif.imag(), t.imag() - dif.real());
139  data[I20] = Complex<T>(t.real() - dif.imag(), t.imag() + dif.real());
140  }
141  void apply(Complex<T>* data, const Complex<T>* w)
142  {
143  Complex<T> t1(data[I10]*w[0]);
144  Complex<T> t2(data[I20]*w[1]);
145 
146  Complex<T> sum(t1 + t2);
147  Complex<T> dif(m_coef * (t1 - t2));
148  Complex<T> t(data[0] - 0.5*sum);
149  data[0] += sum;
150  data[I10] = Complex<T>(t.real() + dif.imag(), t.imag() - dif.real());
151  data[I20] = Complex<T>(t.real() - dif.imag(), t.imag() + dif.real());
152  }
153 };
154 
155 template<int_t M, typename T, int S,
156 template<typename> class Complex>
157 class DFTk_inp<2,M,Complex<T>,S>
158 {
159 public:
160  void apply(Complex<T>* data)
161  {
162  Complex<T> t(data[M]);
163  data[M] = data[0] - t;
164  data[0] += t;
165  }
166  // For decimation-in-time
167  void apply(Complex<T>* data, const Complex<T>* w)
168  {
169  Complex<T> t(data[M] * (*w));
170  data[M] = data[0] - t;
171  data[0] += t;
172  }
173  // For decimation-in-frequency
174 // void apply(const Complex<T>* w, Complex<T>* data)
175 // {
176 // Complex<T> t(data[0] - data[M]);
177 // data[0] += data[M];
178 // data[M] = t * (*w);
179 // }
180 };
181 
183 
193 template<int_t N, int_t SI, int_t DI, typename T, int S,
194 template<typename> class Complex>
195 class DFTk<N,SI,DI,Complex<T>,S>
196 {
197  //GFFT_STATIC_ASSERT(N%2 == 1) // N is assumed odd, otherwise compiler would not come here
198 
199  //typedef typename TempTypeTrait<T>::Result LocalVType;
200  static const int_t K = (N-1)/2;
201  static const int_t NSI = N*SI;
202  static const int_t NDI = N*DI;
203 
204  T m_c[K], m_s[K];
205 
206 public:
207  DFTk()
208  {
209  ComputeTwiddles<T, N, S, K>::apply(m_c, m_s);
210  }
211 
212  void apply(const Complex<T>* src, Complex<T>* dst)
213  {
214  Complex<T> s[K], d[K];
215  for (int_t i=0; i<K; ++i) {
216  const int_t k = (i+1)*SI;
217  s[i] = src[k] + src[NSI-k];
218  d[i] = src[k] - src[NSI-k];
219  }
220 
221  for (int_t i=1; i<K+1; ++i) {
222  Complex<T> t1(0,0), t2(0,0);
223  for (int_t j=0; j<K; ++j) {
224  const bool sign_change = (i*(j+1) % N) > K;
225  const int_t kk = (i+j*i)%N;
226  const int_t k = (kk>K) ? N-kk-1 : kk-1;
227  const T s1 = m_s[k]*d[j].imag();
228  const T s2 = m_s[k]*d[j].real();
229  t1 += m_c[k]*s[j];
230  Complex<T> tt(sign_change ? -s1 : s1, sign_change ? s2 : -s2);
231  t2 += tt;
232  }
233  const int_t k = i*DI;
234  dst[k] = src[0] + t1 + t2;
235  dst[NDI-k] = src[0] + t1 - t2;
236  }
237 
238  dst[0] = src[0];
239  for (int_t i=0; i<K; ++i)
240  dst[0] += s[i];
241  }
242 };
243 
244 template<int_t SI, int_t DI, typename T, int S,
245 template<typename> class Complex>
246 class DFTk<3,SI,DI,Complex<T>,S>
247 {
248  static const int_t SI2 = SI+SI;
249  static const int_t DI2 = DI+DI;
250  T m_coef;
251 
252 public:
253  DFTk() : m_coef(S * Sqrt<3, T>::value() * 0.5) { } // sqrt(3)/2 = sin(2pi/3)
254 
255  void apply(const Complex<T>* src, Complex<T>* dst)
256  {
257  // 4 mult, 12 add
258  Complex<T> s(src[SI] + src[SI2]);
259  Complex<T> d(m_coef * (src[SI] - src[SI2]));
260  Complex<T> t(src[0] - 0.5*s);
261  dst[0] = src[0] + s;
262  dst[DI] = Complex<T>(t.real() + d.imag(), t.imag() - d.real());
263  dst[DI2] = Complex<T>(t.real() - d.imag(), t.imag() + d.real());
264  }
265 };
266 
267 template<int_t SI, int_t DI, typename T, int S,
268 template<typename> class Complex>
269 class DFTk<2,SI,DI,Complex<T>,S>
270 {
271 public:
272  void apply(const Complex<T>* src, Complex<T>* dst)
273  {
274  // the temporaries tr, ti are necessary, because may happen src == dst
275  Complex<T> t(src[0] - src[SI]);
276  dst[0] = src[0] + src[SI];
277  dst[DI] = t;
278  }
279 };
280 
285 template<typename T,
286 template<typename> class Complex>
287 inline void _spec2(Complex<T>* data) {
288  Complex<T> t(data[1]);
289  data[1] = data[0]-t;
290  data[0] += t;
291 }
292 
293 template<typename T,
294 template<typename> class Complex>
295 inline void _spec2(const Complex<T>* src, Complex<T>* dst)
296 {
297  dst[0] = src[0] + src[1];
298  dst[1] = src[0] - src[1];
299 }
300 
301 
302 } //namespace DFT
303 
304 #endif /*__gfftspec_h*/
Generated on Mon Feb 10 2014 for Generative Fast Fourier Transforms (GFFT) by DoxyGen 1.8.3.1