Занимательная электроника

В этом примере разберём, как выполнить разложение сигнала на частотные составляющие с помощью дискретного преобразования Фурье на микроконтроллере STM32F411.

В предыдущем примере мы использовали цифровой фильтр, чтобы убрать помехи из сигнала. Здесь фильтр уже не нужен. Теперь задача другая: взять шумный сигнал и посмотреть, из каких частот он состоит.

Для этого мы программно создадим сигнал, в котором есть три синусоиды:

200 Гц — полезный сигнал
600 Гц — первая помеха
900 Гц — вторая помеха

После этого выполним разложение Фурье и выведем спектр по UART. На графике должны появиться три основных пика: на 200 Гц, 600 Гц и 900 Гц.

---

1. Что такое разложение Фурье

Любой сложный сигнал можно представить как сумму простых синусоид разных частот, амплитуд и фаз. Именно это и показывает преобразование Фурье.

Если посмотреть на сигнал во временной области, мы видим только изменение амплитуды во времени. Например, сигнал может выглядеть как сложная кривая с рябью и помехами.

Но если выполнить преобразование Фурье, мы переходим в частотную область. Там уже видно, какие частоты присутствуют внутри сигнала.

Проще говоря:

Временная область:
смотрим, как сигнал меняется во времени.

Частотная область:
смотрим, из каких частот состоит сигнал.

---

2. Какой сигнал будем анализировать

В программе мы формируем шумный сигнал из трёх составляющих:

input = clean + noise;

Где clean — это полезный сигнал 200 Гц:

clean = sinf(2.0f * PI_F * 200.0f * t);

А noise — это две дополнительные синусоиды, которые имитируют помеху:

noise = 0.35f * sinf(2.0f * PI_F * 600.0f * t)
      \+ 0.25f * sinf(2.0f * PI_F * 900.0f * t);

Итоговый сигнал получается таким:

input = 1.00 · sin(200 Гц)
      \+ 0.35 · sin(600 Гц)
      \+ 0.25 · sin(900 Гц)

То есть в сигнале есть:

200 Гц с амплитудой 1.00
600 Гц с амплитудой 0.35
900 Гц с амплитудой 0.25

После преобразования Фурье мы должны увидеть эти частоты в спектре.

---

3. Частота дискретизации и количество отсчётов

Сигнал внутри STM32 создаётся не непрерывно, а в виде отдельных отсчётов.

В программе используется частота дискретизации:

Fs = 4000 Гц

Это значит, что между соседними отсчётами проходит:

1 / 4000 = 0.00025 с = 250 мкс

Для анализа берём 400 отсчётов:

N = 400

Разрешение по частоте определяется так:

Δf = Fs / N

Подставим наши значения:

Δf = 4000 / 400 = 10 Гц

Это значит, что спектр будет рассчитываться с шагом 10 Гц:

0 Гц
10 Гц
20 Гц
30 Гц
...
1000 Гц

Это удобно, потому что наши частоты точно попадают в эти точки:

200 Гц / 10 Гц = 20
600 Гц / 10 Гц = 60
900 Гц / 10 Гц = 90

То есть пики спектра должны появиться в ячейках 20, 60 и 90.

---

4. Формула дискретного преобразования Фурье

Так как сигнал у нас цифровой, мы используем не непрерывный ряд Фурье, а дискретное преобразование Фурье — DFT.

Пусть у нас есть массив отсчётов:

x[0], x[1], x[2], ..., x[N-1]

Для каждой частотной ячейки k считается комплексное значение:

X[k] = Σ x[n] · e^(-j · 2π · k · n / N)

В программе мы не используем комплексные числа напрямую. Вместо этого отдельно считаем действительную и мнимую части.

Действительная часть:

Re[k] = Σ x[n] · cos(2π · k · n / N)

Мнимая часть:

Im[k] = -Σ x[n] · sin(2π · k · n / N)

После этого амплитуда частоты считается так:

Amplitude[k] = (2 / N) · sqrt(Re[k]^2 + Im[k]^2)

А сама частота для ячейки k:

f[k] = k · Fs / N

В нашем примере:

f[k] = k · 4000 / 400
f[k] = k · 10 Гц

---

5. Почему появляется пик на нужной частоте

DFT работает так: для каждой проверяемой частоты программа умножает входной сигнал на синус и косинус этой частоты.

Если такая частота действительно есть в сигнале, произведения складываются и дают большую сумму.

Если такой частоты в сигнале нет, положительные и отрицательные участки взаимно компенсируются, и сумма получается близкой к нулю.

Поэтому в спектре мы должны увидеть:

200 Гц → большой пик
600 Гц → меньший пик
900 Гц → ещё один пик

