JUK1
dft.hpp
Go to the documentation of this file.
1 #ifndef _FFT_H_INC_
2 #define _FFT_H_INC_
3 #include <vector>
4 #include <complex>
5 #include <math.h>
6 #include <numbers>
7 
8 template<typename T>
9 std::complex<T>
10 nthRootOfUnity(T fraction) {
11  return std::exp(std::complex<T>(0, -2 * std::numbers::pi * fraction));
12 }
13 
14 template<typename T>
15 std::complex<T>
17  return std::exp(std::complex<T>(0, -2 * std::numbers::pi * numerator /
18  static_cast<double>(denominator)));
19 }
20 
21 
22 template<typename T>
23 std::vector<std::complex<T> >
24 dft(const std::vector<T> & inputData) {
25  const size_t length = inputData.size();
26  std::vector<std::complex<T> > toRet(length, std::complex<T>(0.0, 0.0));
27  for (size_t k = 0; k < length; k++) {
28  for (size_t n = 0; n < length; n++) {
29  toRet[k] += inputData[n] * nthRootOfUnity<T>(k * n, length);
30  }
31  }
32  return toRet;
33 }
34 
35 template<typename T>
36 std::vector<std::complex<T> >
37 idft(const std::vector<std::complex<T> > & inputData) {
38  std::vector<std::complex<T> > toRet;
39  const size_t length = inputData.size();
40  toRet.resize(length);
41  for (size_t n = 0; n < length; n++) {
42  for (size_t k = 0; k < length; k++) {
43  toRet[n] += inputData[k] *
44  nthRootOfUnity<T>(-static_cast<int>(k * n), length);
45  }
46  toRet[n] /= length;
47  }
48  return toRet;
49 }
50 
51 
52 template<typename T>
53 std::vector<std::complex<T> >
54 fft(const std::vector<T> & inputData) {
55  std::vector<std::complex<T> > result(inputData.size(), std::complex<T>(0, 0));
56  std::vector<std::complex<T> > scratch(inputData.size(), std::complex<T>(0, 0));
57  _fftHelperRadix2(inputData, result.begin(), scratch.begin(), 0, 0);
58  return result;
59 }
60 
61 template<typename T>
62 std::vector<std::complex<T> >
63 ifft(const std::vector<std::complex<T> > & inputData) {
64  std::vector<std::complex<T> > result(inputData.size(), std::complex<T>(0, 0));
65  std::vector<std::complex<T> > scratch(inputData.size(), std::complex<T>(0, 0));
66  _fftHelperRadix2<std::complex<T>, decltype(result.begin()), T,
67  -1>(inputData, result.begin(), scratch.begin(), 0, 0);
68  for (auto & num : result) {
69  num /= result.size();
70  }
71  return result;
72 }
73 
74 template<typename T, typename Iter, typename U = T, int dir = 1>
75 void
76 _fftHelperRadix2(const std::vector<T> & inputData, Iter result, Iter scratch,
77  size_t offset, size_t stride) {
78  const size_t len = (inputData.size() >> stride);
79  if (len > 2) {
80  _fftHelperRadix2<T, Iter, U, dir>(inputData, scratch, result, offset,
81  stride + 1);
82 
83  _fftHelperRadix2<T, Iter, U, dir>(inputData, scratch + (len / 2),
84  result + (len / 2), offset + (1 << stride),
85  stride + 1);
86 
87  for (size_t i = 0; i < len / 2; i++) {
88  result[i] = scratch[i] +
89  nthRootOfUnity<U>(dir * i, len) * scratch[i + len / 2];
90  result[i + len / 2] = scratch[i] - nthRootOfUnity<U>(dir * i, len) *
91  scratch[i + len / 2];
92  }
93  } else {
94  result[0] = inputData[offset] + inputData[offset + (1 << stride)];
95  result[1] = inputData[offset] - inputData[offset + (1 << stride)];
96  }
97 }
98 
99 #endif
pi
Definition: Freq.m:24
offset
Definition: Freq.m:2
std::vector< std::complex< T > > ifft(const std::vector< std::complex< T > > &inputData)
Definition: dft.hpp:63
void _fftHelperRadix2(const std::vector< T > &inputData, Iter result, Iter scratch, size_t offset, size_t stride)
Definition: dft.hpp:76
std::vector< std::complex< T > > fft(const std::vector< T > &inputData)
Definition: dft.hpp:54
std::vector< std::complex< T > > dft(const std::vector< T > &inputData)
Definition: dft.hpp:24
std::vector< std::complex< T > > idft(const std::vector< std::complex< T > > &inputData)
Definition: dft.hpp:37
std::complex< T > nthRootOfUnity(T fraction)
Definition: dft.hpp:10
for i
for k
Definition: interch.m:15
i **end numerator(col)
denominator(col)=0.0