x 理系フリマ831あるよ

参考

目的
FFTで遊ぶ

結果


      元のデータ    フーリエ変換  逆変換
0|6.00 0.00|0.00 0.00|6.00 0.00
1|8.83 0.00|-0.00 0.00|8.83 -0.00
2|8.91 0.00|-0.00 0.00|8.91 0.00
3|5.69 0.00|3.00 -0.00|5.69 -0.00
4|0.77 0.00|0.00 0.00|0.77 0.00
5|-3.24 0.00|-0.00 0.00|-3.24 0.00
6|-4.50 0.00|0.00 -0.00|-4.50 0.00
7|-3.22 0.00|0.00 0.00|-3.22 -0.00
8|-1.41 0.00|-0.00 0.00|-1.41 0.00
9|-1.31 0.00|-0.00 -2.00|-1.31 -0.00
10|-3.66 0.00|0.00 -0.00|-3.66 -0.00
11|-7.13 0.00|0.00 -0.00|-7.13 -0.00
12|-9.24 0.00|-0.00 -0.00|-9.24 -0.00
13|-8.17 0.00|-0.00 -0.00|-8.17 0.00
14|-4.11 0.00|0.00 -0.00|-4.11 0.00
15|0.80 0.00|0.00 -0.00|0.80 0.00
16|4.00 0.00|0.00 0.00|4.00 0.00
17|4.28 0.00|0.00 -0.00|4.28 0.00
18|2.55 0.00|0.00 0.00|2.55 0.00
19|1.11 0.00|-0.00 0.00|1.11 0.00
20|1.85 0.00|-0.00 0.00|1.85 0.00
21|4.81 0.00|-0.00 -0.00|4.81 0.00
22|8.11 0.00|-0.00 -0.00|8.11 -0.00
23|9.27 0.00|-0.00 -0.00|9.27 -0.00
24|7.07 0.00|0.00 -0.00|7.07 0.00
25|2.44 0.00|0.00 -0.00|2.44 -0.00
26|-2.16 0.00|0.00 -0.00|-2.16 0 0.00|0.77 0.00
5|-3.24 0.00|-0.00 0.00|-3.24 0.00
6|-4.50 0.00|0.00 -0.00|-4.50 0.00
7|-3.22 0.00|0.00 0.00|-3.22 -0.00
8|-1.41 0.00|-0.00 0.00|-1.41 0.00
9|-1.31 0.00|-0.00 -2.00|-1.31 -0.00
10|-3.66 0.00|0.00 -0.00|-3.66 -0.00
11|-7.13 0.00|0.00 -0.00|-7.13 -0.00
12|-9.24 0.00|-0.00 -0.00|-9.24 -0.00
13|-8.17 0.00|-0.00 -0.00|-8.17 0.00
14|-4.11 0.00|0.00 -0.00|-4.11 0.00
15|0.80 0.00|0.00 -0.00|0.80 0.00
16|4.00 0.00|0.00 0.00|4.00 0.00
17|4.28 0.00|0.00 -0.00|4.28 0.00
18|2.55 0.00|0.00 0.00|2.55 0.00
19|1.11 0.00|-0.00 0.00|1.11 0.00
20|1.85 0.00|-0.00 0.00|1.85 0.00
21|4.81 0.00|-0.00 -0.00|4.81 0.00
22|8.11 0.00|-0.00 -0.00|8.11 -0.00
23|9.27 0.00|-0.00 -0.00|9.27 -0.00
24|7.07 0.00|0.00 -0.00|7.07 0.00
25|2.44 0.00|0.00 -0.00|2.44 -0.00
26|-2.16 0.00|0.00 -0.00|-2.16 -0.00
27|-4.42 0.00|0.00 0.00|-4.42 -0.00
28|-3.83 0.00|0.00 0.00|-3.83 0.00
29|-1.92 0.00|-0.00 0.00|-1.92 -0.00
30|-1.07 0.00|-0.00 -0.00|-1.07 -0.00
31|-2.65 0.00|0.00 0.00|-2.65 -0.00
32|-6.00 0.00|-0.00 0.00|-6.00 0.00
33|-8.83 0.00|0.00 -0.00|-8.83 0.00
34|-8.91 0.00|-0.00 0.00|-8.91 0.00
35|-5.69 0.00|-0.00 -0.00|-5.69 0.00
36|-0.77 0.00|0.00 -0.00|-0.77 -0.00
37|3.24 0.00|0.00 -0.00|3.24 0.00
38|4.50 0.00|0.00 0.00|4.50 -0.00
39|3.22 0.00|0.00 0.00|3.22 0.00
40|1.41 0.00|0.00 0.00|1.41 0.00
41|1.31 0.00|-0.00 0.00|1.31 0.00
42|3.66 0.00|-0.00 0.00|3.66 0.00
43|7.13 0.00|-0.00 0.00|7.13 0.00
44|9.24 0.00|-0.00 -0.00|9.24 0.00
45|8.17 0.00|-0.00 -0.00|8.17 -0.00
46|4.11 0.00|0.00 -0.00|4.11 -0.00
47|-0.80 0.00|0.00 0.00|-0.80 -0.00
48|-4.00 0.00|0.00 -0.00|-4.00 0.00
49|-4.28 0.00|0.00 0.00|-4.28 -0.00
50|-2.55 0.00|0.00 0.00|-2.55 -0.00
51|-1.11 0.00|-0.00 0.00|-1.11 -0.00
52|-1.85 0.00|-0.00 0.00|-1.85 0.00
53|-4.81 0.00|0.00 0.00|-4.81 -0.00
54|-8.11 0.00|0.00 0.00|-8.11 -0.00
55|-9.27 0.00|-0.00 2.00|-9.27 -0.00
56|-7.07 0.00|-0.00 -0.00|-7.07 0.00
57|-2.44 0.00|0.00 -0.00|-2.44 0.00
58|2.16 0.00|0.00 0.00|2.16 0.00
59|4.42 0.00|-0.00 -0.00|4.42 0.00
60|3.83 0.00|0.00 -0.00|3.83 -0.00
61|1.92 0.00|3.00 0.00|1.92 0.00
62|1.07 0.00|-0.00 -0.00|1.07 0.00
63|2.65 0.00|-0.00 -0.00|2.65 0.00

