利用DSP进行快速傅里叶变换

流程图

主函数

ADC中断

外部中断10：15线

大致思路

利用ADC采样，并将数据利用DMA搬运至内存
将数据转化为电压，并进行截止频率为100KHz低通滤波
将直流偏置去除后进行快速傅里叶变换，并将变换结果传入串口

初始化

时钟

使能RCC，启用外部晶振
将时钟树安如下配置

中断

中断优先级组设置为0位

按钮

将按钮设置为外部中断触发模式，
将子优先级设置为0

TIM2

TIM2分频设置为84 - 1
设重装载值为10 - 1
触发外部事件

ADC

开启采样通道1
使用TIM2作为外部触发源
启用DMA，方向为外设至内存
使能ADC中断，子优先级为0

串口2

使能串口2，采样异步通讯，波特率9600，位长8，无校验位，1停止位

准备步骤

启用DSP
在编译设置中定义以下内容

USE_HAL_DRIVER,STM32F401xE,ARM_MATH_CM4,ARM_MATH_MATRIX,ARM_MATH_ROUNDING

若不写入以下定义，DSP库内容可能无法调用，具体宏定义需根据Cortex内核型号填写

代码部分

引入以下头文件

main.c

/* USER CODE BEGIN Includes */
#include <stdio.h>
#include <math.h>
#include "arm_math.h"
#include "arm_const_structs.h"
/* USER CODE END Includes */

main.h

/* USER CODE BEGIN Includes */
#include <stdbool.h>//使用标识符
/* USER CODE END Includes */

全局变量

在main.c内定义以下全局变量

/* USER CODE BEGIN PV */
float data[NPT];//转化为电压后的采样数据
uint16_t measured_data[NPT];//ADC采样的原始数据
volatile bool measured_flag = false;//测量标识符，真表示测量完成，假表示未完成
uint8_t measure_times = 0;//采样次数

float data_in[NPT*2]={0};//FFT函数需要的复数数组

float data_out[NPT/2];//存储傅里叶变换结果的数组
float filtered_data[NPT];//存储滤波后电压的数组
int i = 0; 
int j = 0;
#define IIR_NUM_STAGES 2

float iir_state[4 * IIR_NUM_STAGES];

/* 4阶 Butterworth LPF, Fs=100k, Fc=45k */
const float iir_coeffs[5 * IIR_NUM_STAGES] = {
    /* Section 1 */
    0.292893f, 0.585786f, 0.292893f,
    0.171573f, -0.207107f,

    /* Section 2 */
    0.292893f, 0.585786f, 0.292893f,
    0.171573f, -0.207107f
};

arm_biquad_casd_df1_inst_f32 iir_filter;//数字滤波器数据结构体
/* USER CODE END PV */

在完成硬件初始化后，将数字滤波器进行初始化

	arm_biquad_cascade_df1_init_f32(
    &iir_filter,
    IIR_NUM_STAGES,
    (float32_t *)iir_coeffs,
    iir_state
);

重定向

将串口进行重定向

int fputc(int c,FILE*stream)
{
HAL_UART_Transmit(&huart2,(unsigned char *)&c,1,1025);
	return 1;
}

主逻辑

