• API
• FAQ
• Tools
• Archive
SHARE
TWEET

# Untitled

a guest Aug 17th, 2019 65 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. #include <cassert>
2. #include <complex>
3. #include <vector>
4.
5. // To demonstrate runtime
6. #include <chrono>
7. #include <iostream>
8.
9. static std::vector<std::complex<double>> FFTRadix2(const std::vector<std::complex<double>>& x, const std::vector<std::complex<double>>& W);
10. static bool IsPowerOf2(size_t x);
11. static size_t ReverseBits(const size_t x, const size_t n);
12.
13. std::vector<std::complex<double>> FFT(const std::vector<double>& x)
14. {
15.     size_t N = x.size();
16.
17.     // Radix2 FFT requires length of the input signal to be a power of 2
18.     // TODO: Implement other algorithms for when N is not a power of 2
19.     assert(IsPowerOf2(N));
20.
21.     // Taking advantage of symmetry the FFT of a real signal can be computed
22.     // using a single N/2-point complex FFT. Split the input signal into its
23.     // even and odd components and load the data into a single complex vector.
24.     std::vector<std::complex<double>> x_p(N / 2);
25.     for (size_t n = 0; n < N / 2; ++n)
26.     {
27.         // x_p[n] = x[2n] + jx[2n + 1]
28.         x_p[n] = std::complex<double>(x[2 * n], x[2 * n + 1]);
29.     }
30.
31.     // Pre-calculate twiddle factors
32.     std::vector<std::complex<double>> W(N / 2);
33.     std::vector<std::complex<double>> W_p(N / 4);
34.     for (size_t k = 0; k < N / 2; ++k)
35.     {
36.         W[k] = std::polar(1.0, -2 * M_PI * k / N);
37.
38.         // The N/2-point complex DFT uses only the even twiddle factors
39.         if (k % 2 == 0)
40.         {
41.             W_p[k / 2] = W[k];
42.         }
43.     }
44.
45.     // Perform the N/2-point complex FFT
46.     std::vector<std::complex<double>> X_p = FFTRadix2(x_p, W_p);
47.
48.     // Extract the N-point FFT of the real signal from the results
49.     std::vector<std::complex<double>> X(N);
50.     X[0] = X_p[0].real() + X_p[0].imag();
51.     for (size_t k = 1; k < N / 2; ++k)
52.     {
53.         // Extract the FFT of the even components
54.         auto A = std::complex<double>(
55.             (X_p[k].real() + X_p[N / 2 - k].real()) / 2,
56.             (X_p[k].imag() - X_p[N / 2 - k].imag()) / 2);
57.
58.         // Extract the FFT of the odd components
59.         auto B = std::complex<double>(
60.             (X_p[N / 2 - k].imag() + X_p[k].imag()) / 2,
61.             (X_p[N / 2 - k].real() - X_p[k].real()) / 2);
62.
63.         // Sum the results and take advantage of symmetry
64.         X[k] = A + W[k] * B;
65.         X[k + N / 2] = A - W[k] * B;
66.     }
67.
68.     return X;
69. }
70.
71. std::vector<std::complex<double>> FFT(const std::vector<std::complex<double>>& x)
72. {
73.     size_t N = x.size();
74.
75.     // Radix2 FFT requires length of the input signal to be a power of 2
76.     // TODO: Implement other algorithms for when N is not a power of 2
77.     assert(IsPowerOf2(N));
78.
79.     // Pre-calculate twiddle factors
80.     std::vector<std::complex<double>> W(N / 2);
81.     for (size_t k = 0; k < N / 2; ++k)
82.     {
83.         W[k] = std::polar(1.0, -2 * M_PI * k / N);
84.     }
85.
87. }
88.
89. static std::vector<std::complex<double>> FFTRadix2(const std::vector<std::complex<double>>& x, const std::vector<std::complex<double>>& W)
90. {
91.     size_t N = x.size();
92.
93.     // Radix2 FFT requires length of the input signal to be a power of 2
94.     assert(IsPowerOf2(N));
95.
96.     // Calculate how many stages the FFT must compute
97.     size_t stages = static_cast<size_t>(log2(N));
98.
99.     // Pre-load the output vector with the input data using bit-reversed indexes
100.     std::vector<std::complex<double>> X(N);
101.     for (size_t n = 0; n < N; ++n)
102.     {
103.         X[n] = x[ReverseBits(n, stages)];
104.     }
105.
106.     // Calculate the FFT one stage at a time and sum the results
107.     for (size_t stage = 1; stage <= stages; ++stage)
108.     {
109.         size_t N_stage = static_cast<size_t>(std::pow(2, stage));
110.         size_t W_offset = static_cast<size_t>(std::pow(2, stages - stage));
111.         for (size_t k = 0; k < N; k += N_stage)
112.         {
113.             for (size_t n = 0; n < N_stage / 2; ++n)
114.             {
115.                 auto tmp = X[k + n];
116.                 X[k + n] = tmp + W[n * W_offset] * X[k + n + N_stage / 2];
117.                 X[k + n + N_stage / 2] = tmp - W[n * W_offset] * X[k + n + N_stage / 2];
118.             }
119.         }
120.     }
121.
122.     return X;
123. }
124.
125. // Returns true if x is a power of 2
126. static bool IsPowerOf2(size_t x)
127. {
128.     return x && (!(x & (x - 1)));
129. }
130.
131. // Given x composed of n bits, returns x with the bits reversed
132. static size_t ReverseBits(const size_t x, const size_t n)
133. {
134.     size_t xReversed = 0;
135.     for (size_t i = 0; i < n; ++i)
136.     {
137.         xReversed = (xReversed << 1U) | ((x >> i) & 1U);
138.     }
139.
140.     return xReversed;
141. }
142.
143. int main()
144. {
145.     size_t N = 16777216;
146.     std::vector<double> x(N);
147.
148.     int f_s = 8000;
149.     double t_s = 1.0 / f_s;
150.
151.     for (size_t n = 0; n < N; ++n)
152.     {
153.         x[n] = std::sin(2 * M_PI * 1000 * n * t_s)
154.             + 0.5 * std::sin(2 * M_PI * 2000 * n * t_s + 3 * M_PI / 4);
155.     }
156.
157.     auto start = std::chrono::high_resolution_clock::now();
158.     auto X = FFT(x);
159.     auto stop = std::chrono::high_resolution_clock::now();
160.     auto duration = std::chrono::duration_cast<std::chrono::microseconds>(stop - start);
161.
162.     std::cout << duration.count() << std::endl;
163. }
164.
165. import numpy as np
166. import datetime
167.
168. N = 16777216
169. f_s = 8000.0
170. t_s = 1/f_s
171.
172. t = np.arange(0, N*t_s, t_s)
173. y = np.sin(2*np.pi*1000*t) + 0.5*np.sin(2*np.pi*2000*t + 3*np.pi/4)
174.
175. start = datetime.datetime.now()
176. Y = np.fft.fft(y)
177. stop = datetime.datetime.now()
178. duration = stop - start
179.
180. print(duration.total_seconds()*1e6)
181.
182. for (size_t n = 0, rev = 0; n < N; ++n)
183. {
184.     X[n] = x[rev];
185.     size_t change = n ^ (n + 1);
186.     rev ^= change << (__lzcnt64(change) - (64 - stages));
187. }
188.
189. size_t change = n ^ (n + 1);
190. std::bitset<64> bits(~change);
191. rev ^= change << (bits.count() - (64 - stages));
192.
193. std::vector<std::complex<double>> a(10);
194. double *b = reinterpret_cast<double *>(a.data());
195.
196. span<std::complex<double>> x_p(reinterpret_cast<std::complex<double>*>(x.data()), x.size() / 2);
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.
Not a member of Pastebin yet?