勝手な電子工作・・

勝手なオリジナル電子工作に関する記事を書きます

楽々フーリエ変換!(快速マイコンをArduino-IDEで使う)

f:id:a-tomi:20210405120325j:plain

トライアングルの音色は澄んでいるようで意外に複雑です。今回はマイコンでできるフーリエ変換をご紹介したいと思います。

上のグラフはトライアングルの音をマイクで採りつつフーリエ変換した結果で、周波数の成分を表しています。なんどやっても同じグラフになります。たぶんトライアングルには振動の腹が主に3つあるからだと思われます。とはいえ百均で入手したおもちゃなので、プロが使う本物とはそれなりに違うかもしれませんが^^;

f:id:a-tomi:20210405121057j:plain

このトライアングルの音の波形は次です。

f:id:a-tomi:20210405121555j:plain

よい音が長く続きますが、意外に複雑な波形です。

これをマイコンでどうやってフーリエ解析して、周波数成分の分布を知るかを説明します。また、これでいかに正確な分析ができるかについても、ついでにテストしてみたのでご紹介します。

今回の回路は次です。

f:id:a-tomi:20210405171049j:plain



快速マイコンTeensy4.0をArduinoIDEでプログラミングします。ここでの勝手なスケッチは次です。

/*****************************************************************
 *  Fast Foourier Transform & Data Recording, using Teensy 4.0
 *    Initial version V00:  April 3, 2021 (c) Akira Tominaga
 *    Function:
 *      1) Quick sampling and FFT, applying Hamming-window
 *      2) Show FFT graph, with Serial Plotter
 *      3) Record measured data to micro SD with name "FFTnnn.txt"
 *      4)Perform next measurement if tact-switch pushed
 *    Input and output: 
 *      1) Pin A2(16) for signal input (keep within 0 to 3.3V)
 *      2) Pin 10,11,12 and 13 for SPI interface to micro SD
 *      3) Pin 2 for tact-switch to continue to next FFT
 *    Remarks: Arduino-FFT library is used under GNU-GPLicense.
 *    Revisions: (not yet)
 *****************************************************************/
#include "SD.h"       // use SPI micro-SD-card-reader/writer 
#include "SPI.h"      // SPI MOSI=11, MISO=12,CLK=13
const int cS = 10;    // SPI Chip-Select = 10
File sdF;
#include "arduinoFFT.h" // use Fast Fourier Transform
// for 0~20kHz FFT, 40kHz sampling & 1024 samples=512 bins
#define sFreq 40000   // sampling Hz <55000 (ADC 17.5 micro sec)
#define nSamp 1024    // number of samplings 
#define serB 115200   // serial Baud 115.2,57.6,38.4,19.2,14.4,9.6 k
#define aPin  A2      // signal input pin
#define sW 2          // tact switch at D2
arduinoFFT FFT = arduinoFFT();
unsigned int sDur;    // duration micro sec. between samplings
unsigned long micS;   // elapsed time in microsec (max about 70min)
double vR[nSamp];     // Real part of Vector
double vI[nSamp];     // Imaginal part of Vector

void setup() { // ***** Teensy setup() *****
  Serial.begin(serB);
  while (!Serial) {}  // wait for Serial-ready
  pinMode(sW, INPUT_PULLUP); // pull-up tact switch
  sDur = round(1000000 * (1.0 / sFreq)); // calc. sampling dur.μS
  if (!SD.begin(cS)) { // start SD
    Serial.println(F("Init-SD err."));  return;
  }
}

