←
Предыдущий
В этом примере разберём, как выполнить разложение сигнала на частотные составляющие с помощью дискретного преобразования Фурье на микроконтроллере STM32F411.
В предыдущем примере мы использовали цифровой фильтр, чтобы убрать помехи из сигнала. Здесь фильтр уже не нужен. Теперь задача другая: взять шумный сигнал и посмотреть, из каких частот он состоит.
Для этого мы программно создадим сигнал, в котором есть три синусоиды:
200 Гц — полезный сигнал
600 Гц — первая помеха
900 Гц — вторая помехаПосле этого выполним разложение Фурье и выведем спектр по UART. На графике должны появиться три основных пика: на 200 Гц, 600 Гц и 900 Гц.
Любой сложный сигнал можно представить как сумму простых синусоид разных частот, амплитуд и фаз. Именно это и показывает преобразование Фурье.
Если посмотреть на сигнал во временной области, мы видим только изменение амплитуды во времени. Например, сигнал может выглядеть как сложная кривая с рябью и помехами.
Но если выполнить преобразование Фурье, мы переходим в частотную область. Там уже видно, какие частоты присутствуют внутри сигнала.
Проще говоря:
Временная область:
смотрим, как сигнал меняется во времени.
Частотная область:
смотрим, из каких частот состоит сигнал.В программе мы формируем шумный сигнал из трёх составляющих:
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После преобразования Фурье мы должны увидеть эти частоты в спектре.
Сигнал внутри 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.
Так как сигнал у нас цифровой, мы используем не непрерывный ряд Фурье, а дискретное преобразование Фурье — 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 ГцDFT работает так: для каждой проверяемой частоты программа умножает входной сигнал на синус и косинус этой частоты.
Если такая частота действительно есть в сигнале, произведения складываются и дают большую сумму.
Если такой частоты в сигнале нет, положительные и отрицательные участки взаимно компенсируются, и сумма получается близкой к нулю.
Поэтому в спектре мы должны увидеть:
200 Гц → большой пик
600 Гц → меньший пик
900 Гц → ещё один пикТак как амплитуды сигналов заранее заданы как 1.00, 0.35 и 0.25, результат должен быть примерно таким:
200 Гц → 1000
600 Гц → 350
900 Гц → 250Значения выводятся в формате x1000, чтобы не использовать печать чисел с плавающей точкой через printf.
Сначала программа формирует массив 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 потом будет разложен по Фурье.
Функция 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);Программа выводит таблицу:
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 ГцЭто и есть частотный состав шумного сигнала.
В этой программе используется прямой расчёт DFT. Он медленнее, чем FFT, но намного проще для понимания.
Для каждого значения частоты программа вручную считает суммы синуса и косинуса. Поэтому хорошо видно, как именно работает разложение сигнала на частоты.
FFT работает быстрее, но внутри алгоритм сложнее. Когда понятен принцип DFT, можно переходить к библиотеке CMSIS-DSP и использовать готовые функции FFT, например:
arm_rfft_fast_f32()Но для учебного проекта лучше начать именно с прямого DFT.
/* 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После запуска в терминале будет выводиться спектр шумного сигнала. Если построить график по двум столбцам freq_hz и amplitude_x1000, должны быть видны три выраженных пика:
200 Гц — полезный сигнал
600 Гц — первая помеха
900 Гц — вторая помехаЭто показывает, что преобразование Фурье позволяет увидеть частотный состав сигнала и определить, какие именно частоты присутствуют внутри сложного шумного сигнала.
Таким образом, если во временной области сигнал выглядит как непонятная кривая, то после разложения Фурье становится ясно, из каких частот он состоит.