CMSIS-DSP学习-函数

## 函数名组成
## #计算/运算 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
/* 计算平均值 */
    arm_mean_f32(
        input,
        BLOCK_SIZE,
        (float32_t *)&meanValue
    );

/* 计算有效值 */
    arm_rms_f32(
        input,
        BLOCK_SIZE,
        (float32_t *)&rmsValue
    );

/* 查找最大值及其下标 */
    arm_max_f32(
        input,
        BLOCK_SIZE,
        (float32_t *)&maxValue,
        (uint32_t *)&maxIndex
    );
    
/* 整体缩放 */
	arm_scale_f32(
        adcFloat,
		SCALE,
        voltage,
        SAMPLE_COUNT
    );

/* 整体偏移(加运算) */
    arm_offset_f32(
        voltage,
        OFFSET,
        voltageAC,    //结果输出数组
        SAMPLE_COUNT
    );

传入参数

传入数组 数组大小 计算结果（指针方式）

FIR 数字滤波器

Finite Impulse Response 有限冲激响应

一般的 FIR 滤波器公式是： y[n] = b0x[n] + b1x[n − 1] + b2x[n − 2] + ⋯ + bM − 1x[n − M + 1]
或者求和形式 $$ y[n] = \sum_{k=0}^{M-1}b_kx[n-k] $$ 可理解为滑动窗口平均

如何使用

使用浮点 FIR，需要定义： 

1
arm_fir_instance_f32 firInstance;
这个结构体会记录：

抽头数量 系数数组地址 状态数组地址

可以把它理解成新建一个FIR 滤波器对象。 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
#include "arm_math.h"

//有几个用于取平均的变量
#define NUM_TAPS   4U

//每次调用 FIR 函数处理多少个新采样点
#define BLOCK_SIZE 16U

//声明对象
static arm_fir_instance_f32 firInstance;  

//FIR 的系数，也就是每个采样点的权重 
static const float32_t firCoeffs[NUM_TAPS] = {
    0.25f, //b3
    0.25f, //b2
    0.25f, //b1
    0.25f  //b0   权重从旧值到新值
};

//保存当前和过去的输入采样值 临时使用
static float32_t firState[
    NUM_TAPS + BLOCK_SIZE - 1U
];

//输入缓存
static float32_t inputBuffer[BLOCK_SIZE];
//输出缓存
static float32_t outputBuffer[BLOCK_SIZE];

//初始化函数
arm_fir_init_f32(
        &firInstance, //FIR实例
        NUM_TAPS,     //几个抽头
        firCoeffs,    //权重数组
        firState,     //临时变量数组
        BLOCK_SIZE    //一次处理多少个数
    );//初始化函数只需要调用一次
    
//处理函数
arm_fir_f32(
    &firInstance,
    inputBuffer,
    outputBuffer,
    BLOCK_SIZE        //本次处理点数 只要比init中的数小就行
);

NUM_TAPS + BLOCK_SIZE - 1U 的确定：

firState 需要容纳：

1
2
过去的输入：NUM_TAPS - 1 个
当前数据块：BLOCK_SIZE 个

所以：

状态数组长度 = (NUM_TAPS − 1) + BLOCK_SIZE

整理后就是：

NUM_TAPS + BLOCK_SIZE − 1

例如：

1
2
NUM_TAPS   = 4
BLOCK_SIZE = 8

状态数组长度：

4 + 8 − 1 = 11

可以概念性地看成：

1
2
3
4
┌────历史数据────┬────────当前数据块────────┐
│  6   7   8     │  9 10 11 12 13 14 15 16 │
└────────────────┴──────────────────────────┘
   4-1 = 3个              8个

总共：

1
3 + 8 = 11个

CMSIS-DSP 在内部会自动管理这些数据，不需要往 firState 中填值。

缺陷

  1. FIR 会引入延迟

    对于长度为 N 的对称线性相位 FIR，其群延迟通常是： $$ D=\frac{N-1}{2}$$ 单位是采样点

    例如 5 抽头对称 FIR：

    $$D=\frac{5-1}{2}=2\text{个采样点}$$

    若采样率为：

    Fs = 6400 Hz

    对应时间延迟：

    $$t_D=\frac{2}{6400} =0.0003125\text{秒} =0.3125\text{ ms}$$

  2. 不一定适合谐波分析

    滑动平均滤波器不仅抑制高频噪声，也会对 50～500 Hz 产生不同程度的幅值衰减

    这会造成：

     
    1
    滤波后的5次谐波幅值 ≠ 原始5次谐波幅值
    

 进而影响 THD。

    IIR 与 Biquad 滤波器

Infinite Impulse Response 无限冲激响应 #### 与FIR的不同 FIR 只使用：

  • 当前输入
  • 过去的输入

IIR 除了使用输入，还会把过去的输出反馈回来。

1
2
FIR：输入历史 → 当前输出
IIR：输入历史 + 输出历史 → 当前输出