Так как амплитуды сигналов заранее заданы как 1.00, 0.35 и 0.25, результат должен быть примерно таким:

200 Гц → 1000
600 Гц → 350
900 Гц → 250

Значения выводятся в формате x1000, чтобы не использовать печать чисел с плавающей точкой через printf.

---

6. Генерация шумного сигнала

Сначала программа формирует массив inputSignal из 400 отсчётов.

static float inputSignal[DFT_N];

Функция генерации сигнала:

static void Generate_NoisySignal(void)
{
  for (uint16_t n = 0; n < DFT_N; n++)
  {
    float t;
    float clean;
    float noise;
    float input;

    t = (float)n / FS_HZ;

    clean = sinf(2.0f * PI_F * 200.0f * t);

    noise = 0.35f * sinf(2.0f * PI_F * 600.0f * t)
          \+ 0.25f * sinf(2.0f * PI_F * 900.0f * t);

    input = clean + noise;

    inputSignal[n] = input;
  }
}

В этой функции:

n — номер текущего отсчёта
t — текущее время
clean — полезный сигнал 200 Гц
noise — сумма помех 600 Гц и 900 Гц
input — итоговый шумный сигнал

Именно массив inputSignal потом будет разложен по Фурье.

---

7. Расчёт спектра

Функция DFT_SendSpectrumToUart() проходит по частотам от 0 до 1000 Гц и для каждой частоты считает амплитуду.

for (uint16_t k = 0; k <= DFT_MAX_BIN; k++)
{
  float real = 0.0f;
  float imag = 0.0f;

  for (uint16_t n = 0; n < DFT_N; n++)
  {
    float angle;

    angle = 2.0f * PI_F * (float)k * (float)n / (float)DFT_N;

    real += inputSignal[n] * cosf(angle);
    imag -= inputSignal[n] * sinf(angle);
  }
}

Здесь:

k — номер частотной ячейки
n — номер отсчёта сигнала
real — действительная часть X[k]
imag — мнимая часть X[k]

После расчёта действительной и мнимой частей находится амплитуда:

amplitude = (2.0f / (float)DFT_N) * sqrtf((real * real) + (imag * imag));

Затем амплитуда переводится в формат x1000:

amplitude_x1000 = (uint32_t)(amplitude * 1000.0f + 0.5f);

Частота для текущей ячейки:

freq_hz = (uint32_t)k * 10u;

И строка отправляется по UART:

snprintf(line, sizeof(line), "%lu,%lu\r\n",
         (unsigned long)freq_hz,
         (unsigned long)amplitude_x1000);

UART_SendString(line);

---

8. Что выводится по UART

Программа выводит таблицу:

freq_hz,amplitude_x1000

Пример вывода:

0,0
10,0
20,0
30,0
...
200,1000
...
600,350
...
900,250
...
1000,0
END_SPECTRUM

Для построения графика:

ось X — freq_hz
ось Y — amplitude_x1000

В результате на экране должны появиться три основных пика:

200 Гц
600 Гц
900 Гц

Это и есть частотный состав шумного сигнала.

---

9. Почему используется DFT, а не FFT

В этой программе используется прямой расчёт DFT. Он медленнее, чем FFT, но намного проще для понимания.

Для каждого значения частоты программа вручную считает суммы синуса и косинуса. Поэтому хорошо видно, как именно работает разложение сигнала на частоты.

FFT работает быстрее, но внутри алгоритм сложнее. Когда понятен принцип DFT, можно переходить к библиотеке CMSIS-DSP и использовать готовые функции FFT, например:

arm_rfft_fast_f32()

Но для учебного проекта лучше начать именно с прямого DFT.

---

10. Полный код программы

/* USER CODE BEGIN Header */
/**
  ******************************************************************************
  * @file           : main.c
  * @brief          : STM32F411 Fourier decomposition test
  ******************************************************************************
  */
/* USER CODE END Header */

/* Includes ------------------------------------------------------------------*/
#include "main.h"

/* Private includes ----------------------------------------------------------*/
/* USER CODE BEGIN Includes */
#include <math.h>
#include <stdio.h>
#include <string.h>
/* USER CODE END Includes */

/* Private define ------------------------------------------------------------*/
/* USER CODE BEGIN PD */

#define FS_HZ        4000.0f
#define PI_F         3.14159265358979323846f

/*
   N = 400 samples.
   Frequency resolution:
   Fs / N = 4000 / 400 = 10 Hz.
*/
#define DFT_N        400

/*
   We calculate spectrum from 0 Hz to 1000 Hz.
   Step = 10 Hz.
   1000 Hz / 10 Hz = 100 bins.
*/
#define DFT_MAX_BIN  100

/* USER CODE END PD */