在主循环中写入以下内容：

  /* Infinite loop */
  /* USER CODE BEGIN WHILE */
  while (1)
  {
		 if(measured_flag)//判断是否测量完成
    {
        HAL_TIM_Base_Stop(&htim2);//停止及时去计时
        measured_flag = false;//翻转标识符

        /* 1. ADC 原始数据转电压 */
        for(i = 0; i < NPT; i++)
        {
            data[i] = (float)(measured_data[i]) * 3.3f / 4095.0f;
        }
				float mean = 0.0f;//偏置电压
				for(i = 0; i < NPT; i++)
				{
						mean += data[i];
				}
				mean /= NPT; // 计算平均值（也就是直流分量）

				for(i = 0; i < NPT; i++)
				{
						data[i] -= mean; // 每个数据点减去直流分量
				}
        /* 2.  IIR 低通滤波 */
        arm_biquad_cascade_df1_f32(
            &iir_filter,
            data,
            filtered_data,
            NPT
        );

        /* 3. 构造复数输入（实部=滤波后数据，虚部=0） */
        for(j = 0; j < NPT; j++)
        {
            data_in[j*2]   = filtered_data[j];  // 用滤波后的数据
            data_in[j*2+1] = 0;
        }

        /* 4. FFT */
        arm_cfft_f32(&arm_cfft_sR_f32_len1024, data_in, 0, 1);
        arm_cmplx_mag_f32(data_in, data_out, NPT/2);

        /* 5. 串口输出频谱 */
        for(i = 0; i < 256; i++)
        {
            printf("%f\r\n", data_out[i]);
        }
    }

    HAL_Delay(1000);
    /* USER CODE END WHILE */

    /* USER CODE BEGIN 3 */
  }
  /* USER CODE END 3 */

中断

ADC中断

中断服务程序不写入内容，主要内容位于回调函数中

void HAL_ADC_ConvCpltCallback(ADC_HandleTypeDef *hadc){
		if(hadc == &hadc1){//判断中断是否来自ADC1
			measured_flag = true;//将测量标识符设位真
			HAL_ADC_Stop_DMA(&hadc1);///停止ADC采样
		}
}

外部中断10:15线

中断服务程序不写入内容，主要内容位于回调函数中

void HAL_GPIO_EXTI_Callback(uint16_t GPIO_Pin){
	if(GPIO_Pin == B1_Pin){//判断中断是否来自按钮
		measure_times++;
		HAL_TIM_Base_Start(&htim2);//开始TIM2定时
		HAL_ADC_Start_DMA(&hadc1, (uint32_t*)measured_data, 1024);//开启ADC采样
	}
}

实际演示

输入10k方波
串口显示如下
前256分量如下：

10kHz三角波测量结果如下
前256分量数据如下：

代码讲解

数字滤波

数组部分

/* 4阶 Butterworth LPF, Fs=100k, Fc=45k */
const float iir_coeffs[5 * IIR_NUM_STAGES] = {
    /* Section 1 */
    0.292893f, 0.585786f, 0.292893f,
    0.171573f, -0.207107f,

    /* Section 2 */
    0.292893f, 0.585786f, 0.292893f,
    0.171573f, -0.207107f
};

这组系数描述的是 Direct Form I 结构的差分方程：

\( y[n] = b_0 \cdot x[n] + b_1 \cdot x[n-1] + b_2 \cdot x[n-2] - a_1 \cdot y[n-1] - a_2 \cdot y[n-2] \)
在这个公式中：

\(x[n]\) 是当前输入（你的 ADC 数据）
\(y[n]\) 是当前输出（滤波后的数据）
\(b_0, b_1, b_2\) 是前向系数（Feed-Forward / FIR 部分）
\(a_1, a_2\) 是反馈系数（Feedback / IIR 部分）
二、数组的物理含义（逐行拆解）

数组是这样定义的：

const float iir_coeffs[5 * IIR_NUM_STAGES] = {...};

这里 5 * IIR_NUM_STAGES 是关键。每一个二阶节（Biquad Stage）固定占用 5 个 float 空间。

内存布局规则（必须死记）
CMSIS-DSP 规定，这 5 个数的顺序是：

数组下标	变量名	数学含义
`[0]`	`b0`	当前输入的增益
`[1]`	`b1`	前一次输入的增益
`[2]`	`b2`	前两次输入的增益
`[3]`	`a1`	负的前一次输出的增益
`[4]`	`a2`	负的前两次输出的增益