所有影响都一直会持续下去，而不是像FIR一样只会在一定范围内作用

一般形式

一个 N 阶 IIR 滤波器的差分方程通常写为：

$$y[n] = \sum_{k=0}^{M} b_kx[n-k] - \sum_{k=1}^{N} a_ky[n-k]$$

展开就是：

$$\begin{aligned} y[n]={}&b_0x[n]+b_1x[n-1]+\cdots+b_Mx[n-M]\\ &-a_1y[n-1]-a_2y[n-2]-\cdots-a_Ny[n-N] \end{aligned}$$

其中：

  • x[n]：当前输入
  • y[n]：当前输出
  • x[n − k]：过去的输入
  • y[n − k]：过去的输出
  • bk：前馈系数
  • ak：反馈系数

IIR 与 FIR 最大的区别，就是公式中包含过去的输出 y[n − k]。 #### Biquad 是什么

Biquad是CMSIS-DSP实现的一个二阶滤波器，可用于级联拼接出高阶IIR

一个二阶滤波器称为：Biquad 二阶节

CMSIS-DSP Direct Form I 使用的公式是：

$$\begin{aligned} y[n] ={}& b_0x[n]+b_1x[n-1]+b_2x[n-2]\\ &+a_1y[n-1]+a_2y[n-2] \end{aligned}$$

其中：

  • b0, b1, b2：前馈系数
  • a1, a2：反馈系数

如何使用

使用前初始化： 

1
arm_biquad_casd_df1_inst_f32 biquadInstance;

定义系数（实际数量为 NUM_STAGES * 5 每个二阶节使用5个参数）： 此处不像 FIR 那样反序排列

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
//1级
static const float32_t biquadCoeffs[5] = {
    b0,
    b1,
    b2,
    a1,
    a2
};

//2级
static const float32_t biquadCoeffs[10] = {
    /* 第1级 */
    b0_1, b1_1, b2_1, a1_1, a2_1,

    /* 第2级 */
    b0_2, b1_2, b2_2, a1_2, a2_2
};

定义状态数组：

Direct Form I 每个 stage 需要四个状态：

x[n − 1], x[n − 2], y[n − 1], y[n − 2]

所以： stateLength = 4 × numStages

一个 stage：

1
static float32_t biquadState[4];

两个 stage：

1
static float32_t biquadState[8];

通用写法：

1
static float32_t biquadState[4U * NUM_STAGES];

初始化 

1
2
3
4
5
6
7
8
9
arm_biquad_cascade_df1_init_f32(
    &biquadInstance, //实例化的biquad滤波对象地址
    NUM_STAGES,      //级联次数（有几个二阶biquad滤波器）
    biquadCoeffs,    //滤波系数
    biquadState      //临时存储数组(状态数组)
);  //同样只需要执行一次 
//无须BLOCK_SIZE
//因为 Direct Form I Biquad 每一级始终只需保存四个历史状态
//状态数组大小不随数据块长度变化。
调用处理函数： 
1
2
3
4
5
6
arm_biquad_cascade_df1_f32(
    &biquadInstance,
    inputBuffer,
    outputBuffer,
    BLOCK_SIZE        //一次处理多少个数
);
对于多级级联： CMSIS-DSP 会自动完成：

1
第1级输出 → 第2级输入

不需要自己创建中间数组。

不同种类的Biquad

arm biquad casd df1 inst f32
命名空间 滤波器 级联cascade(缩写) 类型 实例 数据类型 
1
2
arm_biquad_casd_df1_inst_f32
Arm 提供的、直接Ⅰ型的、32 位浮点的Biquad 级联滤波器实例

类型有：

1. Direct Form I：直接Ⅰ型

CMSIS-DSP 函数：

1
arm_biquad_cascade_df1_f32();

公式是：

y[n] = b0x[n] + b1x[n − 1] + b2x[n − 2] + a1y[n − 1] + a2y[n − 2]

每个二阶节保存四个状态：

1
2
3
4
x[n-1]
x[n-2]
y[n-1]
y[n-2]

所以每个 stage 的状态数组长度为 4：

1
float32_t state[4 * NUM_STAGES];

特点

  • 结构直观，容易理解
  • 支持浮点和定点
  • 支持 f16、f32、q15、q31
  • 状态变量较多

CMSIS-DSP 的 DF1 使用每级 5 个系数和 4 个状态变量。

  1. Direct Form II：直接Ⅱ型

普通 DF2 会把输入延迟线和输出延迟线合并。

理论上每个二阶节只需要两个状态，而不是四个。

不过普通 DF2 的内部状态动态范围可能较大，数值误差和溢出问题通常比 DF1 更明显。

CMSIS-DSP 没有把普通 DF2 作为主要 Biquad 接口，而是提供更常用的转置结构 DF2T：

  1. DF2T：转置直接Ⅱ型