/* Private variables ---------------------------------------------------------*/
TIM_HandleTypeDef htim2;
UART_HandleTypeDef huart2;

/* USER CODE BEGIN PV */

static float inputSignal[DFT_N];

/* USER CODE END PV */

/* Private function prototypes -----------------------------------------------*/
void SystemClock_Config(void);
static void MX_GPIO_Init(void);
static void MX_USART2_UART_Init(void);

/* USER CODE BEGIN PFP */

static void UART_SendString(const char *s);
static void Generate_NoisySignal(void);
static void DFT_SendSpectrumToUart(void);

/* USER CODE END PFP */

/* Private user code ---------------------------------------------------------*/
/* USER CODE BEGIN 0 */

static void UART_SendString(const char *s)
{
  HAL_UART_Transmit(&huart2, (uint8_t *)s, strlen(s), HAL_MAX_DELAY);
}

static void Generate_NoisySignal(void)
{
  for (uint16_t n = 0; n < DFT_N; n++)
  {
    float t;
    float clean;
    float noise;
    float input;

    t = (float)n / FS_HZ;

    /*
       Useful signal: 200 Hz.
    */
    clean = sinf(2.0f * PI_F * 200.0f * t);

    /*
       Noise components: 600 Hz and 900 Hz.
    */
    noise = 0.35f * sinf(2.0f * PI_F * 600.0f * t)
          \+ 0.25f * sinf(2.0f * PI_F * 900.0f * t);

    /*
       Noisy signal for Fourier decomposition.
    */
    input = clean + noise;

    inputSignal[n] = input;
  }
}

static void DFT_SendSpectrumToUart(void)
{
  char line[80];

  UART_SendString("\r\n");
  UART_SendString("Fourier spectrum of noisy signal\r\n");
  UART_SendString("Fs=4000Hz, N=400, step=10Hz\r\n");
  UART_SendString("Input signal: 200Hz + 600Hz + 900Hz\r\n");
  UART_SendString("freq_hz,amplitude_x1000\r\n");

  for (uint16_t k = 0; k <= DFT_MAX_BIN; k++)
  {
    float real = 0.0f;
    float imag = 0.0f;
    float amplitude;
    uint32_t amplitude_x1000;
    uint32_t freq_hz;

    for (uint16_t n = 0; n < DFT_N; n++)
    {
      float angle;

      angle = 2.0f * PI_F * (float)k * (float)n / (float)DFT_N;

      /*
         DFT:
         X[k] = sum x[n] * e^(-j*2*pi*k*n/N)

         real part = sum x[n] * cos(angle)
         imag part = -sum x[n] * sin(angle)
      */
      real += inputSignal[n] * cosf(angle);
      imag -= inputSignal[n] * sinf(angle);
    }

    /*
       For normal sinus amplitude:
       amplitude = 2 / N * abs(X[k])

       For k = 0 this doubling is not correct, but DC is not important
       in this example because our signal has no constant component.
    */
    amplitude = (2.0f / (float)DFT_N) * sqrtf((real * real) + (imag * imag));

    amplitude_x1000 = (uint32_t)(amplitude * 1000.0f + 0.5f);
    freq_hz = (uint32_t)k * 10u;

    snprintf(line, sizeof(line), "%lu,%lu\r\n",
             (unsigned long)freq_hz,
             (unsigned long)amplitude_x1000);

    UART_SendString(line);
  }

  UART_SendString("END_SPECTRUM\r\n");
}

/* USER CODE END 0 */

/**
  * @brief  The application entry point.
  * @retval int
  */
int main(void)
{
  HAL_Init();

  SystemClock_Config();

  MX_GPIO_Init();
  MX_USART2_UART_Init();

  /* USER CODE BEGIN 2 */

  HAL_Delay(500);

  UART_SendString("\r\n");
  UART_SendString("STM32F411 Fourier decomposition\r\n");

  /* USER CODE END 2 */

  while (1)
  {
    /*
       1. Generate noisy signal.
       2. Decompose it into Fourier spectrum.
       3. Send spectrum to UART.
    */
    Generate_NoisySignal();
    DFT_SendSpectrumToUart();

    HAL_Delay(3000);
  }
}

/**
  * @brief System Clock Configuration
  * @retval None
  */
