音を利用する 4 〜mbed マイコンでの応用：準備編〜
この記事の内容
IoT 方面での応用を想定しマイコン系でも音に関する実験と試作を行いました。今回はその準備編です。ここではプラットフォームとして mbed LPC1768 ボードを使用しています。基本的な考え方は環境に依存しないため要素の多くは他の製品向けに読み替えることも可能でしょう。
準備 1. 部品を用意する
音声データの入出力用に以下の部品を用意しました。
オペアンプ内蔵
マイクモジュールを mbed LPC1768 の p20（AnalogIn）に、オーディオアンプを p18（AnalogOut）に接続しました。これを実験に使います。
準備 2. mbed マイコンでの音声入出力処理を覚える
何はともあれ音声データの入出力方法からです。リファレンスを確認します。
- AnalogIn - アナログ入力 - developer.mbed.org
- AN10974 LPC176x/175x 12-bit ADC design guidelines（PDF） - www.nxp.com
Page 3 of 19 より - AnalogOut - アナログ出力 - developer.mbed.org
val = name.read();
valには0(0V)〜1.0(3.3V)の数値（float）が代入されます。
val = name.read_u16();
valには0(0V)〜65535(0xffff, 3.3V)の数値が代入されます。（上位12ビットが有効）
The LPC175x/6x family is based on the ARM Cortex M3 core, and includes a 12-bit Analog-to-Digital (ADC) module with input multiplexing among eight pins, conversion rates up to 200 kHz, and multiple result registers.
name = 0.5;
0(0V)〜1(3.3V)の数値（float）で出力を指定します。
name.write_u16 = 1000;
0(0V)〜65535(0xffff, 3.3V)の数値（unsigned short）で出力を指定します。（上位10ビットが有効）
試作 1
まず、ごく素朴な録音・再生プログラムを書いてみました。次の内容です。
- サンプリングレートは 4000Hz
- 入力信号を 4000Hz でサンプリングするための時間間隔 w を計算
w = 1sec / 4000Hz = 0.00025sec = 0.25msec = 250μsec - アナログ入力から w 間隔で読み取った値をサンプリング用バッファの終端まで順次格納する
- サンプリング用バッファに格納した値を w 間隔で順次アナログ出力へ書き込む 〜無限ループ
動作の様子
動画：20秒 ※音量注意
ソースコード
sonic01_mic
- main.cpp
#include "mbed.h" #define SAMPLE_RATE 4000 #define BUFSIZE 7680 AnalogIn mic(p20); AnalogOut speaker(p18); DigitalOut led(LED4); float recBuf[BUFSIZE]; int main() { // 当該サンプリングレートでのサンプリング間隔 float w = 1.0 / SAMPLE_RATE; // 録音開始までに LED を 3回点滅させる led = 1; for (int i = 0; i < 6; i++) { wait(0.5); led = !led; } // 録音 for (int i = 0; i < BUFSIZE; i++) { recBuf[i] = mic.read(); wait(w); } led = 0; // 再生 while (1) { for (int i = 0; i < BUFSIZE; i++) { speaker = recBuf[i] ; wait(w); } } return 0; }
試作 1 の問題点について
このようにとても短いコードで簡単に音声データの入出力を実現できました。しかし、録音と再生を試していると微妙な違和感がありました。録音した音と再生される音の高さが違うような気がしたのです。このことを楽器を使って確認してみました。
再生音の高さがおかしい？
動画：21秒 ※音量注意
この動画では楽器の音を録音し再生時に合わせ弾きなどを行っています。なお、前掲の安価なスピーカーではやはり音質が貧弱なので、この動画以降ではアンプ内蔵のスピーカー（PHILIPS 社製 BT50）を AnalogOut ポートにつないで再生に使用しています。
問題の原因と対処方法
このおかしな現象の原因は試作 1 のコードにおいてサンプリング間隔を制御している wait() 関数にあります。mbed プログラミングにおいてきわめてカジュアルに使用されるこの関数のリファレンスをあらためて確認したところ「When you call wait your board's CPU will be busy in a loop waiting for the required time to pass」との記述があり、とりあえず以下の内容で精度の確認を行いました。
インターバル 250μ秒での結果は次の通りです。
このように wait() + AnalogIn.read() で大幅な遅延が発生しています。一方、wait() + AnalogOut.write() では相対的に遅延が小さいため試作 1 のプログラムは結果として「低いレートでサンプリング」したデータを「高いレートで再生」する所作となっていました。上の検証でのインターバル 10〜250μ秒を横軸、所要時間の誤差の割合を縦軸にプロットしたグラフを以下に示します。
（クリックで可読大表示）
サンプリング間隔のコントロールに mbed 標準のタイマ割込み機構である Ticker を使って同じ検証を行いました。
インターバル 250μ秒での結果は次の通りです。
（クリックで可読大表示）
このように結果は wait() 使用時よりも圧倒的に良好です。その一方で、AnalogIn.read() 併用時のインターバル 25μ秒付近の「壁」の存在が目に止まります。手元の検証ではこの結果には明確な再現性がみられました。25μ秒以下の割込みを必要としないサンプリングレートを利用すれば大きな支障はなさそうですが、広く用いられるレート 44100Hz の場合には 1sec / 44100Hz ≒ 23μsec とサンプリング間隔が 25μ秒を割り込むことが微妙に残念でもあります。
一連の調査の過程で Erik Olieman 氏による FastAnalogIn ライブラリの存在を知りました。
所定のアナログ入力ポートに対する mbed 標準の AnalogIn.read() が、呼び出しの都度、当該ポートに対応する ADC チャネル経由で AD 変換を行う仕様（そのため遅い）であるのに対し、FastAnalogIn は ADC の BURST mode を利用して変換をバックグラウンドで最大 200kHzで繰り返し実行し read() が呼び出されると「最後にサンプリングした値」を直ちに返す内容で実装されています。 （※ 前掲の ADC ドキュメント中の "conversion rates up to 200 kHz" のくだりは BURST mode 下のこの所作を指すものと考えられる）
- FastAnalogIn - a mercurial repository | mbed - developer.mbed.org
- FastAnalogIn Class Reference - developer.mbed.org
- ADC_LPC1768 - John Kneen: Microcontrollers - 6. Selecting and triggering the analogue to digital conversion. - sites.google.com
When you read an AnalogIn object it will enable the corresponding ADC channel, depending on the implementation do either one or multiple measurements for more accuracy, and return that value to your program. This way the ADC is only active when it is required, and it is fairly straightforward. However the downside is, is that an ADC is relatively slow.
:
This library uses the 'burst' feature of the microcontroller. This allows the ADC on the background to perform the AD conversions without requiring intervention from the microcontroller's core. Also there are no interrupts used, so also your time-sensitive code is not affected.
:
AnalogIn does a single conversion when you read a value (actually several conversions and it takes the median of that). This library runns the ADC conversion automatically in the background. When read is called, it immediatly returns the last sampled value.
:
Bit 16 : BURST
1 The AD converter does repeated conversions at up to 200 kHz, scanning (if necessary) through the pins selected by bits set to ones in the SEL field.
:
計測用コードの中の AnalogIn をこの FastAnalogIn に変更して同様の計測を行ったところ 25μ秒以下のインターバルでの成績が劇的に改善されました。前掲のふたつのグラフの縦軸最大値が「100%」であったのに対し、下のグラフでは「1.00%」であることに要注目です。
（クリックで可読大表示）
試作 2
上の「対処方法」を反映した形で録音・再生プログラムを書き直しました。検証上の便宜を兼ねて次の内容としています。
- 録音・再生に要素数 14000 の unsigned short 型配列を使用する
- サンプリングレート 8000Hz, 16000Hz, 22000Hz, 30720Hz, 44100Hz の順に録音と再生を行う
- レートが高くなるほど録音可能時間は短くなるが、再生はいずれも 8秒間程度繰り返し行う
動作の様子
動画：1分28秒 ※音量注意
この動画では各サンプリングレートにおいてキーボードの「ファ」の音（F5：698.456Hz）を録音し、その再生音を以前作成したピーク周波数表示プログラムと楽器用のチューナー（KORG 社製 AW-2G）で確認しています。
ソースコード
sonic07_mic3
- main.cpp
// サンプリングと再生 // レート 8000,16000,22000,30720,44100Hz の順にサンプリングを行い // それぞれ 8秒程度ループ再生する // FastAnalogIn ライブラリを使用 #include "mbed.h" #include "FastAnalogIn.h" #define DEVELOP #ifdef DEVELOP #define dbg(...) printf(__VA_ARGS__) #else #define dbg(...) #endif #define BUFSIZE 14000 Ticker tick; Timer timer; int countTick,interval_us, sampleTime; bool recDone; unsigned short *recBufp = NULL; AnalogOut speaker(p18); FastAnalogIn mic(p20); DigitalOut led(LED4); // サンプリング void flipRec() { recBufp[countTick++] = mic.read_u16() & 0xFFFC; if (countTick == BUFSIZE) { tick.detach(); recDone = true; } } // 再生 void flipPlay() { speaker.write_u16(recBufp[countTick++]); if (countTick == BUFSIZE) { countTick = 0; } } void recAndPlay(int sampleRate) { // 当該サンプリングレートでのサンプリング間隔 usec interval_us = 1.0 / sampleRate * 1000000; // 当該サンプリングレートでのサンプリング時間 usec sampleTime = BUFSIZE * interval_us / sizeof(short); dbg("rate=%5dHz irq=%dus time=%dms\r\n", sampleRate, interval_us, sampleTime/1000); // サンプリング開始までに LED を 3回点滅させる led = 1; for (int i = 0; i < 6; i++) { wait(0.4); led =! led; } // サンプリング countTick = 0; recDone = false; tick.attach_us(&flipRec, interval_us); while (!recDone) { wait(0.2); } led = 0; wait(1); // 8秒程度ループ再生 countTick = 0; tick.attach_us(&flipPlay, interval_us); wait(8); tick.detach(); } int main() { recBufp = (unsigned short*)malloc(BUFSIZE); recAndPlay(8000); recAndPlay(16000); recAndPlay(22000); recAndPlay(30720); recAndPlay(44100); dbg("done\r\n"); if (recBufp) { free(recBufp); } return 0; }
準備 3. 任意の周波数の音波を生成・出力する
先般のAndroid 環境での試みと同様に手元の mbed マイコンでサイン波信号を生成・出力してみます。Android 版での到達点を踏襲できるのであまり手はかかりませんでした。
試作
以下の要領でプログラムを作成しました。
- 400Hz〜14000Hz の範囲の 10種類の周波数のサイン波信号を順次生成する
- サンプリングレートは 33330Hz（サンプリング間隔 30μ秒）とし、バッファ上に各周波数 100ミリ秒分の信号データを構築してそれぞれ 4秒間程度 AnalogOut ポートへ出力する
- unsigned short 型配列を引数にとる AnalogOut.write_u16() の引数範囲は 0〜65535 であるため信号生成時に直流成分として 65535 / 2 を加える
動作の様子
動画：1分15秒 ※音量注意
この動画では二台のタブレット上のふたつのプログラムを使って信号の確認を行っています。
- 左側：手元で以前作成したプログラムでピーク周波数を表示
- 右側：Google Play ストアで公開されている AudioUtil Spectrum Analyzer で波形を表示
ソースコード
sonic08_sine2
- main.cpp
// サイン波音声信号を生成・出力する #include "mbed.h" #define DEVELOP #ifdef DEVELOP #define dbg(...) printf(__VA_ARGS__) #else #define dbg(...) #endif #define SAMPLE_RATE 33330 #define PI 3.141592653589793238462 #define MAG 3000 Ticker tick; int count = 0; const int BUFSIZE = SAMPLE_RATE / 10; // バッファサイズ 100ミリ秒分 unsigned short waveBuf[BUFSIZE]; // 当該サンプリングレートでのサンプリング間隔 秒 const float SEC_PER_SAMPLEPOINT = 1.0 / SAMPLE_RATE; const float a0 = 65535 / 2; AnalogOut speaker(p18); // 出力 void flip() { speaker.write_u16(waveBuf[count++]); (count >= BUFSIZE) ? count = 0 : 1; } void outputSineWave(int freq) { for (int i = 0; i < BUFSIZE; i++) { float currentSec = i * SEC_PER_SAMPLEPOINT; double val = a0 + MAG * sin(2.0*PI * freq * currentSec); waveBuf[i] = (unsigned short)val; } count = 0; tick.attach(&flip, SEC_PER_SAMPLEPOINT; wait(4); // 4秒程度出力 tick.detach(); } int main() { dbg("interval=%fs\r\n", SEC_PER_SAMPLEPOINT); outputSineWave(400); outputSineWave(800); outputSineWave(1200); outputSineWave(1600); outputSineWave(2000); outputSineWave(4000); outputSineWave(8000); outputSineWave(10000); outputSineWave(12000); outputSineWave(14000); return 0; }
準備 4. 集音したデータの周波数分析を行う
準備編の最後に入力信号の分析を試みます。
大浦版 FFT について
ここでは高性能かつコンパクトな FFT ライブラリとして広く利用（NASA の例）されている大浦拓哉様（京都大学数理解析研究所 助教）によるライセンスフリーの汎用 FFT パッケージ（通称 大浦版 FFT）を使ってみることにしました。先般 Android 環境下で利用した JTransforms はこのパッケージを母体とするソフトウェアです。
- 大浦 拓哉 (Ooura, Takuya) - 京都大学数理解析研究所 - www.kurims.kyoto-u.ac.jp
- 音声信号は実数データにつき実離散フーリエ変換関数 rdft() を使用する。使い方はシンプル。FFT サイズを 2のべき乗の値とすることに注意。
-------- Real DFT / Inverse of Real DFT --------
[usage]
<case1>
ip[0] = 0; // first time only
rdft(n, 1, a, ip, w);
<case2>
ip[0] = 0; // first time only
rdft(n, -1, a, ip, w);
[parameters]
n :data length (int)
n >= 2, n = power of 2
a[0...n-1] :input/output data (double *)
<case1>
output data
a[2*k] = R[k], 0<=k<n/2
a[2*k+1] = I[k], 0<k<n/2
a[1] = R[n/2]
<case2>
input data
a[2*j] = R[j], 0<=j<n/2
a[2*j+1] = I[j], 0<j<n/2
a[1] = R[n/2]
ip[0...*] :work area for bit reversal (int *)
length of ip >= 2+sqrt(n/2)
strictly,
length of ip >=
2+(1<<(int)(log(n/2+0.5)/log(2))/2).
ip[0],ip[1] are pointers of the cos/sin table.
w[0...n/2-1] :cos/sin table (double *)
w[],ip[] are initialized if ip[0] == 0.
[remark]
Inverse of
rdft(n, 1, a, ip, w);
is
rdft(n, -1, a, ip, w);
for (j = 0; j <= n - 1; j++) {
a[j] *= 2.0 / n;
}
一次元 DFT / DCT / DST
詳細
一次元，２の整数乗個のデータの離散 Fourier 変換， コサイン変換，サイン変換等を行います．このパッケージには C と Fortran の FFT コードが含まれます．
:
内容
fft4g.c : FFT パッケージ in C - 高速版 (基数 4, 2)
:
ファイルの違い
簡易版は作業用配列を一切使いません．高速版は三角関数表などの作業用配列を使います．高速版は簡易版に比べて高速かつ高精度．高速版の仕様はすべて同じ．
パッケージ内のルーチン
cdft() : 複素離散 Fourier 変換 rdft() : 実離散 Fourier 変換
:
ライセンス
Copyright Takuya OOURA, 1996-2001
このソースコードはフリーソフトです．商用目的を含め，このコードの使用， コピー，修正及び配布は自由に行って結構です．ただし，このコードの修正を 行った場合は，そのことを明記してください．
試作
レート 33330Hz でサンプリングした音声信号を mbed マイコン上の FFT 処理にかけてピーク周波数を求めコンソールへ連続出力する内容で試作を行いました。手元の確認では mbed LPC1768 のメモリ容量上 FFT サイズは 210 = 1024 が限界でありその影響らしきものが若干窺えるものの十分に実用水準にあると考えられます。
動作の様子
動画：1分33秒 ※音量注意
この動画では二台のタブレット上のふたつのプログラムを使ってテストを行っています。
- 左側：以前作成したピーク周波数表示プログラム
- 右側：以前作成したサイン波信号出力プログラム
ソースコード
sonic05_FFT2
- main.cpp
// 大浦版 FFT ライブラリで入力音声信号のピーク周波数を得る #include "mbed.h" #define DEVELOP #ifdef DEVELOP #define dbg(...) printf(__VA_ARGS__) #else #define dbg(...) #endif #define PI = 3.141592653589793238462; #define SAMPLE_RATE 33330 // n >= 2, n = power of 2 #define FFTSIZE 1024 // length of ip >= 2+sqrt(n/2) const int IP_LEN = 4 + sqrt((double)(FFTSIZE/2)); Ticker tick; double *recBufp = NULL; int countTick; Timer timer; bool recDone; AnalogIn mic(p20); DigitalOut led(LED4); // 当該サンプリングレートでのサンプリング間隔 秒 const float SEC_PER_SAMPLEPOINT = 1.0 / SAMPLE_RATE; void rdft(int n, int isgn, double *a, int *ip, double *w); int doFFT(double *FFTBuffer) { int ip[IP_LEN]; double *w = (double*)malloc(sizeof(double) * (FFTSIZE/2)); ip[0] = 0; rdft(FFTSIZE, 1, FFTBuffer, ip, w); double maxAmp = 0; int index; for (int i = 0; i < FFTSIZE/2; i++) { double a = FFTBuffer[i*2]; // 実部 double b = FFTBuffer[i*2 + 1]; // 虚部 // a+ib の絶対値 √ a^2 + b^2 = r が振幅値 double r = sqrt(a*a + b*b); if (r > maxAmp) { maxAmp = r; index = i; } } double freq = index * SAMPLE_RATE / FFTSIZE; if (freq > 0 && maxAmp > 1.0) { dbg("freq=%d mag=%f\r\n", (int)freq, maxAmp); } free(w); return (int)freq; } // サンプリング void flip() { // AnalogIn.read() の値範囲は 0〜1.0 につき中心値を 0 に recBufp[countTick++] = mic.read() - 0.5; if (countTick == FFTSIZE) { tick.detach(); recDone = true; } } int main() { recBufp = (double*)malloc(sizeof(double) * FFTSIZE); while (1) { countTick = 0; recDone = false; tick.attach(&flip, SEC_PER_SAMPLEPOINT); while (!recDone) { wait_ms(100); } doFFT(recBufp); } if (recBufp) { free(recBufp); } return 0; }
(tanabe)
