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.  
  86.     return FFTRadix2(x, W);
  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. OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!
 
Top