EasyFFT：Arduinoプロジェクト向け信頼性の高いFFTライブラリ

必要なコンポーネントと用品

パーツ	数量
Arduino Nano R3	1

アプリとオンラインサービス

Arduino IDE

プロジェクト概要

取得した信号から周波数を測定するのは、特に計算能力の低いArduinoでは難しいタスクです。ゼロクロス法（信号がゼロラインを横切る回数を数える）は、複数の周波数が混在する信号では機能しません。

このプロジェクトは、FFT（高速フーリエ変換）をArduinoで簡単に実装できるコードを提供することを目的としています。音楽、信号解析などのプロジェクトで役立ちます。FFTの数学的背景については説明しませんが、参考資料を記載しています。

コードの適用のみ知りたい場合は、ステップ3までスキップしてください。

より高速（約3倍）で精度をわずか（約5%）妥協できる場合は、拙作の ApproxFFT を参照してください。

ステップ1：高速フーリエ変換（FFT）

DFT（離散フーリエ変換）の計算を高速化するため、James CooleyとJohn TukeyによってFFTアルゴリズムが開発されました。これは20世紀で最も重要なアルゴリズムの一つとされています。信号を奇数・偶数列に分割することで計算量を削減し、通常のDFTがN²回の複素数乗算を必要とするのに対し、FFTではN log N回で済みます。サンプル数が多い場合に大きな利点があります。

数学的詳細については以下を参照してください：

ステップ2：コードの解説

高速なサイン/コサイン計算

FFTではサイン/コサインを繰り返し計算しますが、Arduino標準の関数は遅く、64サンプルで実行時間が倍増します。これを解決するため、0〜90度のサイン値を255倍した整数（バイト）として配列 sine_data[91] に格納し、グローバル変数として宣言しています。これによりメモリ使用量を1/4に削減し、高速化しています。

また、FFT実行後に上位5つのピーク周波数を降順で格納するグローバル配列 f_peaks[5] も宣言しています。

byte sine_data[91] = { 
  0, 4, 9, 13, 18, 22, 27, 31, 35, 40, 44, 49, 53, 57, 62, 66, 70, 75, 79, 83, 87, 
  91, 96, 100, 104, 108, 112, 116, 120, 124, 127, 131, 135, 139, 143, 146, 150, 153, 
  157, 160, 164, 167, 171, 174, 177, 180, 183, 186, 189, 192, 195, 198, 201, 204, 206, 
  209, 211, 214, 216, 219, 221, 223, 225, 227, 229, 231, 233, 235, 236, 238, 240, 241, 
  243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 253, 254, 254, 254, 255, 255, 
  255, 255 
};
float f_peaks[5];

このテーブルを用いて、任意の角度のサイン/コサインを高速に求める関数 sine() と cosine() を実装しています。精度は多少落ちますが、64ポイントで約8ms、128ポイントで約20msの高速化が得られます。

ステップ3：FFT関数の解説

FFTはサンプル数が2の累乗（2, 4, 8, 16, 32, 64, 128…）である必要があります。指定数が2の累乗でない場合、直近の低い2の累乗に切り捨てられます（例：70 → 64）。Arduino Nanoではメモリ制約により128サンプル以上は困難です（64推奨）。

処理の流れ：

ビットリバース順序の生成
入力データの並べ替え
FFT実行
複素数の振幅計算
ピーク検出と降順ソート
結果を f_peaks[] から取得

（ピーク以外のデータが必要な場合は、ローカル変数をグローバルにコピーするようコードを修正してください。）

ステップ4：コードのテスト

サンプリング周波数10Hz、信号周波数1.25Hzの三角波を入力としてテストしました。Scilabで計算したFFT結果と概ね一致しますが、高速化のための近似サイン使用により完全一致ではありません。出力される周波数ビン（1.25Hz、3.75Hzなど）は、真の周波数が含まれる範囲を示します。

Arduino Nanoでの実行時間：

16ポイント：約4ms
32ポイント：約10ms
64ポイント：約26ms
128ポイント：約53ms

ステップ5：結論

このFFTコードはリアルタイム用途に十分な速度（約30ms）で動作します。分解能はサンプル数に依存し、Arduinoのメモリが制限要因です。Arduino Megaやより高性能なボードを使えば精度向上が可能です。

質問、提案、修正があればコメントでお知らせください。

コード全文

/* 
// Example data:
int data[64] = {14, 30, 35, 34, 34, 40, 46, 45, 30, 4, -26, -48, -55, -49, -37, -28, -24, -22, -13, 6, 32, 55, 65, 57, 38, 17, 1, -6, -11, -19, -34, -51, -61, -56, -35, -7, 18, 32, 35, 34, 35, 41, 46, 43, 26, -2, -31, -50, -55, -47, -35, -27, -24, -21, -10, 11, 37, 58, 64, 55, 34, 13, -1, -7};
*/
//---------------------------------------------------------------------------//
// グローバル変数として宣言必須
byte sine_data[91] = { 
  0, 4, 9, 13, 18, 22, 27, 31, 35, 40, 44, 49, 53, 57, 62, 66, 70, 75, 79, 83, 87, 
  91, 96, 100, 104, 108, 112, 116, 120, 124, 127, 131, 135, 139, 143, 146, 150, 153, 
  157, 160, 164, 167, 171, 174, 177, 180, 183, 186, 189, 192, 195, 198, 201, 204, 206, 
  209, 211, 214, 216, 219, 221, 223, 225, 227, 229, 231, 233, 235, 236, 238, 240, 241, 
  243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 253, 254, 254, 254, 255, 255, 
  255, 255 
};
float f_peaks[5]; // 上位5つのピーク周波数（降順）
//---------------------------------------------------------------------------//