CMSIS-DSP 函数：

1
arm_biquad_cascade_df2T_f32();

它使用下面的计算过程：

y[n] = b0x[n] + d1d1 = b1x[n] + a1y[n] + d2 d2 = b2x[n] + a2y[n]

每一级只保存 d1 和 d2

状态数组长度：

1
float32_t state[2 * NUM_STAGES];

特点

  • 每级只需要两个状态
  • 比 DF1 节省状态内存
  • 浮点计算中常用
  • 内部状态的动态范围可能较大
  • CMSIS-DSP 的 DF2T 只提供浮点实现

CMSIS-DSP 当前提供 f16、f32 和 f64 的 DF2T 接口，每个 stage 使用两个状态变量。

DF1 和 DF2T 对比

项目 DF1 DF2T
每级系数 5 个 5 个
每级状态 4 个 2 个
浮点支持 是 是
Q15/Q31 定点支持 是 否
RAM 占用 较多 较少

FIR 与 IIR 对比

特性 FIR IIR
使用过去输入 是 是
使用过去输出 否 是
是否反馈 否 是
冲激响应 有限 理论上无限
稳定性 一般天然稳定 系数不当会不稳定
相位 容易实现线性相位 通常是非线性相位
相同滤波要求的阶数 通常较高 通常较低
状态内容 输入历史 输入和输出历史
系数敏感性 相对较低 相对较高 
1
2
重视线性相位 → 更倾向FIR
重视低运算量 → 更倾向IIR

FFT分为两种类型

复数FFT 

1
arm_cfft_f32()
实数FFT（常用） 
1
arm_rfft_fast_f32()

参数计算

FFT 最核心有如下参数：

参数 符号 含义
采样率 Fs ADC每秒采多少次
FFT长度 N 一次分析多少个点
频率分辨率 Δf 能分辨多小的频率
频率范围 0~Fs/2 能看到多少频率

选取Fs的时候要注意奈奎斯特定理： $$F_{max}=\frac{Fs}{2}$$ 

1
要采样一个频率，采样频率至少为这个频率的两倍

频率分辨率Δf： $$\Delta f = \frac{Fs}{N}$$ 其中，Fs为采样率 N为采样点数

如何使用

1. 初始化

实例化fft对象 

1
2
arm_rfft_fast_instance_f32 fft
//实例化fft对象
里面保存： - FFT长度 - twiddle表 - 参数

定义输入输出数组 

1
2
float32_t input[FFT_SIZE]; 
float32_t fftOutput[FFT_SIZE];
选定方向： 
1
ifftFlag = 1 或 0;

0：FFT 1：逆FFT

2. 调用初始化函数

定义： 

1
2
3
4
arm_status arm_rfft_fast_init_f32(
	arm_rfft_fast_instance_f32 * S, //fft实例
	uint16_t fftLen //FFT长度
);

3. 执行FFT

1
2
3
4
5
6
arm_rfft_fast_f32(
	&fft,
	input,
	fftOutput,
	0
);

数据输出后，第0个为DC分量

对于复数FFT

1
2
3
4
5
6
void arm_cfft_f16	(	
	const arm_cfft_instance_f16 * 	S,
	float16_t * 	p1, //数据操作在原地进行
	uint8_t 	ifftFlag,
	uint8_t 	bitReverseFlag //是否进行位反转操作
)

bitReverse 是否进行位反转 一般启用

FFT 使用的是 Cooley-Tukey 蝶形算法。计算过程中，数据访问顺序不是普通顺序，而是按照二进制位翻转后的顺序排列。处理后的原始完成数据的顺序是乱的。 

1
2
3
4
5
6
7
8
9
原始       翻转
000  ---> 000  (0)
001  ---> 100  (4)
010  ---> 010  (2)
011  ---> 110  (6)
100  ---> 001  (1)
101  ---> 101  (5)
110  ---> 011  (3)
111  ---> 111  (7)
启用后可恢复正常顺序 ### 其他后处理函数

计算幅值 FFT过后得到的是 N/2 个实部，虚部对，所以接下来需要计算赋值来判断频域分布。 

1
arm_cmplx_mag_f32()
定义： 
1
2
3
4
5
void arm_cmplx_mag_f32(
	float32_t * pSrc, //处理数据源
	float32_t * pDst, //处理后数据输出
	uint32_t numSamples //有多少数据
);

计算相位 

1
arm_cmplx_phase_f32()
可求得角 θ 
1
2
3
4
5
void arm_cmplx_phase_f32(
    const float32_t *pSrc, //输入复数数组
    float32_t *pDst,       //输出相位数组
    uint32_t numSamples
);

找最大值 

1
arm_max_f32()

1
2
3
4
5
6
void arm_max_f32(
    const float32_t * pSrc,
    uint32_t blockSize,
    float32_t * pResult,
    uint32_t * pIndex
);

结果使用传入指针的方式输出