テストコード
bool flag(true);
uint8_t is_complex(0);
size_t length(0);
float SPS(0);
float freq(0);
const uint8_t ifft_flag(0);
const uint8_t do_bit_reverse(1);
flag = flag && (4 == sscanf(command,
"fft %hhu %u %f %f",
&is_complex,
&length,
&SPS,
&freq));
if (flag)
{
printf("FFT %s %upoint %.1fSPS %.1fHz\n",
is_complex ? "complex" : "real",
length,
SPS,
freq);
}
if (is_complex)
{
const arm_cfft_instance_f32 *S(0);
if (flag)
{ // check length
switch (length)
{
case 16:
S = &arm_cfft_sR_f32_len16;
break;
case 32:
S = &arm_cfft_sR_f32_len32;
break;
case 64:
S = &arm_cfft_sR_f32_len64;
break;
case 128:
S = &arm_cfft_sR_f32_len128;
break;
case 256:
S = &arm_cfft_sR_f32_len256;
break;
case 512:
S = &arm_cfft_sR_f32_len512;
break;
case 1024:
S = &arm_cfft_sR_f32_len1024;
break;
case 2048:
S = &arm_cfft_sR_f32_len2048;
break;
case 4096:
S = &arm_cfft_sR_f32_len4096;
break;
default:
flag = false;
printf("length error\n");
break;
}
}
float *FFT_input(0);
if (flag)
{ // alloc
try
{
FFT_input = new float32_t[length * 2]();
}
catch (std::bad_alloc)
{
flag = false;
printf("bad alloc\n");
}
}
if (flag)
{ //generate wave
for (size_t i(0); i < length; i++)
{
const float phase(i * pi * 2 / SPS * freq);
FFT_input[i * 2 + 0] = cosf(phase);
FFT_input[i * 2 + 1] = sinf(phase);
}
}
float FFT_process_seconds(0);
if (flag)
{ // FFT
const uint32_t start(HAL_GetTick());
arm_cfft_f32(S, FFT_input, ifft_flag, do_bit_reverse);
const uint32_t stop(HAL_GetTick());
FFT_process_seconds = (stop - start) * 0.001f;
}
if (flag)
{ // dump result
printf(" %.3fsec\n\n re im\n", FFT_process_seconds);
for (size_t i(0); i < length; i++)
{
const int j(i - length / 2);
const int k(i < length / 2 ? (i + length / 2) : (i - length / 2));
printf("%f %f %f\n",
j * SPS / length,
FFT_input[k * 2 + 0],
FFT_input[k * 2 + 1]);
}
}
delete[] FFT_input;
FFT_input = 0;
}
else
{
arm_rfft_fast_instance_f32 S;
if (flag)
{ // check length and initialize
flag = (ARM_MATH_SUCCESS == arm_rfft_fast_init_f32(&S, length));
if (!flag)
{
printf("init error\n");
}
}
float *FFT_input(0);
float *FFT_output(0);
if (flag)
{ // alloc
try
{
FFT_input = new float32_t[length]();
FFT_output = new float32_t[length]();
}
catch (std::bad_alloc)
{
flag = false;
printf("bad alloc\n");
}
}
if (flag)
{ //generate wave
for (size_t i(0); i < length; i++)
{
const float phase(i * pi * 2 / SPS * freq);
FFT_input[i] = cosf(phase);
}
}
float FFT_process_seconds(0);
if (flag)
{ // FFT
arm_rfft_fast_f32(&S, FFT_input, FFT_output, ifft_flag);
}
if (flag)
{ // dump result
printf("%.3fsec\n\n re\n", FFT_process_seconds);
for (size_t i(0); i < length; i++)
{
printf("%f %f\n",
i * SPS / length / 2,
FFT_output[i]);
}
}
delete[] FFT_input;
FFT_input = 0;
delete[] FFT_output;
FFT_output = 0;
}
リンカでDrivers/CMSIS/Lib/GCC/libarm_cortexM4lf_math.aを読み込む必要がある。
ARM_MATH_CM4(コアに合わせて選択)をdefine定義した後、arm_math.hとarm_const_structs.hをインクルードする。あとnewのstd::bad_allocをcatchするのに<new>もインクルードする必要がある。
CFFTは複素数、RFFTは実数の処理で、それぞれ呼び方が違う。また引数のとり方も違い、CFFTはinputとoutputを同じ配列で行う。一方、RFFTはinputとoutputは別の配列を使う。
CFFTは負の周波数も解析できる。一方、RFFTは負の周波数は扱えないが、CFFTの分解能の2倍くらいあると思う。1024ポイントでSPS1kHzとした場合、CFFTは1024ポイントで-0.5kHzから+0.5kHzの範囲を取る。RFFTは1024ポイントで0Hzから0.5kHzの範囲を取る。
ちゃんと試してないけど、CFFTは16ポイントから4kポイントまで、RFFTは32ポイントから4kポイントまで解析できる。
CFFT/RFFTを行う上記コードを含めるとバイナリサイズは256kくらいになった。この部分をコメントアウトすると128k位になる。FFT解析を行うには128kほどのFlashが必要になる。もちろん、FFTを行う上でメモリ領域も必要になる。
const char *const commandに文字列を入れておけば、それを解析する。主にUART等からコマンドを受け取って実行する事を想定。もちろんsprintfで作ってもいいけど、そんな事するくらいなら直接呼べという話で。
"fft 1 4096 100 -10" というコマンドを実行すると、複素数で100SPS/-10Hzの波形を生成しFFT解析する。
"fft 0 4096 100 30" というコマンドを実行すると、実部のみで100SPS/30Hzの波形を生成しFFT解析する。
それぞれの結果は以下の通り。
CFFTは左端が-50Hz、右端が+50Hz、RFFTは左端が0Hz、右端が+50Hzになっている。どちらも結果は4096ポイント得られるので、RFFTのほうが倍の分解能になる。
4096ポイントの場合、168MHzのCortex-M4FでCFFTが4msec未満、RFFTが1msec未満だった。CFFTは若干の時間がかかるが、RFFTは恐ろしく早い。
ただ、どちらも4kポイントだと32KiBのRAMが必要になる。自由に使えるヒープ領域は100KiB程度なので、あまり大規模なものを作ろうとするとちょっと心もとない。***
とりあえずFFTが動くようになったので、パルス圧縮の相関処理とかもできるわけだが、時間的・空間的な余裕がかなり少なそう。
超音波距離計もちょっと飽きてきたし、気分転換に別のヤツで遊ぼうかな。
FFT/IFFTが動くならOFDM変調とかできる。
0 件のコメント:
コメントを投稿