void setup() {
  Serial.begin(250000);
}

void loop() {
  /* 使用例:
  FFT(data, 64, 100); // 64サンプル、100Hzサンプリング
  Serial.println(f_peaks[0]);
  Serial.println(f_peaks[1]);
  delay(99999);
  */
}

//----------------------------- FFT Function ----------------------------------------------//
float FFT(int in[], int N, float Frequency) {
  /* ArduinoでFFTを実行する関数
   * セットアップ: sine_data[91]をプログラム冒頭にペースト、この関数を末尾にペースト
   * 引数:
   *   in[]  : データ配列
   *   N     : サンプル数（2,4,8,16,32,64,128...推奨）
   *   Frequency: サンプリング周波数(Hz)
   * 2の累乗でない場合は切り捨て。Arduino Nanoでは128サンプル以上はメモリ不足で不可（64推奨）
   * 2021/03/02: Nの型をfloatからintに変更（256サンプル以上対応）
   */
  unsigned int data_pow2[13] = {1,2,4,8,16,32,64,128,256,512,1024,2048};
  int a, c1, f, o, x;
  a = N;
  for (int i = 0; i < 12; i++) {
    if (data_pow2[i] <= a) { o = i; }
  }
  int in_ps[data_pow2[o]] = {};
  float out_r[data_pow2[o]] = {};
  float out_im[data_pow2[o]] = {};

  x = 0;
  for (int b = 0; b < o; b++) {
    c1 = data_pow2[b];
    f = data_pow2[o] / (c1 + c1);
    for (int j = 0; j < c1; j++) {
      x = x + 1;
      in_ps[x] = in_ps[j] + f;
    }
  }
  for (int i = 0; i < data_pow2[o]; i++) {
    if (in_ps[i] < a) { out_r[i] = in[in_ps[i]]; }
    if (in_ps[i] > a) { out_r[i] = in[in_ps[i] - a]; }
  }

  int i10, i11, n1;
  float e, c, s, tr, ti;
  for (int i = 0; i < o; i++) {
    i10 = data_pow2[i];
    i11 = data_pow2[o] / data_pow2[i+1];
    e = 360 / data_pow2[i+1];
    e = 0 - e;
    n1 = 0;
    for (int j = 0; j < i10; j++) {
      c = cosine(e * j);
      s = sine(e * j);
      n1 = j;
      for (int k = 0; k < i11; k++) {
        tr = c * out_r[i10 + n1] - s * out_im[i10 + n1];
        ti = s * out_r[i10 + n1] + c * out_im[i10 + n1];
        out_r[n1 + i10] = out_r[n1] - tr;
        out_r[n1] = out_r[n1] + tr;
        out_im[n1 + i10] = out_im[n1] - ti;
        out_im[n1] = out_im[n1] + ti;
        n1 = n1 + i10 + i10;
      }
    }
  }

  // 振幅と周波数ビンの計算
  for (int i = 0; i < data_pow2[o-1]; i++) {
    out_r[i] = sqrt(out_r[i]*out_r[i] + out_im[i]*out_im[i]); // 高速化のためにsqrtを削除可能
    out_im[i] = i * Frequency / N;
  }

  // ピーク検出
  x = 0;
  for (int i = 1; i < data_pow2[o-1] - 1; i++) {
    if (out_r[i] > out_r[i-1] && out_r[i] > out_r[i+1]) {
      in_ps[x] = i;
      x++;
    }
  }

  // 振幅順にソート
  float s_val = 0, c_val = 0;
  for (int i = 0; i < x; i++) {
    for (int j = i; j < x; j++) {
      if (out_r[in_ps[i]] < out_r[in_ps[j]]) {
        s_val = in_ps[i];
        in_ps[i] = in_ps[j];
        in_ps[j] = s_val;
      }
    }
  }

  // 上位5つをグローバル配列に格納
  for (int i = 0; i < 5; i++) {
    f_peaks[i] = out_im[in_ps[i]];
  }
}

// 高速サイン関数（0-360度対応）
float sine(int i) {
  int j = i;
  float out;
  while (j < 0) j += 360;
  while (j > 360) j -= 360;
  if (j >= 0 && j < 91) out = sine_data[j];
  else if (j >= 90 && j < 181) out = sine_data[180 - j];
  else if (j >= 180 && j < 271) out = -sine_data[j - 180];
  else out = -sine_data[360 - j];
  return out / 255.0;
}

// 高速コサイン関数
float cosine(int i) {
  int j = i;
  float out;
  while (j < 0) j += 360;
  while (j > 360) j -= 360;
  if (j >= 0 && j < 91) out = sine_data[90 - j];
  else if (j >= 90 && j < 181) out = -sine_data[j - 90];
  else if (j >= 180 && j < 271) out = -sine_data[270 - j];
  else out = sine_data[j - 270];
  return out / 255.0;
}
//------------------------------------------------------------------------------------//

Annoy-O-Bug：ATtiny85で作るコンパクトなLEDブザーおもちゃ高性能Arduinoライブラリ：DCモーター制御用

製造プロセス

3Dプリント

自動制御システム

産業技術