RFFT in CMSIS DSP. Part 1.

RFFT in STM32 using CMSIS DSP. Part 1.

Python warm-up for illustration

I always wanted to use MCU for audio processing. And for my purposes, I need Discrete Fourier Transform(DFT), especially its fast version FFT. If you do not know what is this, read Wikipedia first:) FFT has a huge number examples of usage, for my case I want to build wavetable synthesizer.

Let’s take some real example. We have an audio signal, which is just a set of amplitudes, which you will take from your ADC. Samples are taken with some sampling frequency Fs. For preliminary illustration, we’ll use Python.

import numpy as np
import matplotlib.pyplot as plt

# Our wave settings
N_samples = 32
amplitude = 10
# Here we'll generate our sine wave
t = np.arange(N_samples)
x = amplitude * np.sin(2 * np.pi * t / N_samples)

# Plotting
fig, ax = plt.subplots()
ax.plot(t, x, 'bo-')
ax.set_xlabel('N sample')
ax.set_ylabel('Amplitude');
plt.show()

As a result, we’ll have this:

We have built a sine wave with 32 samples. A sine wave has real values. But FFT results are complex-valued.

from scipy import fftpack
#Two different functions for RFFT, I'll tell you difference later
X = fftpack.rfft(x)
Xc = np.fft.rfft(x) # Result of rfft
spectrum = np.abs(Xc) # Spectrum

fig, ax = plt.subplots()
ax.stem(spectrum)
ax.set_ylabel('Frequency Domain (Spectrum) Magnitude')
ax.set_xlabel('Wavenumber')
plt.grid(True)

print(Xc)

A result of FFT is:

[ 1.89078061e-15+0.00000000e+00j -2.43172857e-14-1.60000000e+02j
 -8.58766605e-15-9.97865745e-15j  2.44724107e-15-9.14741103e-15j
  1.27390193e-15-4.72797206e-15j  9.76953827e-16-8.90236315e-15j
  4.56341236e-15-4.94620943e-15j -7.85174677e-16-1.18395480e-14j
 -1.88397767e-15+3.77475828e-15j  8.99402212e-15-1.18395480e-14j
  1.10746579e-14+3.43377958e-15j  3.67917993e-15+1.75577788e-15j
  4.72810534e-15+5.04199055e-15j  2.20889268e-15+1.51073001e-15j
 -2.15181701e-15+5.50675892e-15j -3.00100364e-15+0.00000000e+00j
 -3.29665435e-16+0.00000000e+00j]

as you see, every value has two parts: real and imaginary. Real part stores the frequency of the sine wave. Imaginary stores info about phase. But to have spectra(dependence of sine wave amplitude vs wavenumber), we need np.abs function, which converts complex values to real amplitude values.

As we see, the first freq is zero. The first frequency is about DC offset. Our sine wave is symmetrical from zero value. The second one is our main frequency.
Also, you can notice, that the number of frequencies is only 17, but a number of samples is 32. This is because our FFT is real-valued and the number of frequencies is limited by Nyquist frequency, which is twice less then number of samples(that is how I understood this, correct me in comments if I am wrong!:)

Real example on STM32

As usual, we’ll use System Workbench for STM32 together with CubeMX for code generation. Our development board is NUCLEO-F401RE. MCU has DSP instructions and the floating point unit, and we want to use the full power of it. To use them, we need CMSIS DSP library, download it from ST site there.
Link to site
Unzip it somewhere.

Project setup

1. Start default project for NUCLEO-F401RE board. There is nothing to explain, we do not change anything. Just generate code for System Workbench IDE.
2. Now we need to modify our project directory tree structure.
   1. Make directory <project path>/Drivers/CMSIS/Lib
   2. Put there ARM and GCC directories from <downloaded cube dir>/Drivers/Lib
   3. Make directory DSP_Lib in the <project path>/Drivers/CMSIS/
   4. Copy <downloaded cube dir>/Drivers/CMSIS/DSP/Source to this directory.
   5. To directory <project path>/Drivers/CMSIS/Includes we need to add additional include files to enable math. Copy arm_common_tables.h, arm_const_structs.h, arm_math.h from <downloaded cube dir>/Drivers/CMSIS/DSP/Include to <project path>/Drivers/CMSIS/Include
3. In System workbench add new library search paths: Settings->Libraries->Library search path(-L). Add Drivers/CMSIS/Lib/ARM and Drivers/CMSIS/Lib/GCC directories. Settings access is through right-click on Project Explorer, then go “Properties”. Then choose C/C++ Build. There is Settings submenu.
4. Add library arm_cortexM4lf_math to the list of librariesSettings->Libraries->Libraries(-l) .
5. Add defined symbols ARM_MATH_CM4 and __FPU_PRESENT=1 (Settings->Preprocessor->Defined symbols)

Generate sine wave in python

Now let’s write some code. Let’s print sine wave from the previous part with python:

import numpy as np

N_samples = 32
N_harmonics = 1
amplitude = 10
t = np.arange(N_samples)

x = amplitude * np.sin(2 * np.pi * t / N_samples)

for sample in x:
    print(sample, ',', sep='')

Output array is:

