【fftw】高速フーリエ変換の例

</p>
<!-- /wp:paragraph -->

<!-- wp:code -->
<pre class="wp-block-code"><code>
#include &lt;iostream>
#include &lt;fftw3.h>
#include &lt;cmath>
#include &lt;vector>
#include &lt;algorithm>

int main() {
 int N = 8;
 std::vector&lt;double> x = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0};
 std::vector&lt;double> y = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};

 fftw_complex *in, *out;
 in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
 out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);

 for (int i = 0; i &lt; N; i++) {
 in&#091;i]&#091;0] = y&#091;i];
 in&#091;i]&#091;1] = 0.0;
 }

 fftw_plan p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
 fftw_execute(p);

 std::vector&lt;double> power_spectrum(N);
 for (int i = 0; i &lt; N; i++) {
 power_spectrum&#091;i] = out&#091;i]&#091;0] * out&#091;i]&#091;0] + out&#091;i]&#091;1] * out&#091;i]&#091;1];
 }

 int n = 3;
 std::vector&lt;int> indices(N);
 for (int i = 0; i &lt; N; i++) {
 indices&#091;i] = i;
 }

 std::sort(indices.begin(), indices.end(), &#091;&amp;power_spectrum](int i1, int i2) {
 return power_spectrum&#091;i1] > power_spectrum&#091;i2];
 });

 std::vector&lt;bool> keep(N, false);
 for (int i = 0; i &lt; n; i++) {
 keep&#091;indices&#091;i]] = true;
 }

 for (int i = 0; i &lt; N; i++) {
 if (!keep&#091;i]) {
 out&#091;i]&#091;0] = 0.0;
 out&#091;i]&#091;1] = 0.0;
 }
 }

 fftw_plan p_inverse = fftw_plan_dft_1d(N, out, in, FFTW_BACKWARD, FFTW_ESTIMATE);
 fftw_execute(p_inverse);

 std::cout &lt;&lt; "Filtered signal (real part):" &lt;&lt; std::endl;
 for (int i = 0; i &lt; N; i++) {
 std::cout &lt;&lt; in&#091;i]&#091;0] / N &lt;&lt; std::endl;
 }

 fftw_destroy_plan(p);
 fftw_destroy_plan(p_inverse);
 fftw_free(in);
 fftw_free(out);

 return 0;
}
</code></pre>
<!-- /wp:code -->

<!-- wp:paragraph -->
<p>