やりたいこと

STM32F303K8 で DMA を介して DAC から出力する
リアルタイムに任意周波数の波形を合成する
オーディオ出力としてその波形を収録する

前提として NUCLEO-F303K8 を使い、PLL によるシステムクロックを 64 MHz に引き上げて使います。

バッファリング

F303K8 は FPU を搭載していますが、fast_math.h の sinf 関数はリアルタイム波形合成に使えるほど高速ではありません。システムクロックを 64 MHz に引き上げてもサンプリング周波数が 2 kHz 程度の速さしか達成できません。

そこで今回は sinf 関数をリアルタイムの波形合成に直接使わず、LUT を使った方法で種々の波形を合成させます。

LUT

計算に時間のかかる値をあらかじめ計算しておき、配列に格納します。値を使うときは位相から配列のインデクスを求め、配列から値を取得します。これを LUT (Lookup Table) と呼びます。値の精度は配列の長さとその補完方法に依存します。

配列に格納する値は実数値である必要はありませんが、インデクスは整数でなくてはならないため、位相をどう整数で表すかが鍵となります。今回は高速化を優先するため、正弦波の位相を uint32_t で表し、配列の長さを 2の冪にします。同様の理由で、LUT に格納する正弦波は 1周期とします。

正確には、位相と配列長は 2の冪にする必要はありません。ただし 2の冪以外の値を使う場合は位相値のクランプや実数演算、余剰演算が必要になります。

以下に正弦波の LUT の実装を示します。

lookup_sin

#include "fastmath.h"

#define BV(x) (1 << (x))
#define PI 3.141592654

#define LUT_BIT 10
#define LUT_SIZE BV(LUT_BIT)
#define LUT_PHASE(x) ((x) / BV(32 - LUT_BIT))

int16_t lut_sine[LUT_SIZE] = { 0 };

int32_t lookup_sin(uint32_t phase) {
  return lut_sine[LUT_PHASE(phase)];
}

void initialize() {
  for (uint32_t i = 0; i < LUT_SIZE; i++) {
    lut_sine[i] = (int16_t)(sinf(i * 2.0 * PI / LUT_SIZE) * (BV(16 - 1) - 1));
  }
}

int main(void) {
  ...
  initialize();
  ...
}

マクロ BV(x) は 2の冪を表します。 BV(10) で 1024 です。

マクロ LUT_PHASE(x) は uint32_t の位相を LUT の位相（インデクス）に変換します。32 ビットから 10 ビットへの除算となり、コンパイルすると 22 ビットの右シフトになります。

配列 lut_sine[LUT_SIZE] は const にしません。初期化の時間コストはかかりますが、実行時に RAM から LUT を予め展開したほうが 20% ほど高速です。

LUT は int16_t ですが、高速化を狙って関数 lookup_sin の返却値は int32_t にしました。

出力周波数の計算

関数 lookup_sin(uint32_t phase) の引数に与える位相を計算します。

位相を uint32_t の範囲で与えるため、最小値は $0$、最大値は $2^{32} - 1$ です。正弦波の 1周期を転送するのにサンプリング周波数の逆数 $\frac{1}{f_s}$ [s] かかるので、$f$ [Hz] の正弦波を出力するのに必要な1サンプルあたりの位相の増分 $\Delta\theta$ は、

$\Delta\theta = \frac{2^{32}}{f_s} \cdot f$

と計算できます。サンプリング周波数 $f_s$ [Hz] は固定とすると、$\frac{2^{32}}{f_s}$ の部分は定数となり、$f$ と乗算すれば $\Delta\theta$ がすぐに求まります。

以下に正弦波のデータをリアルタイムに生成するコードを示します。

#define DAC_BIT 12
#define VALUE_VU12(x) ((x) / (BV(16 - 1) / BV(DAC_BIT - 2)) + BV(DAC_BIT - 1))
#define BUFFER_SIZE BV(10)

uint16_t buffer[BUFFER_SIZE] = { 0 };

const float sampling_frequency = 48.0e3;
const float phase_factor       = 4294967296.0 / sampling_frequency;
float frequency                = 1000.0;

uint32_t phase       = 0;
uint32_t delta_phase = (uint32_t)round(phase_factor * frequency);

void fill_buffer(uint32_t start, uint32_t end) {
  for (uint32_t i = start; i < end; i++) {
    buffer[i] = VALUE_VU12(lookup_sin(phase));
    phase += delta_phase;
  }
}

マクロ VALUE_VU12(x) は 32 ビット整数の出力を DAC 出力に合わせて 12 ビット整数に変換します。

変数 frequency で周波数を変更できます。実行中に変更した場合は変数 delta_phase を再計算してください。