0.0,
1.9509032201612824,
3.826834323650898,
5.555702330196022,
7.071067811865475,
8.314696123025453,
9.238795325112868,
9.807852804032304,
10.0,
9.807852804032304,
9.238795325112868,
8.314696123025454,
7.0710678118654755,
5.555702330196022,
3.826834323650899,
1.9509032201612861,
1.2246467991473533e-15,
-1.9509032201612837,
-3.8268343236508966,
-5.55570233019602,
-7.071067811865475,
-8.314696123025453,
-9.238795325112864,
-9.807852804032303,
-10.0,
-9.807852804032304,
-9.238795325112866,
-8.314696123025454,
-7.071067811865477,
-5.555702330196022,
-3.826834323650904,
-1.9509032201612873,

Write main code

In main.c let’s define global variables:

float32_t testInput[32] =
{
    0.0,
    1.9509032201612824,
    3.826834323650898,
    5.555702330196022,
    7.071067811865475,
    8.314696123025453,
    9.238795325112868,
    9.807852804032304,
    10.0,
    9.807852804032304,
    9.238795325112868,
    8.314696123025454,
    7.0710678118654755,
    5.555702330196022,
    3.826834323650899,
    1.9509032201612861,
    1.2246467991473533e-15,
    -1.9509032201612837,
    -3.8268343236508966,
    -5.55570233019602,
    -7.071067811865475,
    -8.314696123025453,
    -9.238795325112864,
    -9.807852804032303,
    -10.0,
    -9.807852804032304,
    -9.238795325112866,
    -8.314696123025454,
    -7.071067811865477,
    -5.555702330196022,
    -3.826834323650904,
    -1.9509032201612873
};
float32_t testOutput[32];

in main function I’ve deleted unnecessary comments. CMSIS DSP functions are arm_rfft_fast_init_f32 - for initialization FFT structure. arm_rfft_fast_f32 - RFFT itself. It is interesting, it looks like arm_rfft_fast_f32 performs all the changes in place, and it will modify our testInput array. RFFT result will be put in testOuptut array. Also, notice that you need to specify FFT size, in our case, this is 32. All additional documentation you can find there

int main(void)
{
  HAL_Init();
  SystemClock_Config();
  MX_GPIO_Init();
  MX_USART2_UART_Init();
  /* Our code */
  arm_rfft_fast_instance_f32 S; // RFFT instance. 
  arm_status status;
  status = arm_rfft_fast_init_f32(&S, 32);
  arm_rfft_fast_f32(&S, testInput, testOutput, 0);

  if ( status != ARM_MATH_SUCCESS)
  {
    /* Add your own assert */
  }
  printf("N, Freq\r\n");
  for (int i = 0; i < sizeof(testOutput)/sizeof(testOutput[0]); i++)
  {
    printf("%d, %f\r\n", i, testOutput[i]);
  }
  /* Infinite loop */
  /* USER CODE BEGIN WHILE */
  while (1)
  {
  }
}

Miscellaneous - redirect printf to UART

We use printf function. I want to print to UART. To “redirect” printf to UART we need to define our own syscall _write somewhere in our main.c:

int _write(int file, char *ptr, int len)
{
  if ((file != STDOUT_FILENO) && (file != STDERR_FILENO))
  {
     errno = EBADF;
     return -1;
  }

  HAL_StatusTypeDef status =
     HAL_UART_Transmit(&huart2, (uint8_t*)ptr, len, 1000);

  return (status == HAL_OK ? len : 0);
}

Miscellaneous - allow float modifier for printf

By default, printf cannot print floats. To enable this feature we need to add linker flag -u _printf_float. Do this there Properties -> C/C++ Build -> Settings -> MCU GCC Linker -> Miscellaneous -> Linker flags

Compile, download, results

N, Freq
0, 0.000000
1, 0.000000
2, -0.000001
3, -159.999985
4, 0.000000
5, 0.000000
6, 0.000001
7, 0.000002
8, 0.000000
9, 0.000000
10, 0.000001
11, -0.000002
12, 0.000000
13, 0.000000
14, -0.000004
15, -0.000000
16, 0.000000
17, 0.000000
18, 0.000000
19, 0.000002
20, 0.000000
21, 0.000000
22, 0.000001
23, -0.000000
24, 0.000000
25, 0.000000
26, -0.000003
27, 0.000002
28, 0.000000
29, 0.000000
30, 0.000001
31, 0.000012

Ok, what do we need to notice? C has its own type for complex variables, but a lot of libraries still use own ways to represent complex results. That’s why we have more values, compared to python. Values are packed like this:

X = { real[0], imag[0], real[1], imag[1], real[2], imag[2] ... real[(N/2)-1], imag[(N/2)-1 }

Documentation says:

It happens that the first complex number (real[0], imag[0]) is actually all real. real[0] represents the DC offset, and imag[0] should be 0. (real[1], imag[1]) is the fundamental frequency, (real[2], imag[2]) is the first harmonic and so on.

The real FFT functions pack the frequency domain data in this fashion.
The forward transform outputs the data in this form and the inverse
transform expects input data in this form. The function always
performs the needed bitreversal so that the input and output data is
always in normal order. The functions support lengths of [32, 64, 128,
…, 4096] samples.

Similar results in python gives fftpack.rfft() function, but values are packed in some different way, but still without complex-number.

Epilogue

I am not an expert in DSP and this article is based only on practical experience I got. If you have any comments or improvements to this information you are welcome. Write to me directly on a.synytsyn@gmail.com or in comments to this article. Also, grammar fixes are welcome also:)

All code is available: here