プログラム




//fft_c_H743_1



/***********************************************************
    fft.c -- FFT (高速Fourier変換)
***********************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define PI 3.14159265358979323846
/*
  関数{\tt fft()}の下請けとして三角関数表を作る.
*/
static void make_sintbl(int n, double sintbl[])
{
    int i, n2, n4, n8;
    double c, s, dc, ds, t;

    n2 = n / 2;  n4 = n / 4;  n8 = n / 8;
    t = sin(PI / n);
    dc = 2 * t * t;  ds = sqrt(dc * (2 - dc));
    t = 2 * dc;  c = sintbl[n4] = 1;  s = sintbl[0] = 0;
    for (i = 1; i < n8; i++) {
        c -= dc;  dc += t * c;
        s += ds;  ds -= t * s;
        sintbl[i] = s;  sintbl[n4 - i] = c;
    }
    if (n8 != 0) sintbl[n8] = sqrt(0.5);
    for (i = 0; i < n4; i++) {
        sintbl[n2 - i] = sintbl[i];
        sintbl[i + n2] = - sintbl[i];
    }
}
/*
  関数{\tt fft()}の下請けとしてビット反転表を作る.
*/
static void make_bitrev(int n, int bitrev[])
{
    int i, j, k, n2;

    n2 = n / 2;  i = j = 0;
    for ( ; ; ) {
        bitrev[i] = j;
        if (++i >= n) break;
        k = n2;
        while (k <= j) {  j -= k;  k /= 2;  }
        j += k;
    }
}
/*
  高速Fourier変換 (Cooley--Tukeyのアルゴリズム).
  標本点の数 {\tt n} は2の整数乗に限る.
  {\tt x[$k$]} が実部, {\tt y[$k$]} が虚部 ($k = 0$, $1$, $2$,
  \ldots, $|{\tt n}| - 1$).
  結果は {\tt x[]}, {\tt y[]} に上書きされる.
  ${\tt n} = 0$ なら表のメモリを解放する.
  ${\tt n} < 0$ なら逆変換を行う.
  前回と異なる $|{\tt n}|$ の値で呼び出すと,
  三角関数とビット反転の表を作るために多少余分に時間がかかる.
  この表のための記憶領域獲得に失敗すると1を返す (正常終了時
  の戻り値は0).
  これらの表の記憶領域を解放するには ${\tt n} = 0$ として
  呼び出す (このときは {\tt x[]}, {\tt y[]} の値は変わらない).
*/
int fft(int n, double x[], double y[])
{
    static int    last_n = 0;    /* 前回呼出し時の {\tt n} */
    static int    *bitrev = NULL; /* ビット反転表 */
    static double *sintbl = NULL; /* 三角関数表 */
    int i, j, k, ik, h, d, k2, n4, inverse;
    double t, s, c, dx, dy;

    /* 準備 */
    if (n < 0) {
        n = -n;  inverse = 1;  /* 逆変換 */
    } else inverse = 0;
    n4 = n / 4;
    if (n != last_n || n == 0) {
        last_n = n;
        if (sintbl != NULL) free(sintbl);
        if (bitrev != NULL) free(bitrev);
        if (n == 0) return 0;  /* 記憶領域を解放した */
        //sintbl = malloc((n - n4) * sizeof(double));
        sintbl = (double *)malloc(10000);
        //bitrev = malloc(n * sizeof(int));
        bitrev = (int    *)malloc(10000);
        if (sintbl == NULL || bitrev == NULL) {
            fprintf(stderr, "記憶領域不足\n");  return 1;
        }
        make_sintbl(n, sintbl);
        make_bitrev(n, bitrev);
    }
    for (i = 0; i < n; i++) {    /* ビット反転 */
        j = bitrev[i];
        if (i < j) {
            t = x[i];  x[i] = x[j];  x[j] = t;
            t = y[i];  y[i] = y[j];  y[j] = t;
        }
    }
    for (k = 1; k < n; k = k2) {    /* 変換 */
        h = 0;  k2 = k + k;  d = n / k2;
        for (j = 0; j < k; j++) {
            c = sintbl[h + n4];
            if (inverse) s = - sintbl[h];
            else         s =   sintbl[h];
            for (i = j; i < n; i += k2) {
                ik = i + k;
                dx = s * y[ik] + c * x[ik];
                dy = c * y[ik] - s * x[ik];
                x[ik] = x[i] - dx;  x[i] += dx;
                y[ik] = y[i] - dy;  y[i] += dy;
            }
            h += d;
        }
    }
    if (! inverse)    /* 逆変換でないならnで割る */
        for (i = 0; i < n; i++) {  x[i] /= n;  y[i] /= n;  }
    return 0;  /* 正常終了 */
}