DMA 割り込み

DMA が RAM から DAC への転送を完了させると、関数 HAL_DAC_ConvCpltCallbackCh1 が実行されます。この関数は __weak が付与されているため、ユーザ側で上書きできます。

しかし、割り込みが発生してもなお DMA の転送は停止せず、バッファの冒頭にある古いデータを DAC に転送し始めます。つまり、波形合成中のデータを DAC に送ってしまうため波形が崩れてしまいます。

実は DMA 割り込みにはもう一つ、関数 HAL_DAC_ConvHalfCpltCallbackCh1 が存在します。この関数はバッファの半分の転送が終わったら実行されます。つまり同じバッファを2分して、ダブルバッファリングが可能になります。

DMA での転送が開始して前半の転送が終わると、関数 HAL_DAC_ConvHalfCpltCallbackCh1 が呼び出されます。後半データの転送中に、関数内で前半データの生成を行います。後半データの転送が終わると関数 HAL_DAC_ConvCpltCallbackCh1 が呼び出されるので、前半データの転送中に、後半データの生成を行います。こうすることで、切れ目のないデータをリアルタイムに DAC に転送できます。

#include "stm32f3xx_hal.h"

// 前半データ転送終了
void HAL_DAC_ConvHalfCpltCallbackCh1(DAC_HandleTypeDef *hdac) {
  HAL_GPIO_WritePin(GPIOB, GPIO_PIN_4, GPIO_PIN_RESET);
  HAL_GPIO_WritePin(GPIOB, GPIO_PIN_3, GPIO_PIN_SET);

  // 前半データ生成
  fill_buffer(0, BUFFER_SIZE / 2);

  HAL_GPIO_WritePin(GPIOB, GPIO_PIN_3, GPIO_PIN_RESET);
}

// 後半データ転送終了
void HAL_DAC_ConvCpltCallbackCh1(DAC_HandleTypeDef *hdac) {
  HAL_GPIO_WritePin(GPIOB, GPIO_PIN_4, GPIO_PIN_SET);
  HAL_GPIO_WritePin(GPIOB, GPIO_PIN_3, GPIO_PIN_SET);

  // 後半データ生成
  fill_buffer(BUFFER_SIZE / 2, BUFFER_SIZE);

  HAL_GPIO_WritePin(GPIOB, GPIO_PIN_3, GPIO_PIN_RESET);
}

今回は DMA 割り込みがされるごとに GPIO_PIN_4 を反転、データ生成中は GPIO_PIN_3 に H を出力して波形合成の負荷率を計測します。

実装

実行すると、1 kHz の正弦波が PA4 （A3）ピンから、負荷率を表す信号が PB13 (D13) ピンから、DMA の割り込みを表す信号が PB4 (D12) から出力されます。

黄色(ch1)が正弦波出力、水色(ch2)が負荷率、紫色(ch3)が DMA 割り込みです。

サンプリング周波数は 48 kHz、バッファサイズは 1024 なので、10.67 ms の余裕があります。正しく 1 kHz の正弦波が出力され、負荷率は 1.9%、時間にして 200 μs となりました。

正弦波以外の波形を作る

LUT に入れる波形を変えれば、正弦波以外の波形を作れます。配列名 lut_sine からは変わるので、本来ならば別の名前を使うべきですが、今回はこのまま使います。

sine_wave

void initialize() {
  for (uint32_t i = 0; i < LUT_SIZE; i++) {
    lut_sine[i] = (int16_t)(sinf(i * 2.0 * PI / LUT_SIZE) * (BV(16 - 1) - 1));
  }
}

矩形波、三角波、鋸歯波を作ります。

square_wave

void initialize() {
  for (uint32_t i = 0; i < LUT_SIZE; i++) {
    lut_sine[i] = (int16_t)((i < LUT_SIZE / 2 ? 1 : -1) * (BV(16 - 1) - 1));
  }
}

triangle_wave

void initialize() {
  for (uint32_t i = 0; i < LUT_SIZE; i++) {
    lut_sine[i] = (int16_t)((i < LUT_SIZE / 4 ? i :
                             i < (LUT_SIZE / 4) * 3 ? LUT_SIZE / 2 - i :
                                                      i - LUT_SIZE)
                            * (BV(16 - LUT_BIT + 1) - 1));
  }
}

sawtooth_wave

void initialize() {
  for (uint32_t i = 0; i < LUT_SIZE; i++) {
    lut_sine[i] = (int16_t)((i < LUT_SIZE / 2 ? i : i - LUT_SIZE) * (BV(16 - LUT_BIT) - 1));
  }
}