void loop() { // ***** Teensy loop() *****
  // *** perform Fast Fourier Transform
  for (int i = 0; i < nSamp; i++) { // do a cycle of samplings
    micS = micros();  // save sampling-start-time
    vR[i] = analogRead(aPin);
    vI[i] = 0;        // set zero to imaginal part
    while (micros() < (micS + sDur)) {} // wait until cycle-end
  }
  // apply Hamming windowing
  FFT.Windowing(vR, nSamp, FFT_WIN_TYP_HAMMING, FFT_FORWARD);
  FFT.Compute(vR, vI, nSamp, FFT_FORWARD); // calculate FT
  FFT.ComplexToMagnitude(vR, vI, nSamp); // mag. to vR 1st half
  double peak = FFT.MajorPeak(vR, nSamp, sFreq); // get peak freq
  // get peak element number
  int ip = round((nSamp / 2) * peak * 1.0 / (sFreq / 2));
  
  // *** prepare SD to create file=FFTnnn.txt
  char fN[11] =  "FFTnnn.txt"; // file name
  uint16_t fnum = 0;  // file number for character nnn
  char fnumc[4];      // file number in character
  while (true) { // if the name existsts, then set nnn+1
    sprintf(fnumc, "%03d", fnum);
    for (int i = 0; i < 3; i++) {
      fN[i + 3] = fnumc[i];
    }
    if (!SD.exists(fN)) break; // if not exists, go out
    fnum ++; // if it exists, increase nnn
  }
  sdF = SD.open(fN, FILE_WRITE); // create it as output file
  if (!sdF) {		// if error, then
    Serial.println(F("Open-SD err."));
    while (true) {} // if open-err, stop here
  }
  // convert results in vR[] into simple integer hh to record to SD
  for (int i = 0; i < (nSamp / 2); i++) {
    uint16_t ff = round(i * 20000.0 / (nSamp / 2.0)); // get Hz value
    uint16_t hh = round(10000.0 * (vR[i] * 1.0 / vR[ip])); // relative Mag.
    Serial.println(hh); // plot to Serial-plotter
    // *** record( Hz, Mag.) to SD in CSV format
    sdF.print(ff);
    sdF.print(",");
    sdF.println(hh);
  }
  sdF.close();		// close the file when all data recorded
  while (digitalRead(sW) == HIGH) {} // wait for sW pushed
  // *** Peform next measurement and recording when sW pushed
}
// end of program

 

Arduino-FFT高速フーリエ変換)ライブラリーは、もともとGoogle社が開発し、今はArduinoIDEのほかにGCCなどでも使える素晴らしいライブラリーです。ArduinoIDEにこれを追加する方法は普通の入れ方と同じで、「ツール」→「ライブラリーを管理」→「ライブラリーマネジャー画面」で検索バーを使ってFFTを検索し「ArduinoFFT」を選び→該当の「インストール」をクリックします。

今回使っているマイコンTeensy4.xは、フローティング処理も含めて非常に速いので、こういう用途にはぴったりです。もちろん、このFFTはUNOなどでも使えるには使えるわけですが、非常に低い周波数領域しかカバーできません。

サンプリング周波数を40kHzにしている(44kHzもよいかな)理由は、可聴音域の20kHzまでを分析したいからです(FFTではサンプリング周波数の半分までしか分析できません)。

また、Arduino-IDEのシリアルプロッターを使う場合は、横軸を512点しか表示しませんので横軸をそれに合わせるために、サンプリング個数を1024にしています。FFT関数は周波数分解能であるその半分の512個(単位はビン)にわけて答を返してきます。上のプログラムでおわかりかもしれませんが、関数に引き渡すパラメータベクトルのリアル部分用であるvR[ ]の前半分に収めて返してくれます。

 信号測定にはADCを使います。Teensyの迅速AD変換方法はあるようですが、普通のAD変換で実測した結果17.5μS以下でできるので、この用途(40kHzでのサンプリングと処理)用には十分に速いと分かります。

最初に示したトライアングルの測定等では、コンデンサーマイクを用いています。ここでは次のものを使っています。

f:id:a-tomi:20210405135000j:plain

これを信号入力部分につぎのように接続します。

f:id:a-tomi:20210405135214j:plain

マイクは適当な保持具にもたせます。

f:id:a-tomi:20210405135301j:plain

そして動かしシリアルプロッターに表示すると、トライアングルの音の周波数は次のように表示されます。

f:id:a-tomi:20210405135444j:plain

 X軸は20kHzまでの周波数として、等間隔で512点をプロットしますが、横軸の目盛り表示はおしきせで思うに任せません。