#define N 64





// the setup function runs once when you press reset or power the board
void setup() {

  // initialize serial communication at 9600 bits per second:
  Serial.begin(9600);
  Serial.println();
  Serial.println("START");
  Serial.println();

  // initialize digital pin LED_BUILTIN as an output.
  pinMode(LED_BUILTIN, OUTPUT);

}

// the loop function runs over and over again forever
void loop() {

    int i;
    static double x1[N], y1[N], x2[N], y2[N], x3[N], y3[N];

    for (i = 0; i < N; i++) {
        x1[i] = x2[i] = 6 * cos( 6 * PI * i / N)
                      + 4 * sin(18 * PI * i / N);
        y1[i] = y2[i] = 0;
    }
    if (fft(N, x2, y2)){
    //return 1;
    return;
    }
    for (i = 0; i < N; i++) {
        x3[i] = x2[i];  y3[i] = y2[i];
    }
    if (fft(-N, x3, y3)){
    //return 1;
    return;
    }
    printf("      元のデータ    フーリエ変換  逆変換\n");
    //for (i = 0; i < N; i++)
    //    printf("%4d | %6.3f %6.3f | %6.3f %6.3f | %6.3f %6.3f\n",
    //        i, x1[i], y1[i], x2[i], y2[i], x3[i], y3[i]);
    //return 0;

    for (i = 0; i < N; i++){
        Serial.print(i); 
        Serial.print("|");
        Serial.print(x1[i]); Serial.print(' ');
        Serial.print(y1[i]);
        Serial.print("|");
        Serial.print(x2[i]); Serial.print(' ');
        Serial.print(y2[i]);
        Serial.print("|");
        Serial.print(x3[i]); Serial.print(' ');
        Serial.print(y3[i]);
        Serial.println();
    }


  //while(1){} //debug

  
  digitalWrite(LED_BUILTIN, HIGH);   // turn the LED on (HIGH is the voltage level)
  delay(1000);                       // wait for a second
  digitalWrite(LED_BUILTIN, LOW);    // turn the LED off by making the voltage LOW
  delay(1000);                       // wait for a second
}