void SystemClock_Config(void)
{
  RCC_OscInitTypeDef RCC_OscInitStruct = {0};
  RCC_ClkInitTypeDef RCC_ClkInitStruct = {0};

  __HAL_RCC_PWR_CLK_ENABLE();
  __HAL_PWR_VOLTAGESCALING_CONFIG(PWR_REGULATOR_VOLTAGE_SCALE1);

  RCC_OscInitStruct.OscillatorType = RCC_OSCILLATORTYPE_HSI;
  RCC_OscInitStruct.HSIState = RCC_HSI_ON;
  RCC_OscInitStruct.HSICalibrationValue = RCC_HSICALIBRATION_DEFAULT;
  RCC_OscInitStruct.PLL.PLLState = RCC_PLL_ON;
  RCC_OscInitStruct.PLL.PLLSource = RCC_PLLSOURCE_HSI;
  RCC_OscInitStruct.PLL.PLLM = 16;
  RCC_OscInitStruct.PLL.PLLN = 336;
  RCC_OscInitStruct.PLL.PLLP = RCC_PLLP_DIV4;
  RCC_OscInitStruct.PLL.PLLQ = 4;

  if (HAL_RCC_OscConfig(&RCC_OscInitStruct) != HAL_OK)
  {
    Error_Handler();
  }

  RCC_ClkInitStruct.ClockType = RCC_CLOCKTYPE_HCLK | RCC_CLOCKTYPE_SYSCLK
                              | RCC_CLOCKTYPE_PCLK1 | RCC_CLOCKTYPE_PCLK2;
  RCC_ClkInitStruct.SYSCLKSource = RCC_SYSCLKSOURCE_PLLCLK;
  RCC_ClkInitStruct.AHBCLKDivider = RCC_SYSCLK_DIV1;
  RCC_ClkInitStruct.APB1CLKDivider = RCC_HCLK_DIV2;
  RCC_ClkInitStruct.APB2CLKDivider = RCC_HCLK_DIV1;

  if (HAL_RCC_ClockConfig(&RCC_ClkInitStruct, FLASH_LATENCY_2) != HAL_OK)
  {
    Error_Handler();
  }
}

/**
  * @brief USART2 Initialization Function
  * @param None
  * @retval None
  */
static void MX_USART2_UART_Init(void)
{
  huart2.Instance = USART2;
  huart2.Init.BaudRate = 115200;
  huart2.Init.WordLength = UART_WORDLENGTH_8B;
  huart2.Init.StopBits = UART_STOPBITS_1;
  huart2.Init.Parity = UART_PARITY_NONE;
  huart2.Init.Mode = UART_MODE_TX_RX;
  huart2.Init.HwFlowCtl = UART_HWCONTROL_NONE;
  huart2.Init.OverSampling = UART_OVERSAMPLING_16;

  if (HAL_UART_Init(&huart2) != HAL_OK)
  {
    Error_Handler();
  }
}

/**
  * @brief GPIO Initialization Function
  * @param None
  * @retval None
  */
static void MX_GPIO_Init(void)
{
  GPIO_InitTypeDef GPIO_InitStruct = {0};

  __HAL_RCC_GPIOC_CLK_ENABLE();
  __HAL_RCC_GPIOH_CLK_ENABLE();
  __HAL_RCC_GPIOA_CLK_ENABLE();
  __HAL_RCC_GPIOB_CLK_ENABLE();

  HAL_GPIO_WritePin(LD2_GPIO_Port, LD2_Pin, GPIO_PIN_RESET);

  GPIO_InitStruct.Pin = B1_Pin;
  GPIO_InitStruct.Mode = GPIO_MODE_IT_FALLING;
  GPIO_InitStruct.Pull = GPIO_NOPULL;
  HAL_GPIO_Init(B1_GPIO_Port, &GPIO_InitStruct);

  GPIO_InitStruct.Pin = LD2_Pin;
  GPIO_InitStruct.Mode = GPIO_MODE_OUTPUT_PP;
  GPIO_InitStruct.Pull = GPIO_NOPULL;
  GPIO_InitStruct.Speed = GPIO_SPEED_FREQ_LOW;
  HAL_GPIO_Init(LD2_GPIO_Port, &GPIO_InitStruct);
}

/**
  * @brief  This function is executed in case of error occurrence.
  * @retval None
  */
void Error_Handler(void)
{
  __disable_irq();

  while (1)
  {
  }
}

#ifdef USE_FULL_ASSERT
void assert_failed(uint8_t *file, uint32_t line)
{
  (void)file;
  (void)line;
}
#endif

---

11. Что должно получиться

После запуска в терминале будет выводиться спектр шумного сигнала. Если построить график по двум столбцам freq_hz и amplitude_x1000, должны быть видны три выраженных пика:

200 Гц — полезный сигнал
600 Гц — первая помеха
900 Гц — вторая помеха

Это показывает, что преобразование Фурье позволяет увидеть частотный состав сигнала и определить, какие именно частоты присутствуют внутри сложного шумного сигнала.

Таким образом, если во временной области сигнал выглядит как непонятная кривая, то после разложения Фурье становится ясно, из каких частот он состоит.