1
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

@yoshihiro0709(後藤 義宏)

nRF Connect SDK:FFT

1
Last updated at Posted at 2023-05-29

え、FFT？

今まではZephyrやnRF Connectの新しい機能や使い方などがメインでしたが、やや毛色の違うネタとしてFFTを扱ってみたいと思います。

浮動小数点に鳥肌が立つ

僕は非常に古くからの、それこそX68000時代のMC68000からのCPUユーザーなので、マイコン(CPU)に小数点演算をやらせるということ自体に未だ拒否反応があるのは事実です(笑)

昔のCPUの小数点演算はメチャクチャ遅かったのです。扱えないわけではなかったのですが。

今時の若い人(昔のCPUを知らない人)には絶対に分からない感覚なのだろうな、と思います。

nRF SDKにはサンプルがある

おっと、タイトルにnRF Connect SDKとあるのになぜかnRF SDKが出てきました(笑)。だって、nRF Connect SDKにはFFTのサンプルがないんですもの…。いや、もちろん作りますけど、とりあえず簡単に動かせるサンプルがあるnRF SDKからいきましょうよ！

とりあえず動かしてみよう

examples > peripheral > fpu_fftにサンプルプロジェクトがあります。その下はさらに細分化されていて、簡単に言うとIRQを使うか使わないか、BLE(最優先)を使うか使わないかで4つに分かれていますが、正直どれでもよいです(笑)。
実行するとCOMポートにヒストグラムっぽいテキストデータが何度も描画されますが、このFFTの結果が正しいのかどうかは僕には分かりません…っていうか周波数軸が分かんないんだけど…w
まあそれはともかく、なんか知らないけれど動いているっぽいのでどんどん進みます。

ランダムはやりにくいので

オリジナルのサンプルプロジェクトはランダムで周波数を決め、それに基づいてFFTを実行していますが、それだと本当に同じ値になるかどうかの検証が難しいのでランダムを使っている周波数の部分を固定値にします。
当たり前ですが固定値だと何度やっても同じ結果になるので無限ループも削っています。

main.c
/**
 * @brief Function for application main entry.
 */
int main(void)
{
    uint32_t  err_code;
    bool      noise = false;
    float32_t sine_freq;

    err_code = NRF_LOG_INIT(NULL);
    APP_ERROR_CHECK(err_code);

    NRF_LOG_DEFAULT_BACKENDS_INIT();

    // Welcome message.
    NRF_LOG_INFO("FPU FFT example started.");
    NRF_LOG_RAW_INFO("This is FPU usage example with FFT calculation and drawing.\r\n");

#ifdef FPU_INTERRUPT_MODE
    // Enable FPU interrupt
    NVIC_SetPriority(FPU_IRQn, APP_IRQ_PRIORITY_LOWEST);
    NVIC_ClearPendingIRQ(FPU_IRQn);
    NVIC_EnableIRQ(FPU_IRQn);
#endif

    // Generate new samples.
    // Sine wave frequency must be positive number so cast rand to unsigned number.
    // stdlib using with Segger Embedded Studio uses rand() configured to return values up to
    // 32767 instead of 0x7fffffff configured in other compilers. It causes that generated sine
    // frequency is smaller, but example works in the same way.
    sine_freq = (1000 % ((uint32_t)(SINE_WAVE_FREQ_MAX * SIGNALS_RESOLUTION))) / SIGNALS_RESOLUTION;
    fft_generate_samples(m_fft_input_f32,
                         FFT_TEST_COMP_SAMPLES_LEN,
                         FFT_TEST_SAMPLE_FREQ_HZ,
                         sine_freq, noise);

    // Process generated data. 64 pairs of complex data (real, img). It is important to use
    // proper arm_cfft_sR_f32 structure associated with input/output data length.
    // For example:
    //  - 128 numbers in input array (64 complex pairs of samples) -> 64 output bins power data -> &arm_cfft_sR_f32_len64.
    //  - 256 numbers in input array (128 complex pairs of samples) -> 128 output bins power data -> &arm_cfft_sR_f32_len128.
    fft_process(m_fft_input_f32,
                &arm_cfft_sR_f32_len64,
                m_fft_output_f32,
                FFT_TEST_OUT_SAMPLES_LEN);

    // Draw FFT bin power chart.
    draw_fft_header(sine_freq, noise);
    draw_fft_data(m_fft_output_f32, FFT_TEST_OUT_SAMPLES_LEN, GRAPH_WINDOW_HEIGHT);

    NRF_LOG_FLUSH();
}

出力コードを書くのがめんどくさいのでデバッグで変数の中身を覗いて貼り付けます(笑)
image.png
image.png

Pythonで同じものを実装する

当たり前の話ですが、PCであろうがマイコンであろうが同じデータをFFTに入力したら同じ出力結果になるはずです。それを確かめるためにPythonで同じデータを生成してFFTを実行してみます。

fft.py
import math
import numpy as np

# Define constant(s)
SAMPLING_SIZE = 128
SAMPLING_FREQUENCY = 44100
SINE_FREQUENCY = 10

# FFTデータの生成
# nRF52のサンプルと同等になるように生成する
fft_input_data = [];
for i in range(int(SAMPLING_SIZE / 2)):
    # Real part.
    rad = (SINE_FREQUENCY * 2.0 * math.pi * i) / SAMPLING_FREQUENCY
    print("radian = " + str(rad))
    fft_input_data.append(math.sin(rad))
    # Image part.
    fft_input_data.append(0j)

# FFT
fft_output_data = np.fft.fft(fft_input_data)

# 計算結果の表示
for i in range(len(fft_output_data)):
    print("fft_output_data(" + str(i) + ") = " + str(abs(fft_output_data[i])))