⚠️ 极度重要的坑点（针对 a1, a2）：
在数学公式里是减号（\(- a_1 \cdot y[n-1]\)）。
但在 CMSIS-DSP 里，你存入 a1 的值必须是 \(-a_1\)。
也就是说，如果你算出来的理论 \(a_1 = -0.171573\)，那你存进数组的应该是 0.171573。
三、对照代码详解

我们把代码拆开来看：

Section 1（第一级滤波）

0.292893f,  // b0: 当前输入权重
0.585786f,  // b1: 上一个输入权重
0.292893f,  // b2: 上上个输入权重
0.171573f,  // a1: -(-0.171573) -> 反馈项1
-0.207107f  // a2: -(0.207107)  -> 反馈项2

观察：b0 和 b2 相等，b1 是它们的两倍左右。这是 Butterworth 滤波器的典型特征（对称性）。
观察：a2 是负数。这符合 DSP 库的要求（因为 \(a_2\) 理论值通常为正，取负存入）。

Section 2（第二级滤波）

0.292893f,  // b0
0.585786f,  // b1
0.292893f,  // b2
0.171573f,  // a1
-0.207107f  // a2

观察：两级参数一模一样。这是因为这是一个 4 阶 (Order 4) Butterworth 滤波器。
高阶滤波器必须拆分成多个 2 阶节（因为高于 2 阶的滤波器无法保证稳定）。对于偶数阶 Butterworth，拆分后的每个 Biquad 参数通常是一样的。

四、这些数字是怎么来的？（设计溯源）

你用的是 Fs=100kHz, Fc=45kHz, 4阶 Butterworth。

如果用 Python (SciPy) 或 MATLAB 设计，逻辑是这样的：

归一化频率：
数字频率 \(\omega_c = 2 \pi \frac{F_c}{F_s} = 2 \pi \frac{45k}{100k} = 0.9\pi\)。
极点计算：
Butterworth 滤波器的极点在 S 平面均匀分布。映射到 Z 平面后，得到两个共轭复极点对。
系数推导：
由极点推导出传递函数 \(H(z)\)，最终整理成 \(\frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}}\) 的形式。

0.292893 实际上是 \(\frac{1}{2 + \sqrt{2}}\) 的近似值，这是这类滤波器在特定截止频率下的固有常数。

六、总结

项目	解释
数组本质	两个串联的二阶滤波器（Biquad）的系数集合
滤波器类型	4 阶 Butterworth 低通
关键顺序	`{b0, b1, b2, a1, a2}`
a1/a2 符号	存入的是理论值的负数（你的代码是对的）
为何两级相同	4 阶 Butterworth 的标准拆分方式

一句话总结：
这串数字是告诉 STM32 的 DSP 单元：“请用这两个特定的‘漏斗’（\(b\) 系数）去混合新旧输入，再用两个特定的‘回音壁’（\(a\) 系数）去混合新旧输出，从而只允许 45kHz 以下的信号通过。”

函数部分

/**
  @brief         浮点双二阶串联滤波器的处理函数
  @param[in]     S         浮点双二阶级联结构体指针
  @param[in]     pSrc      指向输入数据的地址
  @param[out]    pDst      指向输出数据的地址
  @param[in]     blockSize  数据大小
  @return        none
 */
 void arm_biquad_cascade_df1_f32(
  const arm_biquad_casd_df1_inst_f32 * S,
  const float32_t * pSrc,
        float32_t * pDst,
        uint32_t blockSize);

FFT

        /* 4. FFT */
        arm_cfft_f32(&arm_cfft_sR_f32_len1024, data_in, 0, 1);
        arm_cmplx_mag_f32(data_in, data_out, NPT/2);

一、`arm_cfft_f32` —— 核心 FFT 运算

这行代码是真正执行傅里叶变换的地方。

