LoginSignup
2
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@caa45040(caa45040)

STM32H743で最新アルゴリズム事典のFFTを動かした。(理系フリマで入手)

Last updated at Posted at 2024-08-20

x 理系フリマ831あるよ

参考

目的
FFTで遊ぶ

結果

o_coq336.jpg


      元のデータ    フーリエ変換  逆変換
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
2
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
2
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?