おまけ　オンラインコンパイラで実行


      元のデータ    フーリエ変換  逆変換
   0 |  6.000  0.000 |  0.000  0.000 |  6.000  0.000
   1 |  8.834  0.000 | -0.000  0.000 |  8.834 -0.000
   2 |  8.912  0.000 | -0.000  0.000 |  8.912  0.000
   3 |  5.692  0.000 |  3.000 -0.000 |  5.692 -0.000
   4 |  0.765  0.000 |  0.000  0.000 |  0.765  0.000
   5 | -3.240  0.000 | -0.000  0.000 | -3.240  0.000
   6 | -4.496  0.000 |  0.000 -0.000 | -4.496  0.000
   7 | -3.220  0.000 |  0.000  0.000 | -3.220 -0.000
   8 | -1.414  0.000 | -0.000  0.000 | -1.414  0.000
   9 | -1.311  0.000 | -0.000 -2.000 | -1.311 -0.000
  10 | -3.662  0.000 |  0.000 -0.000 | -3.662 -0.000
  11 | -7.132  0.000 |  0.000 -0.000 | -7.132 -0.000
  12 | -9.239  0.000 | -0.000 -0.000 | -9.239 -0.000
  13 | -8.166  0.000 | -0.000 -0.000 | -8.166  0.000
  14 | -4.114  0.000 |  0.000 -0.000 | -4.114  0.000
  15 |  0.796  0.000 |  0.000 -0.000 |  0.796  0.000
  16 |  4.000  0.000 |  0.000  0.000 |  4.000  0.000
  17 |  4.279  0.000 |  0.000 -0.000 |  4.279  0.000
  18 |  2.553  0.000 |  0.000  0.000 |  2.553  0.000
  19 |  1.110  0.000 | -0.000  0.000 |  1.110  0.000
  20 |  1.848  0.000 | -0.000  0.000 |  1.848  0.000
  21 |  4.810  0.000 | -0.000 -0.000 |  4.810  0.000
  22 |  8.107  0.000 | -0.000 -0.000 |  8.107  0.000
  23 |  9.272  0.000 | -0.000 -0.000 |  9.272 -0.000
  24 |  7.071  0.000 |  0.000 -0.000 |  7.071  0.000
  25 |  2.436  0.000 |  0.000 -0.000 |  2.436 -0.000
  26 | -2.155  0.000 |  0.000 -0.000 | -2.155 -0.000
  27 | -4.416  0.000 |  0.000  0.000 | -4.416 -0.000
  28 | -3.827  0.000 |  0.000  0.000 | -3.827  0.000
  29 | -1.921  0.000 | -0.000  0.000 | -1.921 -0.000
  30 | -1.066  0.000 | -0.000 -0.000 | -1.066 -0.000
  31 | -2.650  0.000 |  0.000  0.000 | -2.650 -0.000
  32 | -6.000  0.000 | -0.000  0.000 | -6.000  0.000
  33 | -8.834  0.000 |  0.000 -0.000 | -8.834  0.000
  34 | -8.912  0.000 | -0.000  0.000 | -8.912  0.000
  35 | -5.692  0.000 | -0.000 -0.000 | -5.692  0.000
  36 | -0.765  0.000 |  0.000 -0.000 | -0.765 -0.000
  37 |  3.240  0.000 |  0.000 -0.000 |  3.240  0.000
  38 |  4.496  0.000 |  0.000  0.000 |  4.496 -0.000
  39 |  3.220  0.000 |  0.000  0.000 |  3.220  0.000
  40 |  1.414  0.000 |  0.000  0.000 |  1.414  0.000
  41 |  1.311  0.000 | -0.000  0.000 |  1.311  0.000
  42 |  3.662  0.000 | -0.000  0.000 |  3.662  0.000
  43 |  7.132  0.000 | -0.000  0.000 |  7.132  0.000
  44 |  9.239  0.000 | -0.000 -0.000 |  9.239  0.000
  45 |  8.166  0.000 | -0.000 -0.000 |  8.166 -0.000
  46 |  4.114  0.000 |  0.000 -0.000 |  4.114 -0.000
  47 | -0.796  0.000 |  0.000  0.000 | -0.796 -0.000
  48 | -4.000  0.000 |  0.000 -0.000 | -4.000  0.000
  49 | -4.279  0.000 |  0.000  0.000 | -4.279 -0.000
  50 | -2.553  0.000 |  0.000  0.000 | -2.553 -0.000
  51 | -1.110  0.000 | -0.000  0.000 | -1.110 -0.000
  52 | -1.848  0.000 | -0.000  0.000 | -1.848  0.000
  53 | -4.810  0.000 |  0.000  0.000 | -4.810 -0.000
  54 | -8.107  0.000 |  0.000  0.000 | -8.107  0.000
  55 | -9.272  0.000 | -0.000  2.000 | -9.272 -0.000
  56 | -7.071  0.000 | -0.000 -0.000 | -7.071  0.000
  57 | -2.436  0.000 |  0.000 -0.000 | -2.436  0.000
  58 |  2.155  0.000 |  0.000  0.000 |  2.155  0.000
  59 |  4.416  0.000 | -0.000 -0.000 |  4.416  0.000
  60 |  3.827  0.000 |  0.000 -0.000 |  3.827 -0.000
  61 |  1.921  0.000 |  3.000  0.000 |  1.921  0.000
  62 |  1.066  0.000 | -0.000 -0.000 |  1.066  0.000
  63 |  2.650  0.000 | -0.000 -0.000 |  2.650  0.000