STM32G4のFMACを使ってみる

　STM32G4を始めとしたいくつかのチップにはFMACと呼ばれるフィルタ計算アクセラレータが搭載されています(FMAC: Filter Math ACcelerator)。簡単な使い方を備忘録代わりに残しておきます。

　FMAC自体の説明に関してはSTマイクロが公式に資料を出しています(要ユーザー登録)。
　[STM32 FMACの使い方｜デザイン/サポート｜STM32, STM8ファミリはSTの32bit／8bit汎用マイクロコントローラ製品](https://www.stmcu.jp/design/document/technical_note/72398/)

## 挙動

　今回はRNG(true Random Number Generator)を使用して乱数(ホワイトノイズ)を作成し、それをFMACで適当にフィルタリングし、DAC(Digital to Analog Converter)でアナログ値として出力します。アナログ値をオシロスコープで観察し、フーリエ変換を行うことでフィルタの周波数特性を確認します。
　今回の使い方は、先の資料に則るとシンクドリブンコントロールというタイプの使い方です。FMAの後ろ(シンク)がDMAの転送タイミングを決めてFMACから値を読み取り、FMACは値が読まれ次第、次の計算に必要な値を読みに行く、という流れです。

## CubeMXでの設定

### Computing - FMAC

* Activatedをチェック
* NVICタブ
    * 必要に応じて設定(基本的にデフォルトで問題なし)
    * 処理で問題が起きると割り込みが飛ぶので、FMACの割り込みを有効化しておくとデバッグで便利かも
* DMA Settingsタブ
    * FMAC_WRITを作成
        * Circularモード
        * Increment AddressはPeripheral/Memory共にチェックを外す
        * Data WidthはPeripheral/Memory共にWord

FMAC自体の設定はソフトウェアで行う必要がある。

### Security - RNG

* Activatedをチェック

特に設定項目は無し。

### Analog - DAC1

* OUT1 modeをConnected to external pin onlyに設定
* Parameter Settingsタブ
    * Signed FormatをEnableに
    * Triggerに適当なタイマを指定(今回はTIM7を使用)
* DMA Settingsタブ
    * DAC1_CH1を作成
        * Circularモード
        * Increment AddressはPeripheral/Memory共にチェックを外す
        * Data WidthはPeripheral/Memory共にWord

### Timers - TIM7

* Activatedをチェック
* Parameter Settings
    * PrescalerとCounter Periodを適切な値に設定
        * 今回はわかりやすく1MHz周期になるように設定した
    * Trigger Event SelectionをUpdate Eventに設定


## プログラム

```c++
// do { } while(false);で囲ってある

HAL_StatusTypeDef status = HAL_OK;

int16_t coeff[] = {
    // LPF (0.2, 31taps, blackman)
    // q15フォーマットで0.5倍
    0,
    -3,
    -8,
    21,
    70,
    0,
    -221,
    -222,
    345,
    843,
    0,
    -1851,
    -1735,
    2852,
    9743,
    13107,
    9743,
    2852,
    -1735,
    -1851,
    0,
    843,
    345,
    -222,
    -221,
    0,
    70,
    21,
    -8,
    -3,
    0,
};

FMAC_FilterConfigTypeDef filter_config = {
    .InputBaseAddress = 0,
    .InputBufferSize = 120,
    .InputThreshold = FMAC_THRESHOLD_1,
    .CoeffBaseAddress = 120,
    .CoeffBufferSize = 120,
    .OutputBaseAddress = 240,
    .OutputBufferSize = 16,
    .OutputThreshold = FMAC_THRESHOLD_1,
    .pCoeffB = coeff,
    .CoeffBSize = sizeof(coeff) / sizeof(coeff[0]),
    .InputAccess = FMAC_BUFFER_ACCESS_DMA,
    .OutputAccess = FMAC_BUFFER_ACCESS_POLLING,
    .Clip = FMAC_CLIP_ENABLED,
    .Filter = FMAC_FUNC_CONVO_FIR,
    .P = sizeof(coeff) / sizeof(coeff[0]),
    .R = 1, // 2^1 = 2倍のゲイン
};
status = HAL_FMAC_FilterConfig(&hfmac, &filter_config);
printf("%d %d %lu\n", __LINE__, status, hfmac.ErrorCode);
if (status)
{
    break;
}

status = HAL_FMAC_FilterStart(&hfmac, nullptr, nullptr);
printf("%d %d %lu\n", __LINE__, status, hfmac.ErrorCode);
if (status)
{
    break;
}

uint16_t input_size = 1000;
status = HAL_FMAC_AppendFilterData(&hfmac, reinterpret_cast<int16_t *>(const_cast<uint32_t *>(&hrng.Instance->DR)), &input_size);
printf("%d %d %lu\n", __LINE__, status, hfmac.ErrorCode);
if (status)
{
    break;
}

status = HAL_DAC_Start_DMA(&hdac1, DAC_CHANNEL_1, const_cast<uint32_t *>(&hfmac.Instance->RDATA), 1000, DAC_ALIGN_12B_L);
printf("%d %d %lu\n", __LINE__, status, hdac1.ErrorCode);
if (status)
{
    break;
}

__HAL_DMA_DISABLE_IT(hfmac.hdmaIn, (DMA_IT_TC | DMA_IT_HT));
__HAL_DMA_DISABLE_IT(hdac1.DMA_Handle1, (DMA_IT_TC | DMA_IT_HT));

HAL_TIM_Base_Start(&htim7);
```

## 周波数応答

　フィルタを使用しない場合(DACとRNGを直結した場合)は以下のようになる。

![2021-09-10_20-27-52 - コピー.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/264452/aee04f95-7ae1-fde1-7d65-389ae7ca8a8c.png)

　DACを1Mspsで駆動し、0Hzから1MHzまでを見ているので、DACのSinc特性の山ひとつ分が見えている。

***

　中心周波数0.2fs、帯域幅0.1fs、77タップのバンドパスフィルタを使用した例。

![2021-09-10_21-30-19 - コピー.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/264452/ceb26dae-7ecd-d5a4-ec78-484743e085ce.png)

　ほぼ設計通りの周波数特性が出ている。

***

　Sinc補償フィルタ(帯域幅0.4fs、80タップ、-6dB)を使用した例。

![2021-09-11_02-45-05 - コピー.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/264452/daa02211-8ece-17ec-2b41-af063d7bd9a5.png)

　Sinc補償フィルタは高周波側を強める特性なので、そのままでは固定小数点で表現できる範囲(DACで出力できる範囲)を超えてしまう可能性がある(サチるように設定はしているが)。対策として係数に-6dBの減衰をあたえている(実際はもう少し弱くてもいいが、-6dBなら直感的なので)。
　前述のRNG-DAC直結では平坦な周波数特性とは言い難い山なりの形状だったが、補正フィルタを通した場合はカットオフ周波数までほぼフラットな特性が実現できる。ただし飽和対策のためにゲインを下げている分、DACの分解能が犠牲になっている。また、ナイキャスト周波数を超えた部分は(相対的に)正の利得なので、DAC後段のアンチエイリアシングフィルタはより強い特性が必要になる。ただし0.5fs付近は-20dBくらい落ち込んでいるので、LPFに要求される急峻さは少し緩和される。このバンドストップ特性(幅の広さ)は係数を作るときに設定できるので、後段のLPFが急峻であればさらに高い帯域まで使用したり、LPFの特性を下げて狭い帯域だけ使ったりもできる。

　今回は予め係数を計算してソースコードに埋め込んでいるが、単純なバンドパスフィルタ(LPFやHPFを含む)はSinc関数で、Sinc補償フィルタは逆フーリエ変換で係数を作成できるので、マイコン内で係数を作成することもできる(Cortex-M用のFFTライブラリを使用してもいいし、せいぜい100ポイント前後だから離散フーリエ変換で処理しても良いはず)。

## スループット

　STの資料にスループットの計算方法が書いてあるが、FMACの計算速度は素晴らしく早いというわけではない。例えば170MHzのコアで1Mspsを計算する場合は85タップが上限となる。とはいえ、オーディオ信号くらいのサンプリングレートであればスループットは特に問題はないはず(FMACのFIRは最大で127タップまでしか設定できないため、660kSPS@170MHz程度までは設定できる)。
　スループットが足りない場合、シンク側がDMA転送を要求するとアンダーフロー割り込みが発生する。

## (参考)動的に係数を作成する

### LPF/BPF/HPF

```c++
// do { } while (false);で囲ってある
enum class Filter_mode
{
    none,
    LPF,
    BPF,
    HPF,
} filter_mode = Filter_mode::none;

float center = 0,
        bandwidth = 0,
        gain_dB = 0;
int32_t taps = 0;

if (false)
{
}
else if (2 == sscanf(linebuff, "filter LPF %f %ld", &bandwidth, &taps))
{
    filter_mode = Filter_mode::LPF;
}
else if (3 == sscanf(linebuff, "filter BPF %f %f %ld", &center, &bandwidth, &taps))
{
    filter_mode = Filter_mode::BPF;
}
else if (2 == sscanf(linebuff, "filter HPF %f %ld", &bandwidth, &taps))
{
    filter_mode = Filter_mode::HPF;
}
else
{
    printf("command error\n");
    break;
}

if (center < 0 || 0.5 < center ||
    bandwidth < 0 || 0.5 < bandwidth ||
    taps < 1 || 110 < taps)
{
    printf("parameter error\n");
    break;
}

switch (filter_mode)
{
case Filter_mode::none:
    break;

case Filter_mode::LPF:
    center = 0;
    bandwidth *= 2;
    break;

case Filter_mode::BPF:
    gain_dB += 6;
    break;

case Filter_mode::HPF:
    center = 0.5f;
    bandwidth *= 2;
    break;
}

// 動的に確保したメモリは開放しない(リーク・フラグメンテーション対策)
static int16_t *coeff_buff = nullptr;

if (!coeff_buff)
{
    coeff_buff = static_cast<int16_t *>(pvPortMalloc(sizeof(int16_t) * 127));
}

if (!coeff_buff)
{
    printf("malloc error\n");
    break;
}

class Window_function
{
public:
    static float32_t blackman(float32_t x)
    {
        return x < 0 || 1 < x ? 0 : 0.42f - 0.5f * cosf(2 * pi * x) + 0.08f * cosf(4 * pi * x);
    }
};

const float32_t taps_f = taps;
const float32_t gain = powf(10, gain_dB / 20.0f);
for (int32_t i = 0; i < taps; i++)
{
    const int32_t j = i - taps / 2;
    const float a = j * pi * bandwidth;
    const float mag = (a == 0 ? 1 : sinf(a) / a) * bandwidth;
    const float rad = j * 2 * pi * center;
    coeff_buff[i] = static_cast<int16_t>(
        roundf(0x8000 * mag * cosf(rad) * gain * Window_function::blackman(i / taps_f)));
}

HAL_TIM_Base_Stop(&htim7);
HAL_FMAC_FilterStop(&hfmac);

HAL_StatusTypeDef status = HAL_OK;
FMAC_FilterConfigTypeDef filter_config = {
    .InputBaseAddress = 0,
    .InputBufferSize = 120,
    .InputThreshold = FMAC_THRESHOLD_1,
    .CoeffBaseAddress = 120,
    .CoeffBufferSize = 120,
    .OutputBaseAddress = 240,
    .OutputBufferSize = 16,
    .OutputThreshold = FMAC_THRESHOLD_1,
    .pCoeffB = coeff_buff,
    .CoeffBSize = static_cast<uint8_t>(taps),
    .InputAccess = FMAC_BUFFER_ACCESS_DMA,
    .OutputAccess = FMAC_BUFFER_ACCESS_POLLING,
    .Clip = FMAC_CLIP_ENABLED,
    .Filter = FMAC_FUNC_CONVO_FIR,
    .P = static_cast<uint8_t>(taps),
};

status = HAL_FMAC_FilterConfig(&hfmac, &filter_config);
if (status)
{
    printf("FMAC config error %d %lu\n", status, hfmac.ErrorCode);
    break;
}

status = HAL_FMAC_FilterStart(&hfmac, nullptr, nullptr);
if (status)
{
    printf("FMAC start error %d %lu\n", status, hfmac.ErrorCode);
    break;
}

osDelay(2); // ノーウエイトで通すとFMACの計算が終わる前に読み出されてアンダーフローが発生する場合がある
HAL_TIM_Base_Start(&htim7);
```

### Sinc補償フィルタの場合

```c++
// do { } while (false);で囲ってある
constexpr uint16_t DFT_N = 128;

float bandwidth = 0,
        gain_dB = 0;
int32_t coeff_N = 0;
if (3 != sscanf(linebuff, "filter sinc %f %f %ld", &bandwidth, &gain_dB, &coeff_N))
{
    printf("command format error\n");
    break;
}

if (bandwidth < 0 || 0.5 < bandwidth ||
    coeff_N < 1 || 110 < coeff_N)
{
    printf("parameter error\n");
    break;
}

// 動的に確保したメモリは開放しない(リーク・フラグメンテーション対策)
static std::complex<float32_t> *DFT_buff = nullptr;
static int16_t *coeff_buff = nullptr;

if (!DFT_buff)
{
    DFT_buff = static_cast<std::complex<float32_t> *>(pvPortMalloc(sizeof(std::complex<float32_t>) * DFT_N));
}
if (!coeff_buff)
{
    coeff_buff = static_cast<int16_t *>(pvPortMalloc(sizeof(int16_t) * DFT_N));
}

if (!DFT_buff ||
    !coeff_buff)
{
    printf("malloc error\n");
    break;
}

const float32_t gain = powf(10, gain_dB / 20.0f);

for (int32_t i = 0; i < DFT_N; i++)
{
    const int32_t j = i < DFT_N / 2 ? i : i - DFT_N;
    const float32_t a = j / static_cast<float32_t>(DFT_N);
    const float32_t sinc = a == 0 ? 1 : sinf(a * pi) / (a * pi);
    DFT_buff[i] = fabsf(a) < bandwidth ? gain / sinc : 0;
}

arm_cfft_f32(&arm_cfft_sR_f32_len128, reinterpret_cast<float32_t *>(DFT_buff), 1, 1);
static_assert(DFT_N == 128);

class Window_function
{
public:
    static float32_t blackman(float32_t x)
    {
        return x < 0 || 1 < x ? 0 : 0.42f - 0.5f * cosf(2 * pi * x) + 0.08f * cosf(4 * pi * x);
    }
};

const float32_t coeff_N_f = coeff_N;
// 時間軸の前後を入れ替えて中央を切り出してcoeffに入れる
for (int32_t i = 0, j = coeff_N / 2, k = DFT_N - coeff_N / 2; i < coeff_N / 2; i++, j++, k++)
{
    coeff_buff[i] = static_cast<int16_t>(
        roundf(0x8000 * DFT_buff[k].real() * Window_function::blackman(i / coeff_N_f)));
    coeff_buff[j] = static_cast<int16_t>(
        roundf(0x8000 * DFT_buff[i].real() * Window_function::blackman(j / coeff_N_f)));
}

HAL_TIM_Base_Stop(&htim7);
HAL_FMAC_FilterStop(&hfmac);

HAL_StatusTypeDef status = HAL_OK;
FMAC_FilterConfigTypeDef filter_config = {
    .InputBaseAddress = 0,
    .InputBufferSize = 120,
    .InputThreshold = FMAC_THRESHOLD_1,
    .CoeffBaseAddress = 120,
    .CoeffBufferSize = 120,
    .OutputBaseAddress = 240,
    .OutputBufferSize = 16,
    .OutputThreshold = FMAC_THRESHOLD_1,
    .pCoeffB = coeff_buff,
    .CoeffBSize = static_cast<uint8_t>(coeff_N),
    .InputAccess = FMAC_BUFFER_ACCESS_DMA,
    .OutputAccess = FMAC_BUFFER_ACCESS_POLLING,
    .Clip = FMAC_CLIP_ENABLED,
    .Filter = FMAC_FUNC_CONVO_FIR,
    .P = static_cast<uint8_t>(coeff_N),
};

status = HAL_FMAC_FilterConfig(&hfmac, &filter_config);
if (status)
{
    printf("FMAC config error %d %lu\n", status, hfmac.ErrorCode);
    break;
}

status = HAL_FMAC_FilterStart(&hfmac, nullptr, nullptr);
if (status)
{
    printf("FMAC start error %d %lu\n", status, hfmac.ErrorCode);
    break;
}

osDelay(2); // ノーウエイトで通すとFMACの計算が終わる前に読み出されてアンダーフローが発生する場合がある
HAL_TIM_Base_Start(&htim7);
```