快速傅里叶变换，IIR数字滤波器使用教程
快速傅里叶变换，IIR数字滤波器使用教程
前置步骤
配置编译环境
- 进入“管理运行环境”
- 保证你的ARMCC版本为五代，并使用微库
![img]()
- 进入
CMSIS，勾选CORE,DSP
![img]()
- 由于后续涉及到的头文件过多，所以需要开启On优化，否则会导致编译结果过大无法下载至烧录器
- 另外需要在宏定义中写入以下内容，否则DSP库中部分函数无法启用:
USE_STDPERIPH_DRIVER,STM32F407xE,ARM_MATH_CM4,ARM_MATH_MATRIX,ARM_MATH_ROUNDING
宏定义
ARM_MATH_CMx中文本取决于你的设备型号，查看方式如下
我的设备为STM32F407VET6，内核为Cortex-M4，因此要再添加：
ARM_MATH_CM4
检查文件
- 请再三检查你的项目里是否包含以下文件，若存在就将他们从项目系统中移除并从自己的项目文件中删除
![img]()
这几个文件已经在CMSIS中包含了，若不删除可能导致重定义
FFT
简介
- FFT是一种DFT的高效算法，称为快速傅里叶变换（fast Fourier transform）。傅里叶变换是时域一频域变换分析中最基本的方法之一。在数字处理领域应用的离散傅里叶变换(DFT：Discrete Fourier Transform)是许多数字信号处理方法的基础
前置步骤
- 引用以下头文件
#include "math.h"
#include "arm_math.h"
#include "arm_const_structs.h"
定义变量
- 需要按照格式输入数据（以时域向频域变换为例）
#define NPT 1024//数据尺寸，要求为2的n次方，这里以2的10次为例
float32_t data[NPT];//数据为双精度浮点，要求为一对复数（偶数项索引值为实部，奇数项索引为虚部）
float32_t data_out[NPT/2];//用于放置频域变换后结果的数组
处理数据
- 在完成数据的初始化后，调用以下函数
arm_cfft_f32(&arm_cfft_sR_f32_len1024,data,0,1);
arm_cmplx_mag_f32(data,data_out,NPT/2);
代码详解
傅里叶变换函数
arm_cfft_f32
/**
@brief 浮点复数快速傅里叶变换（FFT）的处理函数
@param[in] S 指向CFFT结构体的指针
@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);
CFFT结构体
- 该参数取决与原复数函数的尺寸大小，DSP库提供的选择如下，可根据自己数据的尺寸自行选择
定义存于
arm_const_structs.c
/* Floating-point structs */
#if !defined(ARM_MATH_MVEF) || defined(ARM_MATH_AUTOVECTORIZE)
#if !defined(ARM_DSP_CONFIG_TABLES) || defined(ARM_ALL_FFT_TABLES) || (defined(ARM_TABLE_TWIDDLECOEF_F32_16) && defined(ARM_TABLE_BITREVIDX_FLT_16))
const arm_cfft_instance_f32 arm_cfft_sR_f32_len16 = {
16, twiddleCoef_16, armBitRevIndexTable16, ARMBITREVINDEXTABLE_16_TABLE_LENGTH
};
#endif
#if !defined(ARM_DSP_CONFIG_TABLES) || defined(ARM_ALL_FFT_TABLES) || (defined(ARM_TABLE_TWIDDLECOEF_F32_32) && defined(ARM_TABLE_BITREVIDX_FLT_32))
const arm_cfft_instance_f32 arm_cfft_sR_f32_len32 = {
32, twiddleCoef_32, armBitRevIndexTable32, ARMBITREVINDEXTABLE_32_TABLE_LENGTH
};
#endif
#if !defined(ARM_DSP_CONFIG_TABLES) || defined(ARM_ALL_FFT_TABLES) || (defined(ARM_TABLE_TWIDDLECOEF_F32_64) && defined(ARM_TABLE_BITREVIDX_FLT_64))
const arm_cfft_instance_f32 arm_cfft_sR_f32_len64 = {
64, twiddleCoef_64, armBitRevIndexTable64, ARMBITREVINDEXTABLE_64_TABLE_LENGTH
};
#endif
#if !defined(ARM_DSP_CONFIG_TABLES) || defined(ARM_ALL_FFT_TABLES) || (defined(ARM_TABLE_TWIDDLECOEF_F32_128) && defined(ARM_TABLE_BITREVIDX_FLT_128))
const arm_cfft_instance_f32 arm_cfft_sR_f32_len128 = {
128, twiddleCoef_128, armBitRevIndexTable128, ARMBITREVINDEXTABLE_128_TABLE_LENGTH
};
#endif
#if !defined(ARM_DSP_CONFIG_TABLES) || defined(ARM_ALL_FFT_TABLES) || (defined(ARM_TABLE_TWIDDLECOEF_F32_256) && defined(ARM_TABLE_BITREVIDX_FLT_256))
const arm_cfft_instance_f32 arm_cfft_sR_f32_len256 = {
256, twiddleCoef_256, armBitRevIndexTable256, ARMBITREVINDEXTABLE_256_TABLE_LENGTH
};
#endif
#if !defined(ARM_DSP_CONFIG_TABLES) || defined(ARM_ALL_FFT_TABLES) || (defined(ARM_TABLE_TWIDDLECOEF_F32_512) && defined(ARM_TABLE_BITREVIDX_FLT_512))
const arm_cfft_instance_f32 arm_cfft_sR_f32_len512 = {
512, twiddleCoef_512, armBitRevIndexTable512, ARMBITREVINDEXTABLE_512_TABLE_LENGTH
};
#endif
#if !defined(ARM_DSP_CONFIG_TABLES) || defined(ARM_ALL_FFT_TABLES) || (defined(ARM_TABLE_TWIDDLECOEF_F32_1024) && defined(ARM_TABLE_BITREVIDX_FLT_1024))
const arm_cfft_instance_f32 arm_cfft_sR_f32_len1024 = {
1024, twiddleCoef_1024, armBitRevIndexTable1024, ARMBITREVINDEXTABLE_1024_TABLE_LENGTH
};
#endif
#if !defined(ARM_DSP_CONFIG_TABLES) || defined(ARM_ALL_FFT_TABLES) || (defined(ARM_TABLE_TWIDDLECOEF_F32_2048) && defined(ARM_TABLE_BITREVIDX_FLT_2048))
const arm_cfft_instance_f32 arm_cfft_sR_f32_len2048 = {
2048, twiddleCoef_2048, armBitRevIndexTable2048, ARMBITREVINDEXTABLE_2048_TABLE_LENGTH
};
#endif
#if !defined(ARM_DSP_CONFIG_TABLES) || defined(ARM_ALL_FFT_TABLES) || (defined(ARM_TABLE_TWIDDLECOEF_F32_4096) && defined(ARM_TABLE_BITREVIDX_FLT_4096))
const arm_cfft_instance_f32 arm_cfft_sR_f32_len4096 = {
4096, twiddleCoef_4096, armBitRevIndexTable4096, ARMBITREVINDEXTABLE_4096_TABLE_LENGTH
};
复数取模函数
- 由于变换结果依旧为复数，人类难以阅读，所以需要使用取模函数计算绝对值
/**
@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);
IIR
简介
- “递归滤波器”。递归滤波器，也就是IIR数字滤波器，顾名思义，具有反馈。
前置步骤
- 引用以下头文件
#include <math.h>
#include "arm_math.h"
#include "arm_const_structs.h"
- 定义以下数据
// 以Fs=100k, Fc=20k, 4阶 Butterworth为例
#define NPT 1024//数据尺寸为1024
#define IIR_NUM_STAGES 2//2结节
float32_t filtered_data[NPT];//存放滤波后数据的数组
float32_t data[NPT];//存放原始数据的数组
//存放滤波器参数的数组
const float iir_coeffs[5 * IIR_NUM_STAGES] = {
// Section 1
0.136728f, 0.273456f, 0.136728f,
0.451286f, -0.234616f,
// Section 2
0.136728f, 0.273456f, 0.136728f,
0.451286f, -0.234616f
};
arm_biquad_casd_df1_inst_f32 iir_filter;//IIR滤波器结构体
调用函数
arm_biquad_cascade_df1_init_f32
- 初始化IIR滤波器结构体
arm_biquad_cascade_df1_init_f32(
&iir_filter,
IIR_NUM_STAGES,
(float32_t *)iir_coeffs,
iir_state
);
arm_biquad_cascade_df1_f32
- 进行滤波
arm_biquad_cascade_df1_f32(
&iir_filter,
data,
filtered_data,
NPT
);
代码详解
arm_biquad_casd_df1_inst_f32
存于
filter_functions.h
/**
* @brief 浮点双二阶级联滤波器的结构体
*/
typedef struct
{
uint32_t numStages; /**< 二阶节数量 */
float32_t *pState; /**< 指向状态数组. */
const float32_t *pCoeffs; /**< 指向系数数组 */
} arm_biquad_casd_df1_inst_f32;
arm_biquad_cascade_df1_init_f32
/**
@brief 初始化浮点双二阶级联滤波器的结构体
@param[in,out] S 指向浮点双二阶级联滤波器的结构体的指针
@param[in] numStages 浮点双二阶级联滤波器的结构体的阶数
@param[in] pCoeffs 指向系数数组
@param[in] pState 指向状态数组
@return none
@par 系数与状态排序 系数按以下顺序存储在数组<code>pCoeffs</code>中:
<pre>
{b10, b11, b12, a11, a12, b20, b21, b22, a21, a22, ...}
</pre>
@par
其中<code>b1x</code>和<code>a1x</code>是第一阶段的系数，
<code>b2x</code>和<code>a2x</code>是第二阶段的系数，
以此类推。<code>pCoeffs</code>数组总共包含<code>5*numStages</code>个值。
@par
<code>pState</code>是指向状态数组的指针。
每个双二阶环节都有4个状态变量：<code>x[n-1]</code>、<code>x[n-2]</code>、<code>y[n-1]</code>和<code>y[n-2]</code>。
状态变量在<code>pState</code>数组中的排列方式如下
<pre>
{x[n-1], x[n-2], y[n-1], y[n-2]}
</pre>
首先是第一阶段的4个状态变量，然后是第二阶段的4个状态变量，以此类推。
状态数组的总长度为<code>4*numStages</code>个值。
在处理完每个数据块后，状态变量会被更新；系数则保持不变。
@par 对于MVE代码，需要额外增加一个用于存储修改后系数的缓冲区。
其大小为numStages，且该缓冲区的每个元素类型均为arm_biquad_mod_coef_f32。
因此，其总大小为32乘以numStages个float32_t元素。
必须使用的初始化函数是arm_biquad_cascade_df1_mve_init_f32。
*/
void arm_biquad_cascade_df1_init_f32(
arm_biquad_casd_df1_inst_f32 * S,
uint8_t numStages,
const float32_t * pCoeffs,
float32_t * pState);
arm_biquad_cascade_df1_f32
/**
@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);
脚本
- 由于系数的计算过于复杂，所以我写了一个python脚本，可按照需求生成对应系数数据
import numpy as np
from scipy import signal
import sys
def generate_cmsis_iir():
"""交互式生成CMSIS-DSP兼容的IIR滤波器系数"""
print("="*60)
print("STM32 CMSIS-DSP IIR滤波器系数生成工具")
print("支持类型：1.低通(LPF) 2.高通(HPF) 3.带通(BPF) 4.带阻(BSF)")
print("="*60)
# -------------------------- 1. 选择滤波器类型 --------------------------
type_map = {
1: ("lowpass", "低通(LPF)"),
2: ("highpass", "高通(HPF)"),
3: ("bandpass", "带通(BPF)"),
4: ("bandstop", "带阻(BSF)")
}
while True:
try:
ftype_idx = int(input("请选择滤波器类型 (1-4): "))
if ftype_idx in type_map:
ftype_str, ftype_name = type_map[ftype_idx]
break
else:
print("输入无效，请输入1-4之间的数字")
except ValueError:
print("输入无效，请输入数字")
# -------------------------- 2. 输入基础参数 --------------------------
while True:
try:
fs = float(input(f"请输入采样率Fs (Hz，你的项目是100000): "))
if fs <= 0:
print("采样率必须大于0")
continue
break
except ValueError:
print("输入无效，请输入数字")
while True:
try:
order = int(input(f"请输入滤波器阶数 (建议2/4/6/8，默认4): ") or "4")
if order < 2:
print("阶数至少为2")
continue
break
except ValueError:
print("输入无效，请输入整数")
# -------------------------- 3. 输入频率参数 --------------------------
nyq = fs / 2 # 奈奎斯特频率
while True:
try:
if ftype_idx in [1, 2]: # 低通/高通只需要一个截止频率
fc = float(input(f"请输入截止频率Fc (Hz，需小于{nyq:.1f}): "))
if 0 < fc < nyq:
freq_param = fc
break
else:
print(f"截止频率必须在0~{nyq:.1f}Hz之间")
else: # 带通/带阻需要两个频率
flo = float(input(f"请输入低频截止Flo (Hz，需小于{nyq:.1f}): "))
fhi = float(input(f"请输入高频截止Fhi (Hz，需大于Flo且小于{nyq:.1f}): "))
if 0 < flo < fhi < nyq:
freq_param = [flo, fhi]
break
else:
print(f"频率必须满足：0 < Flo < Fhi < {nyq:.1f}Hz")
except ValueError:
print("输入无效，请输入数字")
# -------------------------- 4. 生成SOS系数（核心步骤） --------------------------
print("\n正在生成滤波器系数...")
try:
# 生成级联二阶节(SOS)格式系数，这是CMSIS-DSP最稳定的格式
sos = signal.butter(
order,
freq_param,
btype=ftype_str,
analog=False,
fs=fs,
output='sos'
)
except Exception as e:
print(f"系数生成失败：{e}")
sys.exit(1)
n_sections = sos.shape[0] # 二阶节的数量，等于 ceil(阶数/2)
print(f"生成成功！共{n_sections}个二阶节，请将IIR_NUM_STAGES定义为{n_sections}")
# -------------------------- 5. 转换为CMSIS-DSP格式 --------------------------
# SciPy SOS格式：[b0, b1, b2, a0, a1, a2]，其中a0恒为1
# CMSIS-DSP格式：每个二阶节按 [b0, b1, b2, -a1, -a2] 排列（注意a1/a2取反！）
cmsis_coeffs = []
for i in range(n_sections):
b0, b1, b2 = sos[i, 0], sos[i, 1], sos[i, 2]
a1, a2 = sos[i, 4], sos[i, 5]
cmsis_coeffs.extend([b0, b1, b2, -a1, -a2]) # 关键：a1/a2取反
# -------------------------- 6. 输出C语言数组 --------------------------
print("\n" + "="*60)
print("复制到你的C代码中的内容：")
print("="*60)
# 输出宏定义
print(f"/* {ftype_name} 滤波器参数：Fs={fs}Hz, 阶数={order}")
if ftype_idx in [1,2]:
print(f" 截止频率={freq_param}Hz */")
else:
print(f" 通带/阻带：{freq_param[0]}Hz - {freq_param[1]}Hz */")
print(f"#define IIR_NUM_STAGES {n_sections}")
print(f"#define IIR_COEFF_LEN (5 * IIR_NUM_STAGES)")
print("\n/* CMSIS-DSP IIR系数数组（已自动处理a1/a2符号） */")
print("const float iir_coeffs[IIR_COEFF_LEN] = {")
# 按二阶节分组打印，方便阅读
for i in range(n_sections):
start = i * 5
end = start + 5
coeff_str = ", ".join([f"{x:.6f}f" for x in cmsis_coeffs[start:end]])
if i == n_sections -1:
print(f" {coeff_str} // Section {i+1}")
else:
print(f" {coeff_str}, // Section {i+1}")
print("};")
# 输出状态数组定义提示
print(f"\n/* 别忘了定义状态数组（大小必须为4*IIR_NUM_STAGES） */")
print(f"float iir_state[4 * IIR_NUM_STAGES];")
# 输出初始化代码示例
print("\n/* 初始化滤波器（放在main函数初始化部分） */")
print("arm_biquad_cascade_df1_init_f32(")
print(" &iir_filter,")
print(" IIR_NUM_STAGES,")
print(" (float32_t*)iir_coeffs,")
print(" iir_state")
print(");")
# -------------------------- 7. 绘制幅频响应（可选） --------------------------
try:
import matplotlib.pyplot as plt
plot_flag = input("\n是否绘制幅频响应图？(Y/n，默认Y): ").lower() or "y"
if plot_flag == "y":
# 计算频率响应
w, h = signal.sosfreqz(sos, worN=2000, fs=fs)
# 转换为dB
h_db = 20 * np.log10(np.maximum(abs(h), 1e-10))
plt.figure(figsize=(10, 6))
plt.plot(w, h_db, linewidth=2)
plt.title(f"{ftype_name} 幅频响应 (Fs={fs}Hz, Order={order})")
plt.xlabel("频率 (Hz)")
plt.ylabel("幅值 (dB)")
plt.grid(True, alpha=0.7)
plt.axhline(-3, color='r', linestyle='--', label="-3dB点")
# 标记关键频率
if ftype_idx in [1,2]:
plt.axvline(freq_param, color='g', linestyle='--', label=f"Fc={freq_param}Hz")
else:
plt.axvline(freq_param[0], color='g', linestyle='--', label=f"Flo={freq_param[0]}Hz")
plt.axvline(freq_param[1], color='g', linestyle='--', label=f"Fhi={freq_param[1]}Hz")
plt.axvline(nyq, color='gray', linestyle=':', label=f"奈奎斯特频率({nyq}Hz)")
plt.legend()
plt.tight_layout()
plt.show()
except ImportError:
print("\n未安装matplotlib，跳过绘图。如需绘图请运行：pip install matplotlib")
except Exception as e:
print(f"\n绘图失败：{e}")
if __name__ == "__main__":
generate_cmsis_iir()
浙公网安备 33010602011771号