また、左端の12点(つまり0~469Hz)は縦軸がかぶさり表示されません(ただし、左端には商用電源のノイズなどが載りやすいため、表示が要らないといえば要らないわけなのですが^^;)余談ですが、上のプログラムでは、peak = FFT.MajorPeak(vR, nSamp, sFreq);により支配的な周波数を取り出しています。試しに50Hzをまぜてみた限り、この(左端に近い)あたりのBinの値がたとえ一番大きくなっても、FFT.MajorPeak関数ではそこはピークとみなされない仕様になっている感じです。

じっくり分析するには、シリアルプロッターではなくシリアルから信号を表示するプログラムを作るか、測定結果をSDにでも記録してExcelで分析するのがよいかということで、とりあえずは各回の測定記録をマイクロSDに書くようにしたわけです。

さて、このフーリエ変換の精度はどうかを調べたくなりますね。

そこで、ファンクションジェネレータの正確な波形を信号として入れてみることにします。まずは2.000kHzの正確な正弦波(sin波)です。

f:id:a-tomi:20210405142152j:plain

オシロで表示すると次。

f:id:a-tomi:20210405142231j:plain

上の写真の右下をごらんになると、正確に2.0000kHzであるのがわかりますね。最後の0.02Hzの半端はおそらくこのオシロ側での有効数字6桁末尾の許容誤差だと思います。

これを測定した結果は次です。

f:id:a-tomi:20210405142604j:plain

高調波の全くない、実に正確なサイン波であることが納得できますし、逆に今回作ってみたFFT装置の精度も高いものであることが理解できます。 

では、同じ周波数のブロック波(square波)をいれてみます。

f:id:a-tomi:20210405142807j:plain

オシロの波形は次です。

f:id:a-tomi:20210405142842j:plain

これを測定すると次の結果になりました。

f:id:a-tomi:20210405142928j:plain

理屈通りの3倍、5倍、7倍、9倍の周波数成分が正確に出ていることがわかります!

(2021年4月8日追記)上のグラフの横軸の分解能はサンプリング数÷2の512ビンで、この場合はビン幅約39Hzへの確率分布で、いわば少しならされた値になります。論理的には幅のない成分の高さは、そのため若干下に表示されるわけです。

例えば、3倍の6kHzなら論理的には3.33、5倍の10kHzなら2.00ですが、そういうわけで表示は少しだけ低くなります。ビン幅を狭めるために、ためしにサンプリング数を16倍の16384回にしてみます。上のプログラムなら#define中の1行を、#define nSamp 16384 に変えるだけのことですのでやってみましたら次のグラフになりました。

f:id:a-tomi:20210408151501j:plain

横軸の分解能は2.44Hzになり波形は前より尖り、高さも理論値に少し近づくのがわかりますね。(以上2021.4.8追記)

 どっちにしてもファンクションジェネレータが正確な事と、今回のFFTの精度が高い点の両方がよくわかりました^^ 

それでは・・・、ということで気をよくして、今度はトイピアノを持ち出してきて測ってみます。

f:id:a-tomi:20210405143342j:plain

一番上の「ド」を鳴らします。

オシロでみると次です。

f:id:a-tomi:20210405143500j:plain

なんだか頼りない音ですね。でも測定結果は次です。

f:id:a-tomi:20210405143544j:plain

おお、あるべき基本周波数1046.5Hzを実に正確に出しています!驚き。高調波も1オクターブ上の2093Hz、そして2オクターブ上の4186Hzに!おもちゃピアノといっても馬鹿にできませんね。

今や玩具でもクリスタル付きマイコンで電子音を出してるわけですから、むしろあたりまえですかね^^

 

これに気をよくして、次はファンクションジェネレータからこれぞホワイトノイズのAWGNを出してこれにいれてみます。

f:id:a-tomi:20210405144507j:plain

f:id:a-tomi:20210405144524j:plain

そして測ってみると;

f:id:a-tomi:20210405144553j:plain

ありゃ、平たんに近いのではなくおおいに癖があるのですね。まあ、人間が数式でこしらえている波である限りは、一部を取り出せばこうなるのがむしろ当たり前かもしれませんね。しかしこれはちょっと拍子抜け・・^^;

 

では今回の記事はこのへんにしたいと思います。皆様のお役に少しでも立てば幸いです。

 

©2021 Akira Tominaga, All rights reserved.