/**
  @brief         浮点复数快速傅里叶变换（FFT）的处理函数
  @param[in]     S              指向FFT结构体的指针
  @param[in,out] p1            指向大小为<code>2*fftLen</code>的复杂数据缓冲区。处理过程在数据执行
  @param[in]     ifftFlag       选择变换方向的标识符
                   - value = 0: 正向变换
                   - value = 1: 逆向变换
  @param[in]     bitReverseFlag 启用/禁用输出位反转的标志
                   - value = 0: 禁用输出位翻转
                   - value = 1: 启用输出位翻转
  @return        none
 */

void arm_cfft_f32(
  const arm_cfft_instance_f32 * S,
        float32_t * p1,
        uint8_t ifftFlag,
        uint8_t bitReverseFlag);

1. 函数功能

它计算复数快速傅里叶变换 (Complex FFT)。
它会读取 data_in 中的复数数据，经过一系列蝴蝶运算，将时域信号转换为频域信号。
⚠️ 关键点：它是原地操作 (In-place)。也就是说，计算完成后，结果会覆盖掉原来的 data_in 数组。

参数详解

参数 1：&arm_cfft_sR_f32_len1024

类型：const arm_cfft_instance_f32*
作用：FFT 的配置结构体（控制块）。
详解：
- len1024：告诉 DSP 库这次要做 1024 点 FFT。
- sR：代表使用 位反转 (Bit-Reversal) 查找表。这是为了提高效率。
- 这个结构体里包含了旋转因子（Twiddle Factors）表和位反转表。你不需要手动修改它，只需根据你的点数选择正确的宏（_len256, _len512, _len1024 等）。

参数 2：data_in

类型：float32_t*
作用：输入/输出缓冲区。

内存布局（极其重要）：
它必须是一个长度为 2 * NPT 的数组。因为 FFT 处理的是复数，需要交替存放实部和虚部：

data_in[0] = 实部(sample0)  // Re[0]
data_in[1] = 虚部(sample0)  // Im[0]
data_in[2] = 实部(sample1)  // Re[1]
data_in[3] = 虚部(sample1)  // Im[1]
...

在你的代码中，你做了这一步：

data_in[j*2]   = filtered_data[j]; // 实部 = 采样值
data_in[j*2+1] = 0;                // 虚部 = 0 (因为ADC是实数信号)

参数 3：ifftFlag

类型：uint8_t
作用：指定正向 FFT 还是反向 FFT (IFFT)。
取值：
- 0：正向 FFT（时域 -> 频域）。你这里用的是这个。
- 1：反向 FFT（频域 -> 时域）。

参数 4：bitReverseFlag

类型：uint8_t
作用：是否执行位反转。
取值：
- 1：执行位反转。这是最常用的设置。FFT 算法为了效率，计算出来的结果顺序是乱的（X(0), X(8), X(4)...），位反转会把这些数据重新排列成正常的频率顺序（X(0), X(1), X(2)...）。
- 0：不执行。除非你对接的底层协议有特殊要求，否则永远填 1。

二、`arm_cmplx_mag_f32` —— 求模（幅值）

FFT 算完之后，data_in 里存的是复数（实部+虚部）。人眼无法直接看出“某个频率的能量有多大”，除非我们计算它的模（Magnitude）。这个函数就是干这个的。

/**
  @brief         浮点复数模值
  @param[in]     pSrc        指向输入向量的指针
  @param[out]    pDst        指向输出向量的指针
  @param[in]     numSamples  向量的元素数量
  @return        none
 */
void arm_cmplx_mag_f32(
  const float32_t * pSrc,
        float32_t * pDst,
        uint32_t numSamples);

1. 函数功能

计算复数数组中每个元素的模（幅值）。
计算公式：\(Mag = \sqrt{Re^2 + Im^2}\)

2. 参数详解

参数 1：data_in

类型：float32_t*
作用：输入源。
注意：它读取的是 刚刚被 arm_cfft_f32 覆盖过的 data_in。此时 data_in 里已经不再是时域波形，而是频域的复数谱。

参数 2：data_out