これを実行してみると…

実行してみた結果がこれです。全部は多いので一部抜粋。

fft_output_data(0) = 2.8703547309877373
fft_output_data(1) = 0.9281097203899827
fft_output_data(2) = 0.4645410409839072
fft_output_data(3) = 0.3103080754130883
fft_output_data(4) = 0.23338515235259694
・・・
fft_output_data(59) = 0.1873855250140315
fft_output_data(60) = 0.23338515235259694
fft_output_data(61) = 0.3103080754130884
fft_output_data(62) = 0.4645410409839072
fft_output_data(63) = 0.9281097203899826

ちゃんと一致していることが確認できました。

nRF Connect SDK(以下NCS)に実装

さて、いよいよNCSでの実装ですが、先ほども書いたように残念ながらNCSにはサンプルプロジェクトは付いていません。
まずARMの浮動小数点演算ですが、CMSIS(しーえむしす)というライブラリで提供されているようです。

基本的にはここにある関数を使ってFFTを実装することになります。って言われても！！！どうやって書けばいいんだよ！！！
ということで実装してみましたが、nRF SDKよりはやっていることは分かりやすいです。

main.c
/*
 * Copyright (c) 2012-2014 Wind River Systems, Inc.
 *
 * SPDX-License-Identifier: Apache-2.0
 */

#include <arm_math.h>
#include <stdio.h>
#include <zephyr/kernel.h>
#include <zephyr/logging/log.h>

LOG_MODULE_REGISTER(main, LOG_LEVEL_INF);

// Define constant(s)
#define SAMPLING_SIZE		128
#define SAMPLING_FREQUENCY	44100
#define SINE_FREQUENCY		10
#define FORWARD_TRANSFORM	0
#define INVERSE_TRANSFORM	1
#define BIT_REVERSE_NO		0
#define BIT_REVERSE_YES		1

void main(void)
{
	arm_cfft_instance_f32 cfft_f32;
	arm_status status = ARM_MATH_SUCCESS;
	float32_t fft_input_data[SAMPLING_SIZE];
	float32_t fft_output_data[SAMPLING_SIZE / 2];
	float32_t rad;
	float32_t max_value;
	uint32_t max_value_index = 0;
	char string_buffer[256];

	// Make FFT data
	memset(fft_input_data, 0, sizeof(fft_input_data));
	for (int i = 0; i < SAMPLING_SIZE; i+=2) {
		rad = (SINE_FREQUENCY * 2.0 * PI * (i / 2)) / SAMPLING_FREQUENCY;
		fft_input_data[i] = sin(rad);
		// For confirmation
		// sprintf(string_buffer, "input_data(%d) = %f", i, fft_input_data[i]);
		// LOG_INF("%s", string_buffer);
		// k_sleep(K_MSEC(10));
	}

	// Arm Complex FFT
	status = arm_cfft_init_f32(&cfft_f32, (SAMPLING_SIZE / 2));
	if (status == ARM_MATH_SUCCESS) {
		LOG_INF("Initialize done.");
	} else {
		LOG_ERR("Failed to initialize.");
		return;
	}
	arm_cfft_f32(&cfft_f32, fft_input_data, FORWARD_TRANSFORM, BIT_REVERSE_YES);
	arm_cmplx_mag_f32(fft_input_data, fft_output_data, (SAMPLING_SIZE / 2));
	arm_max_f32(fft_output_data, (SAMPLING_SIZE / 2), &max_value, &max_value_index);

	for (int i = 0; i < SAMPLING_SIZE / 2; i++) {
		sprintf(string_buffer, "output_data(%d) = %f", i, fft_output_data[i]);
		LOG_INF("%s", string_buffer);
		k_sleep(K_MSEC(10));
	}
	sprintf(string_buffer, "Max value is %f", max_value);
	LOG_INF("%s", string_buffer);
	LOG_WRN("waiting for reboot...");
}
prj.conf
CONFIG_MAIN_STACK_SIZE=2048
CONFIG_CMSIS_DSP=y
CONFIG_CMSIS_DSP_COMPLEXMATH=y
CONFIG_CMSIS_DSP_STATISTICS=y
CONFIG_CMSIS_DSP_TRANSFORM=y
CONFIG_NEWLIB_LIBC=y
CONFIG_NEWLIB_LIBC_MIN_REQUIRED_HEAP_SIZE=2048
CONFIG_NEWLIB_LIBC_FLOAT_PRINTF=y
CONFIG_DEBUG_THREAD_INFO=y
CONFIG_LOG=y
CONFIG_DEBUG_OPTIMIZATIONS=y

メインメモリを2048 Bytes以上取らないとメモリ不足でFatal Errorを引き起こします。

実行して出力を見てみると…
image.png
おお、これも前述のnRF SDKやPythonの計算と一致しています。

NCSで実装したほうがよさそう？

実際にFFTをプロジェクトレベルで使ったわけではないのでこれは推測になりますが、やはりRTOSでスイッチできるのとFPU自体は非同期で動作している(と思われる)ことを考えるとNCSで実装するほうがいいんじゃないかと思います。
ここらへんは今後も検証していく予定です。

そりゃ必要だから動作検証しているわけですよ、もちろん…

余談

元々の素材となったnRF SDKのサンプルプロジェクトですが、冷静にデータを見てみると入力データで生成しているのが1波形の1/20くらいしかない(44100Hzなので1周期の2πまでプロットするのに2205データ必要)ので、FFTも何もあったものではないような気はします(笑)。
が、本題はそれぞれのFFTが一致するかの検証なので元データの質は問わないものとしましょう。

1
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?