类型：float32_t*
作用：输出缓冲区。
内容：存放计算出的幅值。
大小：你只需要 NPT/2 个点。为什么？见下文。

参数 3：numSamples

类型：uint32_t
作用：要处理的复数对数量。
为什么是 NPT/2？
这是一个非常经典的 DSP 知识点：
1. 你输入了 NPT 个实数采样点（虚部补0）。
2. 经过 FFT 后，产生了 NPT 个复数输出。
3. 但是，对于实数输入信号（你的 ADC 采样就是实数），其频谱是关于中心点 共轭对称 的。
4. 因此，后半部分的数据是前半部分的镜像，不包含新信息。
5. 我们只关心 0Hz 到 Fs/2 (Nyquist频率) 这段频谱。
6. 这段频谱恰好包含 NPT/2 个有效点（对应 DC 到 50kHz）。
7. 所以你只需要计算并输出前 NPT/2 个点的幅值。

三、数据流全景图（帮助你建立直觉）

让我们把你这两行代码串起来，看看数据是怎么变的：

准备阶段：
filtered_data = [1.0, -1.0, 1.0, -1.0, ...] (1024个实数)
构造复数：
data_in = [1.0, 0, -1.0, 0, 1.0, 0, -1.0, 0, ...] (2048个数)
执行 arm_cfft_f32：
data_in (原地覆盖) = [Mag_DC, Phase_DC, Mag_F1, Phase_F1, ..., Mag_Nyq, Phase_Nyq]
(此时数组里是复数，人类很难直接读)
执行 arm_cmplx_mag_f32：
data_out = [Mag_DC, Mag_F1, Mag_F2, ..., Mag_Nyq]
(此时数组里是纯幅值，就是你 printf 打印出来的那些数据)

四、总结表格

函数	角色	输入	输出	关键参数
`arm_cfft_f32`	引擎	复数时域信号	复数频域信号	`0`(正向), `1`(位反转)
`arm_cmplx_mag_f32`	翻译官	复数频域信号	实数幅值谱	`NPT/2`(只看一半)

posted @ 2026-07-06 19:47 奶龙大王阅读(81) 评论(0) 收藏举报

刷新页面返回顶部

奶龙科技

技术奶龙拯救世界 Tech Nailongs Save the World

利用DSP进行快速傅里叶变换

利用DSP进行快速傅里叶变换

流程图

主函数

ADC中断

外部中断10：15线

大致思路

初始化

时钟

中断

按钮

TIM2

ADC

串口2

准备步骤

代码部分

引入以下头文件

全局变量

重定向

主逻辑

中断

ADC中断

外部中断10:15线

实际演示

代码讲解

数字滤波

数组部分

函数部分

FFT

一、`arm_cfft_f32` —— 核心 FFT 运算

1. 函数功能

参数详解

二、`arm_cmplx_mag_f32` —— 求模（幅值）

1. 函数功能

2. 参数详解

三、数据流全景图（帮助你建立直觉）

四、总结表格

公告

奶龙科技

技术奶龙拯救世界 Tech Nailongs Save the World

利用DSP进行快速傅里叶变换

利用DSP进行快速傅里叶变换

流程图

主函数

ADC中断

外部中断10：15线

大致思路

初始化

时钟

中断

按钮

TIM2

ADC

串口2

准备步骤

代码部分

引入以下头文件

全局变量

重定向

主逻辑

中断

ADC中断

外部中断10:15线

实际演示

代码讲解

数字滤波

数组部分

函数部分

FFT

一、arm_cfft_f32 —— 核心 FFT 运算

1. 函数功能

参数详解

二、arm_cmplx_mag_f32 —— 求模（幅值）

1. 函数功能

2. 参数详解

三、数据流全景图（帮助你建立直觉）

四、总结表格

公告

一、`arm_cfft_f32` —— 核心 FFT 运算

二、`arm_cmplx_mag_f32` —— 求模（幅值）