
充值网站分销站怎么做,广州做网站公司,荷兰网站域名,做网站如何赚钱知乎CMSIS DSP库函数学习笔记 从向量运算到数字滤波的实践指南CMSIS DSP 库函数学习笔记
从向量运算到数字滤波的实践指南
author: hachimi | date: 2026/7/4向量(数组)运算向量运算基本数据类型向量初始化方法向量运算函数…CMSIS DSP库函数学习笔记 从向量运算到数字滤波的实践指南CMSIS DSP 库函数学习笔记
从向量运算到数字滤波的实践指南
author: hachimi | date: 2026/7/4
向量(数组)运算
向量运算基本数据类型向量初始化方法
向量运算函数
加减乘乘 (!?乘乘?!)线性变换和符号操作常用统计函数(aka loss function)
三角函数和开方
正余弦函数反正切函数开方函数和幂函数复数数组和复数小工具
FFT
FFT基本数据类型FFT初始化方法
FFT运算函数
float32_t实数fft与ifftarm_rfft_f32(被弃用)float32_t复数fft与ifft定点数rfft(不好用)定点数cfft(可以用)运行速度测试
数字滤波IIR
概况
代码
浮点数IIR代码定点数IIR代码
FIR
概况
代码
浮点数FIR代码
定点数FIR代码
比较M0上的fast_FIR和原版FIR以及直接用卷积函数测速
矩阵运算
矩阵运算基本数据类型矩阵初始化方法
矩阵运算函数
加减乘法转置与逆向量(数组)运算
向量运算基本数据类型 typedef int8_t q7_t;typedef int16_t q15_t;typedef int32_t q31_t;typedef float float32_t;向量初始化方法
直接q7_t A[10]就行了
向量运算函数
加减乘乘 (!?乘乘?!)
加法、减法、逐元素乘法的接口均一致,pSrcA 操作数A,pSrcB 操作数B,pDst 目标数组C,blockSize 进行运算的长度
但! 注意点积运算的接口,甚至顺序都不一样
pSrcA 操作数A,pSrcB 操作数B,blockSize 进行运算的长度,result 是一个指针,把计算的结果变量的地址传进来
//加法
void arm_add_f32 (float32_t *pSrcA, float32_t *pSrcB, float32_t *pDst, uint32_t blockSize)
void arm_add_q15 (q15_t *pSrcA, q15_t *pSrcB, q15_t *pDst, uint32_t blockSize)
void arm_add_q31 (q31_t *pSrcA, q31_t *pSrcB, q31_t *pDst, uint32_t blockSize)
void arm_add_q7 (q7_t *pSrcA, q7_t *pSrcB, q7_t *pDst, uint32_t blockSize)//减法
void arm_sub_f32 (float32_t *pSrcA, float32_t *pSrcB, float32_t *pDst, uint32_t blockSize)
void arm_sub_q15 (q15_t *pSrcA, q15_t *pSrcB, q15_t *pDst, uint32_t blockSize)
void arm_sub_q31 (q31_t *pSrcA, q31_t *pSrcB, q31_t *pDst, uint32_t blockSize)
void arm_sub_q7 (q7_t *pSrcA, q7_t *pSrcB, q7_t *pDst, uint32_t blockSize)//元素乘法
void arm_mult_f32 (float32_t *pSrcA, float32_t *pSrcB, float32_t *pDst, uint32_t blockSize)
void arm_mult_q15 (q15_t *pSrcA, q15_t *pSrcB, q15_t *pDst, uint32_t blockSize)
void arm_mult_q31 (q31_t *pSrcA, q31_t *pSrcB, q31_t *pDst, uint32_t blockSize)
void arm_mult_q7 (q7_t *pSrcA, q7_t *pSrcB, q7_t *pDst, uint32_t blockSize)//点积
void arm_dot_prod_f32 (float32_t *pSrcA, float32_t *pSrcB, uint32_t blockSize, float32_t *result)
void arm_dot_prod_q15 (q15_t *pSrcA, q15_t *pSrcB, uint32_t blockSize, q63_t *result)
void arm_dot_prod_q31 (q31_t *pSrcA, q31_t *pSrcB, uint32_t blockSize, q63_t *result)
void arm_dot_prod_q7 (q7_t *pSrcA, q7_t *pSrcB, uint32_t blockSize, q31_t *result)示例代码,还是得注意定点数的乘法运算
q7_t A[10] = {11,21,31,41,51,61,71,81,91,101};
q7_t B[10] = {7,8,9,10,11,12,13,14,15,16};
q7_t C[10] = {0};arm_add_q7(A, B, C, 10);
printf("A + B = C\n");
for(int i = 0; i 10; i++) {printf("%4d ", C[i]);
}
printf("\n");arm_sub_q7(A, B, C, 10);
printf("A - B = C\n");
for(int i = 0; i 10; i++) {printf("%4d ", C[i]);
}
printf("\n");arm_mult_q7(A, B, C, 10);
printf("A ⊙ B = C\n");
for(int i = 0; i 10; i++) {printf("%4d ", C[i]);
}
printf("\n");q31_t result = 0;
arm_dot_prod_q7(A, B,10,result);
printf("A · B = %d\n", result);运行结果
A + B = C18 29 40 51 62 73 84 95 106 117
A - B = C4 13 22 31 40 49 58 67 76 85
A ⊙ B = C0 1 2 3 4 5 7 8 10 12
A · B = 7265线性变换和符号操作
直接缩放的函数有arm_scale_q15,arm_shift_q15,实现对向量进行缩放,即对向量vvv进行v′=Cvv' = Cvv′=Cv,CCC为常数
arm_negate_q15和arm_abs_q15分别进行相反数和绝对值的操作
我将常数乘向量 向量整体左移右移、绝对值、相反数都归为缩放函数
接口的规律是,只有定点数才有左移右移的功能,因此在scale函数里面会多一个shift参数,正数实现左移,负数实现右移
最后一个就是给数组加一个直流量的函数,arm_offset_q15
这个线性变换的函数就不想试了,懒得给示例代码了
后来用了两下,发现这五种线性变换函数都可以把源数组和目标数组设置成同一个数组,节省空间
//相乘
void arm_scale_f32 (float32_t *pSrc, float32_t scale, float32_t *pDst, uint32_t blockSize)
void arm_scale_q15 (q15_t *pSrc, q15_t scaleFract, int8_t shift, q15_t *pDst, uint32_t blockSize)
void arm_scale_q31 (q31_t *pSrc, q31_t scaleFract, int8_t shift, q31_t *pDst, uint32_t blockSize)
void arm_scale_q7 (q7_t *pSrc, q7_t scaleFract, int8_t shift, q7_t *pDst, uint32_t blockSize)
//定点数限定,正数左移,负数右移
void arm_shift_q15 (q15_t *pSrc, int8_t shiftBits, q15_t *pDst, uint32_t blockSize)
void arm_shift_q31 (q31_t *pSrc, int8_t shiftBits, q31_t *pDst, uint32_t blockSize)
void arm_shift_q7 (q7_t *pSrc, int8_t shiftBits, q7_t *pDst, uint32_t blockSize)
//相反数
void arm_negate_f32 (float32_t *pSrc, float32_t *pDst, uint32_t blockSize)
void arm_negate_q15 (q15_t *pSrc, q15_t *pDst, uint32_t blockSize)
void arm_negate_q31 (q31_t *pSrc, q31_t *pDst, uint32_t blockSize)
void arm_negate_q7 (q7_t *pSrc, q7_t *pDst, uint32_t blockSize)
//绝对值
void arm_abs_f32 (float32_t *pSrc, float32_t *pDst, uint32_t blockSize)
void arm_abs_q15 (q15_t *pSrc, q15_t *pDst, uint32_t blockSize)
void arm_abs_q31 (q31_t *pSrc, q31_t *pDst, uint32_t blockSize)
void arm_abs_q7 (q7_t *pSrc, q7_t *pDst, uint32_t blockSize)
//偏置
void arm_offset_f32 (float32_t *pSrc, float32_t offset, float32_t *pDst, uint32_t blockSize)
void arm_offset_q15 (q15_t *pSrc, q15_t offset, q15_t *pDst, uint32_t blockSize)
void arm_offset_q31 (q31_t *pSrc, q31_t offset, q31_t *pDst, uint32_t blockSize)
void arm_offset_q7 (q7_t *pSrc, q7_t offset, q7_t *pDst, uint32_t blockSize)常用统计函数(aka loss function)
CMSIS DSP库提供极值查找,均值、有效值、功率计算,甚至还有方差标准差计算(太冷门了,懒得写)
这些函数支持f32 q31 q15 q7的数据 类型,接口都是一致的
pSrc 待统计的数组
blockSize 数组大小
pResult 存放结果的变量的地址(一定是地址)
pIndex 结果在数组中的下标 的地址(一定要用一个指针存)
C语言没办法返回俩变量,所以用了指针这种方式输出变量
这几个函数还是挺有用的,免得总是自己for循环自己写
//最大值
void arm_max_f32 (float32_t *pSrc, uint32_t blockSize, float32_t *pResult, uint32_t *pIndex)
void arm_max_q15 (q15_t *pSrc, uint32_t blockSize, q15_t *pResult, uint32_t *pIndex)
void arm_max_q31 (q31_t *pSrc, uint32_t blockSize, q31_t *pResult, uint32_t *pIndex)
void arm_max_q7 (q7_t *pSrc, uint32_t blockSize, q7_t *pResult, uint32_t *pIndex)
//最小值
void arm_min_f32 (float32_t *pSrc, uint32_t blockSize, float32_t *pResult, uint32_t *pIndex)
void arm_min_q15 (q15_t *pSrc, uint32_t blockSize, q15_t *pResult, uint32_t *pIndex)
void arm_min_q31 (q31_t *pSrc, uint32_t blockSize, q31_t *pResult, uint32_t *pIndex)
void arm_min_q7 (q7_t *pSrc, uint32_t blockSize, q7_t *pResult, uint32_t *pIndex)
//平均值
void arm_mean_f32 (float32_t *pSrc, uint32_t blockSize, float32_t *pResult)
void arm_mean_q15 (q15_t *pSrc, uint32_t blockSize, q15_t *pResult)
void arm_mean_q31 (q31_t *pSrc, uint32_t blockSize, q31_t *pResult)
void arm_mean_q7 (q7_t *pSrc, uint32_t blockSize, q7_t *pResult)
//功率
void arm_power_f32 (float32_t *pSrc, uint32_t blockSize, float32_t *pResult)
void arm_power_q15 (q15_t *pSrc, uint32_t blockSize, q63_t *pResult)
void arm_power_q31 (q31_t *pSrc, uint32_t blockSize, q63_t *pResult)
void arm_power_q7 (q7_t *pSrc, uint32_t blockSize, q31_t *pResult)
//有效值
void arm_rms_f32 (float32_t *pSrc, uint32_t blockSize, float32_t *pResult)
void arm_rms_q15 (q15_t *pSrc, uint32_t blockSize, q15_t *pResult)
void arm_rms_q31 (q31_t *pSrc, uint32_t blockSize, q31_t *pResult)有点无聊,让AI写了个测试代码
float32_t data[10] = {1.2f, -3.4f, 5.6f, 7.8f, -9.0f,2.3f, 4.5f, -6.7f, 8.9f, 0.1f
};
float32_t result_f32 = 0.0f;
uint32_t index = 0;arm_max_f32(data, 10, result_f32, index);
printf("最大值:%.3f,位置:%u\n", result_f32, (unsigned int)index);arm_min_f32(data, 10, result_f32, index);
printf("最小值:%.3f,位置:%u\n", result_f32, (unsigned int)index);arm_mean_f32(data, 10, result_f32);
printf("平均值:%.3f\n", result_f32);arm_power_f32(data, 10, result_f32);
printf("功率:%.3f\n", result_f32);arm_rms_f32(data, 10, result_f32);
printf("有效值:%.3f\n", result_f32);输出结果
最大值:8.900,位置:8
最小值:-9.000,位置:4
平均值:1.130
功率:335.850
有效值:5.795三角函数和开方
正弦余弦和 开方概念挺简单的,但是考虑到定点数的抽象特性,这里还是单开一小节专门说一下,练一下
CMSISdsp库里面的正余弦都挺有特色的,会内置一个粗略的查找表,对你输入的相位值查找最近的两个表内的相位,然后插值
正余弦函数
float32_t arm_cos_f32 (float32_t x)
q15_t arm_cos_q15 (q15_t x)
q31_t arm_cos_q31 (q31_t x)float32_t arm_sin_f32 (float32_t x)
q15_t arm_sin_q15 (q15_t x)
q31_t arm_sin_q31 (q31_t x)浮点数没啥好说的,弧度制,输入[0,2π][0,2\pi][0,2π],不在范围内的自动折叠到这个范围
主要是怎么用定点数,q31_t定点数的三角函数的输入[0x80000000,0x7fffffff][0x80000000,0x7fffffff][0x80000000,0x7fffffff]对应真实的[−2π,+2π][-2\pi,+2\pi][−2π,+2π],输出的范围也是[0x80000000,0x7fffffff][0x80000000,0x7fffffff][0x80000000,0x7fffffff],对应真实的[−1,+1][-1,+1][−1,+1]。对比一下生成某一频率的浮点正弦数组和定点数组的方法
设置采样频率是48000Hz,信号频率是2400,那么在数字域看来,20点出现一次完整的正弦周期
#define sig_len 128
q31_t q_sin[sig_len];
float32_t f_sin[sig_len];
q31_t fs = 48000,f0 = 2400;
//生成浮点正弦数组
for(int i = 0; i sig_len; i++)f_sin[i] = arm_sin_f32(2 * PI * f0 * i / fs);
//生成定点正弦数组
q31_t phase = 0,step = (q31_t)(f0 * (q63_t)0x80000000 / fs);
for(int i = 0; i sig_len; i++){q_sin[i] = arm_sin_q31(phase);phase += step;
}
for(int i = 0; i sig_len; i++)printf("%f,%d\n",f_sin[i],q_sin[i]); 输出结果区别就是总长度的常量,一个是2*pi,一个是0x80000000。如果把q31_t换成q15_t,总长度要写成0x8000
需要注意的是,定点的正弦输出很大,有时候为了完成相加的任务,需要适量右移,不然很容易溢出。
反正切函数
CMSIS DSP库里面我只找到了反正切的函数,函数接口如下
很清楚,一个y,一个x,算出一个结果,放在result里面,注意数据的格式float32_t的函数返回的数据是弧度制,定点数的也是,但是定点数的格式不再是q31_t和q15_t,毕竟范围不包括[−2π,2π][-2\pi,2\pi][−2π,2π]。
对于q31_t的函数输出的值,是q2_29,q15_t返回的是q2_13。
arm_status arm_atan2_f32(float32_t y,float32_t x,float32_t *result)
arm_status arm_atan2_q31(q31_t y,q31_t x,q31_t *result)
arm_status arm_atan2_q15(q15_t y,q15_t x,q15_t *result)测试代码
float32_t result;
arm_atan2_f32(1.0f, 1.0f, result);
printf("atan2(%f, %f) = %f\n", 1.0f, 1.0f, result);q31_t result_q = 0;
arm_atan2_q31(0x40000000, 0x20000000, result_q);
printf("atan2_q31(%d, %d) = %d\n", 0x40000000, 0x20000000, result_q);输出结果
atan2(1.000000, 1.000000) = 0.785398
atan2_q31(1073741824, 536870912) = 594395943计算一下594395943229=1.107....=arctan(2)\frac{594395943}2^{29} = 1.107.... = \arctan(2)259439594329=1.107....=arctan(2)
挺好
开方函数和幂函数
arm_status arm_sqrt_q15 (q15_t in, q15_t *pOut)
arm_status arm_sqrt_q31 (q31_t in, q31_t *pOut)
arm_status arm_sqrt_f32 (float32_t in, float32_t *pOut)void arm_power_f32(const float32_t * pSrc,uint32_t blockSize,float32_t * pResult)
void arm_power_q31(const q31_t * pSrc,uint32_t blockSize,q63_t * pResult)这里需要注意定点数的开方,由于定点数的范围是[-1,1],所以会出现定点数开方越开越大的现象
float32_t result;
arm_sqrt_f32(4.0f, result);
printf("sqrt(%f) = %f\n", 4.0f, result);
q31_t result_q = 0;
arm_sqrt_q31(114514, result_q);
printf("sqrt_q31(%d) = %d\n", 114514, result_q);输出结果
sqrt(4.000000) = 2.000000
sqrt_q31(114514) = 15681739计算114514231=15681739231\sqrt{\frac{114514}{2^{31}}} = \frac{15681739}{2^{31}}231114514=23115681739
幂函数也是一个道理,不多做解释了(不过函数的接口换成了数组)
复数数组和复数小工具
复数主要是用在FFT里面,FFT的输出是复数,因此需要借助复数计算工具。有时候,一些电路参数用复数计算也会很方便,因此这里加一节关于复数计算的章节
复数计算有求模、共轭、数组乘法(包括乘实数、虚数)、点积result=Σn=0N−1a[n]×b[n]result = \Sigma_{n = 0}^{N-1}a[n]\times b[n]result=Σn=0N−1a[n]×b[n]
/* 复数共轭 */
void arm_cmplx_conj_f32(const float32_t *pSrc, float32_t *pDst, uint32_t numSamples);
void arm_cmplx_conj_q31(const q31_t *pSrc, q31_t *pDst, uint32_t numSamples);
void arm_cmplx_conj_q15(const q15_t *pSrc, q15_t *pDst, uint32_t numSamples);/* 复数求模(幅度) */
void arm_cmplx_mag_f32(const float32_t *pSrc, float32_t *pDst, uint32_t numSamples);
void arm_cmplx_mag_q31(const q31_t *pSrc, q31_t *pDst, uint32_t numSamples);
void arm_cmplx_mag_q15(const q15_t *pSrc, q15_t *pDst, uint32_t numSamples);/* 复数乘复数 */
void arm_cmplx_mult_cmplx_f32(const float32_t *pSrcA, const float32_t *pSrcB,float32_t *pDst, uint32_t numSamples);
void arm_cmplx_mult_cmplx_q31(const q31_t *pSrcA, const q31_t *pSrcB,q31_t *pDst, uint32_t numSamples);
void arm_cmplx_mult_cmplx_q15(const q15_t *pSrcA, const q15_t *pSrcB,q15_t *pDst, uint32_t numSamples);/* 复数乘实数 */
void arm_cmplx_mult_real_f32(const float32_t *pSrcCmplx, const float32_t *pSrcReal,float32_t *pDst, uint32_t numSamples);
void arm_cmplx_mult_real_q31(const q31_t *pSrcCmplx, const q31_t *pSrcReal,q31_t *pDst, uint32_t numSamples);
void arm_cmplx_mult_real_q15(const q15_t *pSrcCmplx, const q15_t *pSrcReal,q15_t *pDst, uint32_t numSamples);/* 复数点积 */
void arm_cmplx_dot_prod_f32(const float32_t *pSrcA, const float32_t *pSrcB,uint32_t numSamples, float32_t *realResult, float32_t *imagResult);
void arm_cmplx_dot_prod_q31(const q31_t *pSrcA, const q31_t *pSrcB,uint32_t numSamples, q63_t *realResult, q63_t *imagResult);
void arm_cmplx_dot_prod_q15(const q15_t *pSrcA, const q15_t *pSrcB,uint32_t numSamples, q31_t *realResult, q31_t *imagResult);代码如下,注意numSamples是复数的数量。
float32_t A[sig_len * 2] = {},B[sig_len * 2] = {},result_mul[sig_len * 2] = {},mag[sig_len] = {};
float32_t fs = 48000;
float32_t f = 1000;
for(int i = 0;i sig_len * 2;i++){A[i*2] = arm_cos_f32(2 * PI * f / fs * i);A[i*2 + 1] = arm_sin_f32(2 * PI * f / fs * i);
}printf("A:\n");
for(int i = 0;i sig_len;i++){printf("%f + j%f\n", A[i*2], A[i*2 + 1]);
}arm_cmplx_conj_f32(A, B, sig_len);
printf("A's conjugate:\n");
for(int i = 0;i sig_len;i++){printf("%f + j%f\n", B[i*2], B[i*2 + 1]);
}
arm_cmplx_mag_f32(B, mag, sig_len);
printf("Magnitude of A's conjugate:\n");
for(int i = 0;i sig_len;i++){printf("%f\n", mag[i]);
}
arm_cmplx_mult_cmplx_f32(A, B, result_mul, sig_len);
printf("A * B:\n");
for(int i = 0;i sig_len;i++){printf("%f + j%f\n", result_mul[i*2], result_mul[i*2 + 1]);
}printf("A * B:\n");
arm_cmplx_mult_real_f32(A, mag, result_mul, sig_len);
printf("A * |B| 复数×实数:\n");
for(int i = 0;i sig_len;i++){printf("%f + j%f\n", result_mul[i*2], result_mul[i*2 + 1]);
}float32_t result_real = 0,result_imag = 0;
printf("A dot B:\n");
arm_cmplx_dot_prod_f32(A, B, sig_len, result_real,result_imag);
printf("Real: %f, Imaginary: %f\n", result_real, result_imag);输出结果
A:
1.000000 + j0.000000
0.991428 + j0.130524
0.965910 + j0.258815
0.923880 + j0.382683
0.866011 + j0.499992
0.793340 + j0.608751
0.707107 + j0.707107
0.608751 + j0.793340
A's conjugate:
1.000000 + j0.000000
0.991428 + j-0.130524
0.965910 + j-0.258815
0.923880 + j-0.382683
0.866011 + j-0.499992
0.793340 + j-0.608751
0.707107 + j-0.707107
0.608751 + j-0.793340
Magnitude of A's conjugate:
1.000000
0.999983
0.999983
1.000000
0.999983
0.999983
1.000000
0.999983
A * B:
1.000000 + j0.000000
0.999967 + j0.000000
0.999967 + j0.000000
1.000000 + j0.000000
0.999966 + j0.000000
0.999967 + j0.000000
1.000000 + j0.000000
0.999967 + j0.000000
A * B:
A * |B| 复数×实数:
1.000000 + j0.000000
0.991412 + j0.130522
0.965894 + j0.258810
0.923879 + j0.382683
0.865996 + j0.499983
0.793327 + j0.608741
0.707107 + j0.707107
0.608741 + j0.793327
A dot B:
Real: 7.999833, Imaginary: 0.000000FFT
FFT基本数据类型
//复数fft实例对象
arm_cfft_instance_q15;
arm_cfft_instance_q31;
arm_cfft_instance_f32;
arm_cfft_instance_f64;//实数fft实例对象(定点数的由于一些原因用不了)
arm_rfft_fast_instance_f32;
arm_rfft_fast_instance_f64;
arm_rfft_instance_f32;
arm_rfft_instance_f64;FFT初始化方法
用前面的6个函数就够了,其他的函数固定住了点数,防止初始化太多的时候,会开启很多很多的查找表把flash占满
也可以开编译优化,就不会有这个问题
arm_status arm_cfft_init_q15(arm_cfft_instance_q15 *S, uint32_t fftLen);
arm_status arm_cfft_init_q31(arm_cfft_instance_q31 *S, uint32_t fftLen);
arm_status arm_cfft_init_f32(arm_cfft_instance_f32 *S, uint32_t fftLen);
arm_status arm_cfft_init_f64(arm_cfft_instance_f64 *S, uint32_t fftLen);void arm_rfft_fast_f32(const arm_rfft_fast_instance_f32 *S,float32_t *p,float32_t *pOut,uint8_t ifftFlag);
void arm_rfft_fast_f64(const arm_rfft_fast_instance_f64 *S,float64_t *p,float64_t *pOut,uint8_t ifftFlag);arm_status arm_cfft_init_4096_q15(arm_cfft_instance_q15 * S);
arm_status arm_cfft_init_2048_q15(arm_cfft_instance_q15 * S);
arm_status arm_cfft_init_1024_q15(arm_cfft_instance_q15 * S);
arm_status arm_cfft_init_512_q15(arm_cfft_instance_q15 * S);
arm_status arm_cfft_init_256_q15(arm_cfft_instance_q15 * S);
arm_status arm_cfft_init_128_q15(arm_cfft_instance_q15 * S);
arm_status arm_cfft_init_64_q15(arm_cfft_instance_q15 * S);
arm_status arm_cfft_init_32_q15(arm_cfft_instance_q15 * S);
arm_status arm_cfft_init_16_q15(arm_cfft_instance_q15 * S);arm_status arm_cfft_init_4096_q31(arm_cfft_instance_q31 * S);
arm_status arm_cfft_init_2048_q31(arm_cfft_instance_q31 * S);
arm_status arm_cfft_init_1024_q31(arm_cfft_instance_q31 * S);
arm_status arm_cfft_init_512_q31(arm_cfft_instance_q31 * S);
arm_status arm_cfft_init_256_q31(arm_cfft_instance_q31 * S);
arm_status arm_cfft_init_128_q31(arm_cfft_instance_q31 * S);
arm_status arm_cfft_init_64_q31(arm_cfft_instance_q31 * S);
arm_status arm_cfft_init_32_q31(arm_cfft_instance_q31 * S);
arm_status arm_cfft_init_16_q31(arm_cfft_instance_q31 * S);arm_status arm_cfft_init_4096_f32(arm_cfft_instance_f32 * S);
arm_status arm_cfft_init_2048_f32(arm_cfft_instance_f32 * S);
arm_status arm_cfft_init_1024_f32(arm_cfft_instance_f32 * S);
arm_status arm_cfft_init_512_f32(arm_cfft_instance_f32 * S);
arm_status arm_cfft_init_256_f32(arm_cfft_instance_f32 * S);
arm_status arm_cfft_init_128_f32(arm_cfft_instance_f32 * S);
arm_status arm_cfft_init_64_f32(arm_cfft_instance_f32 * S);
arm_status arm_cfft_init_32_f32(arm_cfft_instance_f32 * S);
arm_status arm_cfft_init_16_f32(arm_cfft_instance_f32 * S);arm_status arm_cfft_init_4096_f64(arm_cfft_instance_f64 * S);
arm_status arm_cfft_init_2048_f64(arm_cfft_instance_f64 * S);
arm_status arm_cfft_init_1024_f64(arm_cfft_instance_f64 * S);
arm_status arm_cfft_init_512_f64(arm_cfft_instance_f64 * S);
arm_status arm_cfft_init_256_f64(arm_cfft_instance_f64 * S);
arm_status arm_cfft_init_128_f64(arm_cfft_instance_f64 * S);
arm_status arm_cfft_init_64_f64(arm_cfft_instance_f64 * S);
arm_status arm_cfft_init_32_f64(arm_cfft_instance_f64 * S);
arm_status arm_cfft_init_16_f64(arm_cfft_instance_f64 * S);FFT运算函数
float32_t实数fft与ifft
void arm_rfft_fast_f32(arm_rfft_fast_instance_f32 *S,float32_t *p,float32_t *pOut,uint8_t ifftFlag) S:rFFT实例化,p: 原数组,pOut: 目标数组,ifftFlag:为0时,表示FFT,为1时表示iFFT
使用感觉:
支持32到4096点,但是G3507极限也就1024点,栈开到0x2000
输入64点实数序列,计算得到64点的复数单边谱,实部和虚部交替排布的形式。然后再转成幅度谱相位谱,就只剩32点了。fft完了之后,需要对数组除以总点数,这样频域上的信号幅度就能够对上时域信号的幅度(当然如果你不在意幅度,也可以不缩放。如果用的是定点数,直接右移更快)。由于只有单边谱,幅度谱的大小会小0.5倍,基波不变。
逆变换很方便,例如,如果你已经使用了
arm_rfft_fast_f32(S, sig, fft_out, 0);
那么你只用
arm_rfft_fast_f32(S,fft_out,sig, 1)
就可以把频域信息变成时域信息
示例代码
float32_t sig[sig_len] = {0};
float32_t fft_out[sig_len] = {0};
float32_t fft_mag[sig_len1] = {0};
float32_t fft_phase[sig_len1] = {0};
//合成一个测试的数组,将在频域上有两个频率分量,分别是17和3,以及基波,预计幅度分别为0.5,0.5,3
//实数fft是复数偶对称谱,dsp库只输出一半长度的复数频谱,实部虚部间隔
for (int i = 0; i sig_len; i++) {sig[i] = arm_cos_f32(2 * PI * i / sig_len * 17) + arm_sin_f32(2 * PI * i / sig_len * 3) + 3;
}
arm_rfft_fast_instance_f32 S;
arm_rfft_fast_init_f32(S, sig_len);//初始化fast rfft实例
arm_rfft_fast_f32(S, sig, fft_out, 0);//计算fft
arm_scale_f32(fft_out, 1.0 / sig_len, fft_out, sig_len);//64点实数fft实际需要除以64,得到的数组是半边的复数频谱,为了得到幅度谱还需要求模,幅度谱的长度是32点
arm_cmplx_mag_f32(fft_out,fft_mag,sig_len1);//幅度谱
for(int i = 0; i (sig_len1); i++) {//相位谱arm_atan2_f32(fft_out[i*2+1], fft_out[i*2], fft_phase[i]);
}for(int i = 0; i (sig_len1); i++) {printf("%f,%f\n", fft_mag[i], fft_phase[i]);
}运行结果
3.000000,0.000000 - 基波
0.000000,0.549134
0.000000,2.663975
0.500000,-1.570796 - 3
0.000000,2.454789
0.000000,-0.037095
0.000000,0.155140
0.000000,-2.711984
0.000000,1.254453
0.000000,-1.898380
0.000000,-0.093135
0.000000,1.703919
0.000000,-1.740963
0.000000,2.905744
0.000000,0.198857
0.000000,2.097739
0.000000,2.761086
0.500000,-0.000001 - 17
0.000000,-0.592633
0.000001,1.091771
0.000000,-2.990585
0.000000,-0.841667
0.000000,0.698781
0.000000,0.324189
0.000000,1.143796
0.000000,-3.099275
0.000000,2.357303
0.000000,-0.005645
0.000000,-2.013420
0.000000,2.677945
0.000000,1.004309
0.000000,-2.248130
VOFA收到的幅度谱和相位谱arm_rfft_f32(被弃用)
就不写了,这一系列的就不写了float32_t复数fft与ifft
复数fft接口
arm_cfft_f32(const arm_cfft_instance_f32 * S,float32_t * pSrc,uint8_t ifftFlag,uint8_t bitReverseFlag)S :复数fft实例指针,用arm_cfft_init_f32初始化
pSrc: 需要fft的数组,执行之后,频域的数据会覆盖这个数组,就是说这个数组充当数据输出与输入的角色
ifftFlag:0表示正向fft,1表示iFFT
bitReverseFlag:尽量选1,填0之后数据我没看懂,所以不做解释了
实例代码如下
float32_t sig[sig_len*2] = {0};
float32_t fft_mag[sig_len] = {0};
float32_t fft_phase[sig_len] = {0};
//合成复数数组,将在频域上有两个频率分量,分别是17和3,以及基波,预计幅度分别为0.5,0.5,3
//实数fft是复数偶对称谱,dsp库只输出一半长度的复数频谱,实部虚部间隔
for (int i = 0; i sig_len; i++) {sig[i*2] = arm_cos_f32(2 * PI * i / sig_len * 17) + 3*arm_sin_f32(2 * PI * i / sig_len * 3) + 3; // 实部sig[i*2+1] = 0; // 虚部
}
arm_cfft_instance_f32 S;
arm_cfft_init_f32(S, sig_len);
arm_cfft_f32(S, sig, 0, 1);//64点arm_scale_f32(sig, 1.0 / sig_len, sig, sig_len*2);//64点实数fft实际需要除以64,得到的数组是半边的复数频谱,为了得到幅度谱还需要求模,幅度谱的长度是32点
arm_cmplx_mag_f32(sig,fft_mag,sig_len);//幅度谱
for(int i = 0; i (sig_len); i++) {//相位谱arm_atan2_f32(sig[i*2+1], sig[i*2], fft_phase[i]);
}for(int i = 0; i (sig_len); i++) {printf("%f,%f\n", fft_mag[i], fft_phase[i]);
}输出结果
3.000000,0.000000
0.000000,0.551655
0.000000,2.332161
1.500000,-1.570796
0.000000,2.662353
0.000000,0.550166
0.000000,1.407460
0.000000,-2.155067
0.000000,0.819699
0.000000,-1.701072
0.000000,-0.252930
0.000000,1.625230
0.000000,-2.470974
0.000000,0.000000
0.000000,0.467486
0.000000,2.254616
0.000000,3.141593
0.500000,-0.000001
0.000000,-0.563053
0.000001,1.210391
0.000000,-3.035063
0.000000,-0.861556
0.000000,0.549403
0.000000,0.380506
0.000000,1.494440
0.000000,-2.957358
0.000000,1.663686
0.000000,-0.339818
0.000000,-2.066445
0.000000,0.000000
0.000000,1.370536
0.000000,-2.259967
0.000000,0.000000
0.000000,2.260596
0.000000,-1.370536
0.000000,0.834539
0.000000,2.066445
0.000000,0.380506
0.000000,-1.663686
0.000000,2.959717
0.000000,-1.494440
0.000000,-0.100453
0.000000,-0.549403
0.000000,0.948625
0.000000,3.035063
0.000001,-1.159644
0.000000,0.563053
0.500000,0.000001
0.000000,3.141593
0.000000,-2.246648
0.000000,-0.467486
0.000000,0.495552
0.000000,2.470974
0.000000,-1.716227
0.000000,0.252930
0.000000,1.695151
0.000000,-0.819699
0.000000,2.121056
0.000000,-1.407460
0.000000,-0.609969
0.000000,-2.662353
1.500000,1.570796
0.000000,-2.332161
0.000000,-0.521964看看幅度谱相位谱,非常标准的对称谱,这表明cfft输出的是对称谱,和FFT的定义一致定点数rfft(不好用)
本着定点数计算更快的想法,我想,fft并不关心实际的幅度的大小,仅仅关心dB,是不是应该用定点数的fft,于是我赤了这依托⑩
下面是定点数fft测试的代码,这里应用了32定点正弦余弦合成,然后fft,发现
q31_t sig[sig_len] = {0};
q31_t fft_out[sig_len] = {0};
q31_t fft_mag[sig_len1] = {0};
q31_t fft_phase[sig_len1] = {0};q31_t i_phase = 0, q_phase = 0;
q31_t i_step = (q31_t)(((q63_t)5 * 0x7fffffff * 2) / sig_len); // 注意精度
q31_t q_step = (q31_t)(((q63_t)8 * 0x7fffffff * 2) / sig_len);for (int i = 0; i sig_len; i++) {sig[i] = (arm_cos_q31(i_phase) 8) + (arm_sin_q31(q_phase) 7);i_phase += i_step; // 自然溢出,完美对应周期延拓q_phase += q_step;
}
arm_rfft_instance_q31 S;
arm_rfft_init_128_q31(S,0, 1);
arm_rfft_q31(S, sig, fft_out);//128点
arm_cmplx_mag_q31(fft_out,fft_mag,sig_len1);//幅度谱
for(int i = 0; i (sig_len1); i++) {//相位谱arm_atan2_q31(fft_out[i*2+1], fft_out[i*2], fft_phase[i]);
}
for(int i = 0; i (sig_len1); i++) {printf("%10d,%10d\n", fft_mag[i], fft_phase[i]);
}
跑出来的频谱是可以看到,这个明明是单边谱的fft,在频谱的另一边居然有两个小峰
然后我把代码中生成正弦函数的代码修改了一下,把正弦数组的幅度改了一下,发现又能用了,小峰看不到了
for (int i = 0; i sig_len; i++) {sig[i] = (arm_cos_q31(i_phase) 12) + (arm_sin_q31(q_phase) 11);i_phase += i_step; // 自然溢出,完美对应周期延拓q_phase += q_step;
}我把信号的幅度再调小,直接全部没有信号了
for (int i = 0; i sig_len; i++) {sig[i] = (arm_cos_q31(i_phase) 15) + (arm_sin_q31(q_phase) 16);i_phase += i_step; // 自然溢出,完美对应周期延拓q_phase += q_step;
}我感觉是dsp库的定点乘法直接把精度损失掉了
大信号有奇怪的镜像小峰
小信号直接出不来信号
我建议是不用定点数的fft
定点数cfft(可以用)
我又试了一下cfft,居然是正常的
18bit到31bit运算都是正常的,小于18bit的数据看不到尖峰
然后换成了16bit,考虑到绝大多数情况下,进行运算的数据都是16bit的,这里贴出来16bit定点数cfft的代码
q15_t数据运算的时候,如果数据仅有6位,输出的结果就有点异常了,所以尽量让输入的数据占据13-15的位宽
q15_t sig[sig_len1] = {0};
q15_t fft_out[sig_len1] = {0};
q15_t fft_mag[sig_len] = {0};
q15_t fft_phase[sig_len] = {0};q15_t i_phase = 0, q_phase = 0;
q15_t i_step = (q15_t)(((q63_t)5 * 0x7fff * 2) / sig_len); // 注意精度
q15_t q_step = (q15_t)(((q63_t)8 * 0x7fff * 2) / sig_len);for (int i = 0; i sig_len; i++) {sig[i*2] = (arm_cos_q15(i_phase) 1) + (arm_sin_q15(q_phase) 1);sig[i*2+1] = 0;i_phase += i_step; // 自然溢出,完美对应周期延拓q_phase += q_step;
}
arm_cfft_instance_q15 S;
arm_cfft_init_q15(S, sig_len);
arm_cfft_q15(S, sig, 0, 1);
arm_cmplx_mag_q15(sig,fft_mag,sig_len);//幅度谱
for(int i = 0; i (sig_len); i++) {//相位谱arm_atan2_q15(sig[i*2+1], sig[i*2], fft_phase[i]);
}
for(int i = 0; i (sig_len); i++) {printf("%10d,%10d\n", fft_mag[i], fft_phase[i]);
}运行速度测试
测试512点fft速度(float32_t float32_t_cfft q31 q15)
包括计算fft以及转成幅度谱的时间
测试结果运算类型消耗时间(us)float32_t cfft
23195.700000float32_t rfft
18373.137500q31_t cfft
6109.600000q15_t cfft
4883.525000
哦,考虑用cfft的时候,还要给数据填充0,测一下整体的时间
测试结果运算类型消耗时间(us)float32_t cfft
23318.925000float32_t rfft
18373.137500q31_t cfft
6444.512500q15_t cfft
5123.850000
测试后还发现,IFFT的速度与正向的速度一致
数字滤波
dsp库提供了直接I型、直接二型转置、格型IIR,FIR有正经的直接型和格型滤波器,还有FIR抽取器
本文档只用最简单最实用的直接I型IIR和直接型FIR
数字滤波器这一节,我按照FIR IIR的章节结构进行编写,方便复现
IIR
概况
CMSIS Dsp库不提供直接型的IIR实例,提供的方法是级联biquad IIR滤波器与lattice型IIR实现可变阶数的IIR计算
有经验的读者会知道,lattice型这种实现的IIR,计算效率很低,仅仅适用于自适应滤波器的应用(只需要检验反射系数是否大于1就可以判断滤波器的稳定性)
所以我们这里仅以直接I型IIR为例,介绍如何联合matlab/python工具,设计2*N阶的数字滤波器
以防读者不清楚级联biquad IIR滤波器是什么,这里结合CMSIS DSP库的图解简单说一下
一个N阶的直接型滤波器的递推公式可以表示为双二阶滤波器的递推公式是
y[n]=b0∗x[n]+b1∗x[n−1]+b2∗x[n−2]+a1∗y[n−1]+a2∗y[n−2]y[n] = b_0 * x[n] + b_1 * x[n-1] + b_2 * x[n-2] + a_1 * y[n-1] + a_2 * y[n-2]y[n]=b0∗x[n]+b1∗x[n−1]+b2∗x[n−2]+a1∗y[n−1]+a2∗y[n−2]在CMSIS dsp库中,为了减少直接型滤波器的有限字长效应,做了一个对象arm_biquad_casd_df1_inst_f32(f32可以换成q31 q15),用来级联这样的双二阶滤波器,实现更长,更大阶数的IIR,等效直接型的IIR这样的实例需要用这样的函数初始化
ARM_DSP_ATTRIBUTE void arm_biquad_cascade_df1_init_f32(arm_biquad_casd_df1_inst_f32 * S,uint8_t numStages,const float32_t * pCoeffs,float32_t * pState)其中
S是实例化对象
numStages是级联的Biquad IIR的个数,IIR的总阶数为2∗numStages2*numStages2∗numStages
pCoeffs是滤波器系数,用一维数组排列,排列规则如下对于f32,q31,f64的计算,数组大小必须为5∗numStages5*numStages5∗numStages,系数排列顺序为
b10,b11,b12,−a11,−a12,b20,b21,b22,−a21,−a22,...{b_{10}, b_{11}, b_{12}, -a_{11}, -a_{12}, b_{20}, b_{21}, b_{22}, -a_{21}, -a_{22}, ...}b10,b11,b12,−a11,−a12,b20,b21,b22,−a21,−a22,...
对于q15的计算,数组大小必须为6∗numStages6*numStages6∗numStages,系数排列顺序为
b10,0,b11,b12,−a11,−a12,b20,0,b21,b22,a21,a22,...{b_{10}, 0, b_{11}, b_{12}, -a_{11}, -a_{12}, b_{20}, 0, b_{21}, b_{22}, a_{21}, a_{22}, ...}b10,0,b11,b12,−a11,−a12,b20,0,b21,b22,a21,a22,... ,就是得填一个0,这是为了满足SIMD运算
pState是中间变量存储的数组,大小必须为4∗numStages4*numStages4∗numStages,用来存储IIR运算的x[n−1],x[n−2],y[n−1],y[n−2]{x[n-1], x[n-2], y[n-1], y[n-2]}x[n−1],x[n−2],y[n−1],y[n−2]
非常容易弄出问题的一个点,IIR系数是归一化的,a0a_0a0默认为1,而且a1a_1a1和a2a_2a2必须取相反数输进去!!!
IIR滤波的主体函数是
ARM_DSP_ATTRIBUTE void arm_biquad_cascade_df1_f32(const arm_biquad_casd_df1_inst_f32 * S,const float32_t * pSrc,float32_t * pDst,uint32_t blockSize)其中S是IIR实例指针,pSrc是源数组指针,pDst是目标数组指针,blocksize是处理数组的大小
定点数和浮点数的滤波函数的接口是一样的
需要特别注意的是pSrc和pDst是可以用同一个数组的
代码
这里以matlab和python代码为例子,假设采样率是48000Hz,设计一个8阶的椭圆IIR滤波器,截止频率2000Hz
matlab可以用自带的滤波器设计工具 filter designer完成,我这的C代码的系数就是用它生成的,但是他总是把通带增益弄得很大,我很不喜欢,还是手写代码更好掌控
fs = 48000;
f0 = 2000;
[z,p,k] = ellip(8,1,80,f0/fs*2,"low");% 椭圆低通滤波器,通带波纹1,阻带波纹-80,截止频率2000Hz
sos = zp2sos(z,p,k,"up");
figure(1)
freqz(sos)num_sections = size(sos, 1);
biquad_params = struct('b0', [], 'b1', [], 'b2', [], 'a1', [], 'a2', [], 'gain', []);for i = 1:num_sections% sos格式: [b0 b1 b2 a0 a1 a2]% 归一化使a0 = 1a0 = sos(i, 4);biquad_params(i).b0 = sos(i, 1) / a0;biquad_params(i).b1 = sos(i, 2) / a0;biquad_params(i).b2 = sos(i, 3) / a0;biquad_params(i).a1 = sos(i, 5) / a0;biquad_params(i).a2 = sos(i, 6) / a0;biquad_params(i).gain = 1; % 初始增益为1,总增益在g中
endfor i = 1:length(biquad_params)% 输出适配dsp库的格式数据fprintf('%f,%f,%f,%f,%f,\n',...biquad_params(i).b0, biquad_params(i).b1, biquad_params(i).b2, ...-biquad_params(i).a1, -biquad_params(i).a2);
end
fprintf('========================================\n\n');为了满足像唐宇坤这样爱用python的孩子,老大这里还提供了一个简单的python代码示例
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt# 参数设置
fs = 48000
f0 = 2000# 椭圆低通滤波器设计
z, p, k = signal.ellip(8, 1, 80, f0/(fs/2), 'low', output='zpk')
sos = signal.zpk2sos(z, p, k)# 绘制频率响应
w, h = signal.sosfreqz(sos, fs=fs)
plt.figure(1)
plt.semilogx(w, 20*np.log10(abs(h)))
plt.grid(True)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude (dB)')
plt.title('Frequency Response')
plt.show()# 提取并输出双二阶系数
num_sections = sos.shape[0]
for i in range(num_sections):# sos格式: [b0 b1 b2 a0 a1 a2]# 归一化使a0 = 1a0 = sos[i, 3]b0 = sos[i, 0] / a0b1 = sos[i, 1] / a0b2 = sos[i, 2] / a0a1 = sos[i, 4] / a0a2 = sos[i, 5] / a0# 输出适配dsp库的格式数据(注意a1和a2取负)print(f'{b0:.10f},{b1:.10f},{b2:.10f},{-a1:.10f},{-a2:.10f},')print('========================================')运行结果
0.000144,-0.000101,0.000144,1.885565,-0.892758,
1.000000,-1.747297,1.000000,1.895284,-0.925868,
1.000000,-1.861904,1.000000,1.907743,-0.962513,
1.000000,-1.888198,1.000000,1.921517,-0.989109,
========================================把这个结果贴在数组里面就能用了
浮点数IIR代码
#define sig_len 512
#define NUM_STAGES 4arm_biquad_casd_df1_inst_f32 S;
float32_t biquadCoeffs[5*NUM_STAGES] = {0.000144,-0.000101,0.000144,1.885565,-0.892758,1.000000,-1.747297,1.000000,1.895284,-0.925868,1.000000,-1.861904,1.000000,1.907743,-0.962513,1.000000,-1.888198,1.000000,1.921517,-0.989109
};
float32_t biquadState[4*NUM_STAGES] = {0};
arm_biquad_cascade_df1_init_f32(S, NUM_STAGES, (float32_t *)biquadCoeffs, (float32_t *)biquadState);
float32_t input[sig_len] = {0},output[sig_len] = {0},fs = 48000.0;
for(int i = 0; i sig_len; i++){input[i] = arm_sin_f32(2*PI*i * 1000/fs) + arm_cos_f32(2*PI*i *2400/fs) + arm_cos_f32(2*PI*i * 4000/fs);
}
arm_biquad_cascade_df1_f32(S, input, output, sig_len);
for(int i = 0; i sig_len; i++){printf("%f\n", output[i]);
}输出结果
看到多个频率混合的正弦信号只剩下一个频率了,中间中断出是两次串口发送的间断点。
IIR本身有瞬态响应的,所以开头有点变形定点数IIR代码
做两手准备,我们还是给出q31_t,q15_t数据类型的IIR代码
需要注意的是初始化函数的输入
先给出q31_t的代码
arm_biquad_cascade_df1_init_q31是31点的定点直接I型IIR初始化函数,注意他的输入参数是
ARM_DSP_ATTRIBUTE void arm_biquad_cascade_df1_init_q31(arm_biquad_casd_df1_inst_q31 * S,uint8_t numStages,const q31_t * pCoeffs,q31_t * pState,int8_t postShift)多了一个postshift参数,这个参数是为了兼容数值大于1的定点数。
因为CMSIS DSP库中,IIR使用的32位定点数数据类型是q311q^{1}_{31}q311,范围是[-1,1]
有时候,IIR计算得出来的参数是大于一的,这个时候没办法转成q311q^{1}_{31}q311输进C代码里面,就需要在matlab/python里面修改小数位,比如原来的小数位是31位,在matlab里面改成30位之后,系数整体的真实数值会缩小2倍(右移1),因此计算出来的y[n]y[n]y[n]也会缩小2倍,这个时候就需要用到postshift参数,把计算结果再重新把计算结果左移1位(正数代表左移)
滤波器生成的matlab代码如下
输出:fs = 48000;
f0 = 2000;
% 设计滤波器
[z,p,k] = ellip(8,0.1,80,f0/fs*2,"low");
sos = zp2sos(z,p,k,"up");
sos_q = zeros(size(sos));
% 提取系数
num_sections = size(sos, 1);
word_length = 32;
fraction_length = 30;
for i = 1:num_sectionsa0 = sos(i,4);sos_q(i,1) = int(fi(sos(i, 1) / a0, true, word_length, fraction_length));sos_q(i,2) = int(fi(sos(i, 2) / a0, true, word_length, fraction_length));sos_q(i,3) = int(fi(sos(i, 3) / a0, true, word_length, fraction_length));sos_q(i,4) = int(fi(1, true, word_length, fraction_length));sos_q(i,5) = int(fi(sos(i, 5) / a0, true, word_length, fraction_length));sos_q(i,6) = int(fi(sos(i, 6) / a0, true, word_length, fraction_length));
end% 绘制频率响应
[h,w] = freqz(sos);
[h_q,w_q] = freqz(sos_q);figure(1)
title('幅频响应')
subplot(2,1,1)
plot(w/pi/2,20*log10(abs(h)),'Color','r','LineWidth',2,'LineStyle','--')
hold on
plot(w_q/pi/2,20*log10(abs(h_q)),'Color','b','LineWidth',2,'LineStyle',':')
legend("目标幅频响应","定点数实现幅频响应")
subplot(2,1,2)
plot(w/pi/2,angle(h),'Color','r','LineWidth',2,'LineStyle','--')
hold on
plot(w_q/pi/2,angle(h_q),'Color','b','LineWidth',2,'LineStyle',':')
legend("目标相频响应","定点数实现相频响应")fprintf('定点格式: Q%d.%d (字长=%d, 小数位=%d)\n', ...word_length - fraction_length, fraction_length, word_length, fraction_length);
fprintf('系数范围: [%f, %f]\n', -2^(word_length-fraction_length-1), 2^(word_length-fraction_length-1)-2^(-fraction_length));
fprintf('\n');
fprintf("结果:\n")
for i = 1:num_sectionsfprintf('%d,%d,%d,%d,%d,\n', ...sos_q(i,1), sos_q(i,2), sos_q(i,3), -sos_q(i,5),-sos_q(i,6));
end定点格式: Q2.30 (字长=32, 小数位=30)
系数范围: [-2.000000, 2.000000]结果:
161226,-68859,161226,1939138878,-879493246,
1073741824,-1799818739,1073741824,1969882589,-932549757,
1073741824,-1963501837,1073741824,2009600067,-997174934,
1073741824,-2001994486,1073741824,2048176666,-1050294896,q31_t的C代码如下,
#define sig_len 512
#define NUM_STAGES 4
arm_biquad_casd_df1_inst_q31 S;
q31_t biquadCoeffs[5*NUM_STAGES] = {161226,-68859,161226,1939138878,-879493246,1073741824,-1799818739,1073741824,1969882589,-932549757,1073741824,-1963501837,1073741824,2009600067,-997174934,1073741824,-2001994486,1073741824,2048176666,-1050294896,
};
q31_t biquadState[4*NUM_STAGES] = {0};
arm_biquad_cascade_df1_init_q31(S, NUM_STAGES, (q31_t *)biquadCoeffs, (q31_t *)biquadState, 1);
q31_t input[sig_len] = {0},output[sig_len] = {0},fs = 48000.0;q31_t phase1 = 0, phase2 = 0, phase3 = 0,step1 = (q63_t)0x80000000*1000/fs, step2 = (q63_t)0x80000000*2400/fs, step3 = (q63_t)0x80000000*4000/fs;
for(int i = 0; i sig_len; i++){input[i] = (arm_sin_q31(phase1)3) + (arm_sin_q31(phase2)4) + (arm_sin_q31(phase3)5);phase1 += step1;phase2 += step2;phase3 += step3;
}
arm_biquad_cascade_df1_q31(S, input, output, sig_len);
for(int i = 0; i sig_len; i++){printf("%d,%d\n", output[i], input[i]);
}输出结果
左边是输出,右边是输入q15_t的C代码如下,注意biquadCoeffs必须是6*NUM_STAGES个数字,第二个系数填充一个零
#define sig_len 512
#define NUM_STAGES 4
arm_biquad_casd_df1_inst_q15 S;
q15_t biquadCoeffs[6*NUM_STAGES] = {2,0,-1,2,29589,-13420,16384,0,-27463,16384,30058,-14230,16384,0,-29961,16384,30664,-15216,16384,0,-30548,16384,31253,-16026
};
q15_t biquadState[4*NUM_STAGES] = {0};
arm_biquad_cascade_df1_init_q15(S, NUM_STAGES, (q15_t *)biquadCoeffs, (q15_t *)biquadState, 1);
q15_t input[sig_len] = {0},output[sig_len] = {0};
q31_t fs = 48000;q15_t phase1 = 0, phase2 = 0, phase3 = 0,step1 = (q63_t)0x8000*1000/fs, step2 = (q63_t)0x8000*3000/fs, step3 = (q63_t)0x8000*4000/fs;
for(int i = 0; i sig_len; i++){input[i] = (arm_sin_q15(phase1)3) + (arm_sin_q15(phase2)2) + (arm_sin_q15(phase3)5);phase1 += step1;phase2 += step2;phase3 += step3;
}
arm_biquad_cascade_df1_q15(S, input, output, sig_len);
for(int i = 0; i sig_len; i++){printf("%d,%d\n", output[i], input[i]);
}输出结果
左边是输出,右边是输入FIR
概况
FIR的初始化代码的输入参数是是实例化、抽头数、系数数组、计算缓存区、处理数组长度,其中数组长度要求是
Coeffs长度:numTaps
State长度:blockSize + numTaps - 1
numTaps长度挺长的,跟输入数组差不多长,这就导致FIR很耗资源,好巧不巧FIR这种滤波器就是得上大阶数,不然不好用
void arm_fir_init_f32 (arm_fir_instance_f32 *S, uint16_t numTaps, float32_t *pCoeffs, float32_t *pState, uint32_t blockSize)
arm_status arm_fir_init_q15 (arm_fir_instance_q15 *S, uint16_t numTaps, q15_t *pCoeffs, q15_t *pState, uint32_t blockSize)
void arm_fir_init_q31 (arm_fir_instance_q31 *S, uint16_t numTaps, q31_t *pCoeffs, q31_t *pState, uint32_t blockSize)
void arm_fir_init_q7 (arm_fir_instance_q7 *S, uint16_t numTaps, q7_t *pCoeffs, q7_t *pState, uint32_t blockSize)代码
浮点数FIR代码
这里不再提供FIR的matlab和python设计代码,累了,直接用matlab的filter designer算了,点几下的事情,还能看误差,需要注意的是这个小app设计40阶FIR的时候会导出41个系数,正常
这里的代码设计的一个40阶的FIR滤波器,汉明窗,具体参数如图左上角点文件,把变量导出到工作区,然后打印出来。由于窗函数设计法得出的滤波器是线性相位的,所以这里看不出来FIR的顺序和CMSIS dsp库里面的输入参数的顺序,按照官网的说明,FIR系数必须是按照顺序b[numTaps−1],h[numTaps−2],h[N−2],...,h[1],h[0]{b[numTaps-1], h[numTaps-2], h[N-2], ..., h[1], h[0]}b[numTaps−1],h[numTaps−2],h[N−2],...,h[1],h[0],而matlab的输出顺序一般是[h[0],h[1],h[2],...,h[N]][h[0], h[1], h[2], ..., h[N]][h[0],h[1],h[2],...,h[N]]。总之就是得注意一下,虽然窗函数设计的FIR滤波器是对称的,万一用到Hilbert变换呢,这个时候就要注意滤波器系数的顺序
#define sig_len 512
#define NUM_STAGES 40
arm_fir_instance_f32 S;
float32_t sig[sig_len],sigout[sig_len],state[sig_len + NUM_STAGES - 1],coeffs[NUM_STAGES] = {
-0.001378,-0.002530,-0.003679,-0.004043,-0.002206,0.003199,0.012242,0.022602,0.029528,0.027366,0.012320,-0.014758,-0.047130,-0.073560,-0.082272,-0.066011,-0.025975,0.027407,0.077586,0.107875,0.107875,0.077586,0.027407,-0.025975,-0.066011,-0.082272,-0.073560,-0.047130,-0.014758,0.012320,0.027366,0.029528,0.022602,0.012242,0.003199,-0.002206,-0.004043,-0.003679,-0.002530,-0.001378
};
float32_t fs = 48000,f1 = 1000,f2 = 1500,f3 = 3876,f4 = 7920;
for(int i = 0; i sig_len; i++){sig[i] = arm_sin_f32(2*PI*f1/fs*i) + arm_sin_f32(2*PI*f2/fs*i) + arm_sin_f32(2*PI*f3/fs*i) + arm_sin_f32(2*PI*f4/fs*i);
}
arm_fir_init_f32(S, NUM_STAGES, coeffs, state, sig_len);
arm_fir_f32(S, sig, sig, sig_len);printf("Filtered signal:\n");
for(int i = 0; i sig_len; i++){printf("%f\n", sig[i]);
}输出结果基本滤除,我不是很喜欢FIR因为阶数要求太高了...四十阶的FIR效果都不咋地,40阶的IIR能毁天灭地了。但是不能否认FIR在线性相位和稳定性上,不可挑剔
定点数FIR代码
没什么特别的,接口和浮点数的一样,不过filter_designer不能直接输出定点数据,得加个matlab转化一下定点数。FIR甚至连输出左移右移的接口都没给,看来跟IIR的接口不一样,不知道如果数值过大怎么办,估计得自己加一个数组整体右移
hn = Num; % Num是filter designer输出的变量
for i = 1:length(hn)hn(i) = int(fi(hn(i),1,32,31));fprintf("%d,",hn(i));
end
freqz(hn,1,512)#define sig_len 512
#define NUM_STAGES 40
arm_fir_instance_q31 S;
q31_t sig[sig_len],sigout[sig_len],state[sig_len + NUM_STAGES - 1],coeffs[NUM_STAGES] = {
-2959417,-5433665,-7899954,-8681870,-4736935,6869835,26288983,48536942,63410362,58767596,26458059,-31692771,-101210454,-157969537,-176678208,-141757258,-55779836,58856511,166614231,231659914,231659914,166614231,58856511,-55779836,-141757258,-176678208,-157969537,-101210454,-31692771,26458059,58767596,63410362,48536942,26288983,6869835,-4736935,-8681870,-7899954,-5433665,-2959417
};
q31_t fs = 48000,f1 = 1000,f2 = 1500,f3 = 3876,f4 = 7920;
q31_t phase1 = 0,phase2 = 0,phase3 = 0,phase4 = 0;
q31_t step1 = (q31_t)((q63_t)0x80000000*f1/fs),step2 = (q31_t)((q63_t)0x80000000*f2/fs),step3 = (q31_t)((q63_t)0x80000000*f3/fs),step4 = (q31_t)((q63_t)0x80000000*f4/fs);for(int i = 0; i sig_len; i++){sig[i] = arm_sin_q31(phase1)/4 + arm_sin_q31(phase2)/8 + arm_sin_q31(phase3)/8 + arm_sin_q31(phase4)/4;phase1 += step1;phase2 += step2;phase3 += step3;phase4 += step4;
}
arm_fir_init_q31(S, NUM_STAGES, coeffs, state, sig_len);
arm_fir_q31(S, sig, sigout, sig_len);printf("Filtered signal:\n");
for(int i = 0; i sig_len; i++){printf("%d,%d\n", sigout[i],sig[i]);
}输出结果依旧加上去q15_t的代码
Q15_t好啊,计算快,跟现成的数据流配合好
arm_fir_instance_q15 S;
q15_t sig[sig_len],sigout[sig_len],state[sig_len + NUM_STAGES - 1],coeffs[NUM_STAGES] = {
-45,-83,-121,-132,-72,105,401,741,968,897,404,-484,-1544,-2410,-2696,-2163,-851,898,2542,3535,3535,2542,898,-851,-2163,-2696,-2410,-1544,-484,404,897,968,741,401,105,-72,-132,-121,-83,-45
};
q31_t fs = 48000,f1 = 1000,f2 = 1500,f3 = 3876,f4 = 7920;
q15_t phase1 = 0,phase2 = 0,phase3 = 0,phase4 = 0;
q15_t step1 = (q15_t)((q31_t)0x8000*f1/fs),step2 = (q15_t)((q31_t)0x8000*f2/fs),step3 = (q15_t)((q31_t)0x8000*f3/fs),step4 = (q15_t)((q31_t)0x8000*f4/fs);for(int i = 0; i sig_len; i++){sig[i] = arm_sin_q15(phase1)/4 + arm_sin_q15(phase2)/8 + arm_sin_q15(phase3)/8 + arm_sin_q15(phase4)/4;phase1 += step1;phase2 += step2;phase3 += step3;phase4 += step4;
}
arm_fir_init_q15(S, NUM_STAGES, coeffs, state, sig_len);
arm_fir_q15(S, sig, sigout, sig_len);
printf("Filtered signal:\n");
for(int i = 0; i sig_len; i++){printf("%d,%d\n", sigout[i],sig[i]);
}计算结果比较M0上的fast_FIR和原版FIR以及直接用卷积函数
突然发现定点数有一个FIR函数arm_fir_fast_q31,文档里说了,这个函数针对M3 M4内核做了 优化,G3507是M0+内核,估计没什么效果
还有一个卷积函数arm_conv_q31,接口是
ARM_DSP_ATTRIBUTE void arm_conv_q31(const q31_t * pSrcA,uint32_t srcALen,const q31_t * pSrcB,uint32_t srcBLen,q31_t * pDst)本身FIR就是卷积实现的,这个函数甚至不需要初始化FIR实例化,也不需要缓冲区,缺点是不能把输入输出作为同一个数组
还是比较了一下他和原版的速率
DL_SYSTICK_resetValue();
arm_fir_q31(S, sig, sigout, sig_len);
tick_end = DL_SYSTICK_getValue();
printf("Time taken to arm_fir_q31 %f us\n",(tick_start - tick_end)/80.0);DL_SYSTICK_resetValue();
arm_fir_fast_q31(S, sig, sigout, sig_len);
tick_end = DL_SYSTICK_getValue();
printf("Time taken to arm_fir_fast_q31 %f us\n",(tick_start - tick_end)/80.0);输出结果
Time taken to arm_fir_q31 14878.025000 us
Time taken to arm_fir_fast_q31 16375.462500 us
Time taken to arm_conv_q31 16648.775000 us好吧,arm_fir_q31胜出!
算个破卷积算这么久....
测速
这里把优化设置成了O0,得到的是比较原始的速度,测试之后发现修改编译优化并不能提速,Dsp库已经力竭了
四个2阶节的IIR滤波器,输入数据为512点,测试不同数据类型的运行速度数据类型运行时间/usfloat32_t
22879.012usq31_t
7915.175usq15_t
2898.375us
40阶的FIR滤波器,输入数据为512点,测试不同数据类型的运行速度数据类型运行时间/usfloat32_t
48269.602usq31_t
16132.425usq15_t
4646.075us
Q15_t好啊,计算快,跟现成的数据流配合好
这一节不提供代码,代码太长了,影响篇幅,都是垃圾代码,没营养
矩阵运算
矩阵运算比较冷门,用的比较少,一般解决简单的机器学习(ML)任务如最小二乘法、逆变识问题,这里放在最后一节
矩阵运算基本数据类型 typedef struct{uint16_t numRows; /** number of rows of the matrix. */uint16_t numCols; /** number of columns of the matrix. */float32_t *pData; /** points to the data of the matrix. */} arm_matrix_instance_f32;/*** @brief Instance structure for the floating-point matrix structure.*/typedef struct{uint16_t numRows; /** number of rows of the matrix. */uint16_t numCols; /** number of columns of the matrix. */float64_t *pData; /** points to the data of the matrix. */} arm_matrix_instance_f64;/*** @brief Instance structure for the Q7 matrix structure.*/typedef struct{uint16_t numRows; /** number of rows of the matrix. */uint16_t numCols; /** number of columns of the matrix. */q7_t *pData; /** points to the data of the matrix. */} arm_matrix_instance_q7;/*** @brief Instance structure for the Q15 matrix structure.*/typedef struct{uint16_t numRows; /** number of rows of the matrix. */uint16_t numCols; /** number of columns of the matrix. */q15_t *pData; /** points to the data of the matrix. */} arm_matrix_instance_q15;/*** @brief Instance structure for the Q31 matrix structure.*/typedef struct{uint16_t numRows; /** number of rows of the matrix. */uint16_t numCols; /** number of columns of the matrix. */q31_t *pData; /** points to the data of the matrix. */} arm_matrix_instance_q31;结构体的成员顾名思义,numRows,numCols,pData分别是行、列以及对应的数组指针
矩阵初始化方法
有现成的初始化函数,如下
void arm_mat_init_f32 (arm_matrix_instance_f32 *S, uint16_t nRows, uint16_t nColumns, float32_t *pData)
void arm_mat_init_q15 (arm_matrix_instance_q15 *S, uint16_t nRows, uint16_t nColumns, q15_t *pData)
void arm_mat_init_q31 (arm_matrix_instance_q31 *S, uint16_t nRows, uint16_t nColumns, q31_t *pData)S是矩阵实例,nRows是行数,nColumns是列数,pData指向二维数组
F12翻开函数一看,其实就是对矩阵实例的成员进行定义
ARM_DSP_ATTRIBUTE void arm_mat_init_q15(arm_matrix_instance_q15 * S,uint16_t nRows,uint16_t nColumns,q15_t * pData)
{/* Assign Number of Rows */S-numRows = nRows;/* Assign Number of Columns */S-numCols = nColumns;/* Assign Data pointer */S-pData = pData;
}所以可以直接定义矩阵结构体,然后对结构体成员定义,如下
q15_t a[4][10] = {0, 18, 28, 38, 48, 58, 68, 78, 88, 98, 108, 118, 128, 138, 148, 158, 168, 178, 188, 198,0, 18, 28, 38, 48, 58, 68, 78, 88, 98, 108, 118, 128, 138, 148, 158, 168, 178, 188, 198};
q15_t b[10][4] = {0, 18, 28, 38, 48, 58, 68, 78, 88, 98, 108, 118, 128, 138, 148, 158, 168, 178, 188, 198,0, 18, 28, 38, 48, 58, 68, 78, 88, 98, 108, 118, 128, 138, 148, 158, 168, 178, 188, 198};
q15_t c[4][4] = {0};arm_matrix_instance_q15 mat_a = {4, 10, (q15_t *)a};
arm_matrix_instance_q15 mat_b = {10, 4, (q15_t *)b};
arm_matrix_instance_q15 mat_c = {4, 4, (q15_t *)c};这意味着,不仅仅需要这个arm_matrix_instance_q15类型的实例化,还需要保证传进这个实例的数组A,B,C不被其他的代码修改,更不能多个矩阵实例共用同一个数组
矩阵运算函数
CMSIS DSP库的矩阵二元运算函数一般的接口都是
arm_status arm_mat_mult_f32 (const arm_matrix_instance_f32 *pSrcA, const arm_matrix_instance_f32 *pSrcB, arm_matrix_instance_f32 *pDst)返回值arm_status返回计算的状态,一般矩阵运算的行列对上号,就不用管返回值
source A B就是参与运算的两个参数 A×B=CA \times B = CA×B=C
加减
函数声明如下,常用的矩阵加减法法运算函数如下,注意保证矩阵的维度一致
这些函数的接口都是A±B=CA ± B = CA±B=C,没有特殊的接口
//加法
arm_status arm_mat_add_f32 (const arm_matrix_instance_f32 *pSrcA, const arm_matrix_instance_f32 *pSrcB, arm_matrix_instance_f32 *pDst)
arm_status arm_mat_add_q15 (const arm_matrix_instance_q15 *pSrcA, const arm_matrix_instance_q15 *pSrcB, arm_matrix_instance_q15 *pDst)
arm_status arm_mat_add_q31 (const arm_matrix_instance_q31 *pSrcA, const arm_matrix_instance_q31 *pSrcB, arm_matrix_instance_q31 *pDst)
//减法
arm_status arm_mat_sub_f32 (const arm_matrix_instance_f32 *pSrcA, const arm_matrix_instance_f32 *pSrcB, arm_matrix_instance_f32 *pDst)
arm_status arm_mat_sub_q15 (const arm_matrix_instance_q15 *pSrcA, const arm_matrix_instance_q15 *pSrcB, arm_matrix_instance_q15 *pDst)
arm_status arm_mat_sub_q31 (const arm_matrix_instance_q31 *pSrcA, const arm_matrix_instance_q31 *pSrcB, arm_matrix_instance_q31 *pDst)乘法
乘法分为矩阵乘法A×B=CA\times B = CA×B=C和常数乘法C∗A=CAC * A = CAC∗A=CA
需要注意的是特别唐的一个问题,arm_mat_mult_fast_q15 和arm_mat_mult_q15 居然有第四个参数,有些教程描述为UNUSE,这个库函数描述为pState,我试了一下,第四个参数填NULL,arm_mat_mult_q15可以使用,但是arm_mat_mult_fast_q15用不了,计算有问题,这个问题我没有解决,好像计算能快个1/6
第二个注意的是,定点数的乘法运算,CMSIS DSP库里面的q15意味着一位符号位,15位小数位,0位整数位,我一般把q15标记成q150q^0_{15}q150,q31标记成q310q^0_{31}q310,按照这个标记方法,我们来模拟一下dsp库进行定点计算的流程
假设a=256,形式是q150,b=512,形式是q150q150∗q150=q301为了保证小数位以及整体的位宽对齐,需要约掉一位整数位,15位小数位,因此c=a×b15,再手动忽略掉多出来的一位整数位c=256∗51215=4假设a = 256, 形式是q^0_{15},b = 512,形式是q^0_{15} \\
q^0_{15} * q^0_{15} = q^1_{30}\\
为了保证小数位以及整体的位宽对齐,需要约掉一位整数位,15位小数位,因此\\
c = a\times b 15,再手动忽略掉多出来的一位整数位\\
c = 256*51215 = 4假设a=256,形式是q150,b=512,形式是q150q150∗q150=q301为了保证小数位以及整体的位宽对齐,需要约掉一位整数位,15位小数位,因此c=a×b15,再手动忽略掉多出来的一位整数位c=256∗51215=4
矩阵乘法的接口
pSrcA A矩阵实例指针,pSrcB B矩阵实例指针,pDst C矩阵实例指针
A×B=CA\times B = CA×B=C
矩阵尺度缩放的接口
对于浮点数乘法,pSrc A矩阵实例指针,scaleFract 缩放因子C,pDst 目标矩阵的地址
对于定点数乘法,pSrc A矩阵实例指针,scaleFract 缩放因子C,shift,计算结果左移的数量,pDst 目标矩阵的地址
注意,pSrc和pDst不能是同一个矩阵,不然会坠机 helicopter helicopter
对于定点数scale乘法,shift参数可以额外将运算结果左移,也可以把这个参数设为负数实现右移
//矩阵乘法
arm_status arm_mat_mult_f32 (const arm_matrix_instance_f32 *pSrcA, const arm_matrix_instance_f32 *pSrcB, arm_matrix_instance_f32 *pDst)
arm_status arm_mat_mult_fast_q15 (const arm_matrix_instance_q15 *pSrcA, const arm_matrix_instance_q15 *pSrcB, arm_matrix_instance_q15 *pDst, q15_t *pState)
arm_status arm_mat_mult_fast_q31 (const arm_matrix_instance_q31 *pSrcA, const arm_matrix_instance_q31 *pSrcB, arm_matrix_instance_q31 *pDst)
arm_status arm_mat_mult_q15 (const arm_matrix_instance_q15 *pSrcA, const arm_matrix_instance_q15 *pSrcB, arm_matrix_instance_q15 *pDst, q15_t *pState CMSIS_UNUSED)
arm_status arm_mat_mult_q31 (const arm_matrix_instance_q31 *pSrcA, const arm_matrix_instance_q31 *pSrcB, arm_matrix_instance_q31 *pDst)//矩阵缩放
arm_status arm_mat_scale_f32 (const arm_matrix_instance_f32 *pSrc, float32_t scale, arm_matrix_instance_f32 *pDst)
arm_status arm_mat_scale_q15 (const arm_matrix_instance_q15 *pSrc, q15_t scaleFract, int32_t shift, arm_matrix_instance_q15 *pDst)
arm_status arm_mat_scale_q31 (const arm_matrix_instance_q31 *pSrc, q31_t scaleFract, int32_t shift, arm_matrix_instance_q31 *pDst)下面给出一段实例代码,以说明scale和mult的使用方法
q15_t a[4][10] = {0, 18, 28, 38, 48, 58, 68, 78, 88, 98, 108, 118, 128, 138, 148, 158, 168, 178, 188, 198,0, 18, 28, 38, 48, 58, 68, 78, 88, 98, 108, 118, 128, 138, 148, 158, 168, 178, 188, 198};
q15_t b[10][4] = {0, 18, 28, 38, 48, 58, 68, 78, 88, 98, 108, 118, 128, 138, 148, 158, 168, 178, 188, 198,0, 18, 28, 38, 48, 58, 68, 78, 88, 98, 108, 118, 128, 138, 148, 158, 168, 178, 188, 198};
q15_t c[4][4] = {0};
q15_t d[4][4] = {0};arm_matrix_instance_q15 mat_a = {4, 10, (q15_t *)a};
arm_matrix_instance_q15 mat_b = {10, 4, (q15_t *)b};
arm_matrix_instance_q15 mat_c = {4, 4, (q15_t *)c};
arm_matrix_instance_q15 mat_d = {4, 4, (q15_t *)d};uint32_t tick_start = 0X01000000;
uint32_t tick_end = 0;
DL_SYSTICK_setPeriod(tick_start);
DL_SYSTICK_resetValue();
arm_mat_mult_q15(mat_a, mat_b, mat_c, NULL);
tick_end = DL_SYSTICK_getValue();
for(int i = 0; i 4; i++) {for(int j = 0; j 4; j++) {printf("%4d ",c[i][j]);}printf("\n");
}
printf("使用arm_dsp库函数 计算4*4矩阵 Time: %.3fus\n",(tick_start-tick_end)/80.0);arm_mat_scale_q15(mat_c, 19, 1, mat_d);
for(int i = 0; i 4; i++) {for(int j = 0; j 4; j++) {printf("%4d ",d[i][j]);}printf("\n");
}输出结果880 970 1060 1150
2480 2770 3060 3350 880 970 1060 1150
2480 2770 3060 3350
使用arm_dsp库函数 计算4*4矩阵 Time: 60.300us27 30 33 35 77 86 95 104 27 30 33 35 77 86 95 104 测试过一次乘法运算的速度,如果是A:10×1,B:1×10,A×B=CA:10\times 1,B:1\times 10,A\times B=CA:10×1,B:1×10,A×B=C的运算基本没有加法的,这个时候,dsp库的效率没有手搓的运算高,其他情况下,dsp库的计算效率应该都会更高
转置与逆
这种单元运算的函数接口依然是pSrc和pDest,直接输就行了
需要注意的是,逆矩阵运算只支持浮点数
//转置
arm_status arm_mat_trans_f32 (const arm_matrix_instance_f32 *pSrc, arm_matrix_instance_f32 *pDst)
arm_status arm_mat_trans_q15 (const arm_matrix_instance_q15 *pSrc, arm_matrix_instance_q15 *pDst)
arm_status arm_mat_trans_q31 (const arm_matrix_instance_q31 *pSrc, arm_matrix_instance_q31 *pDst)//逆
arm_status arm_mat_inverse_f32 (const arm_matrix_instance_f32 *pSrc, arm_matrix_instance_f32 *pDst)float32_t a[4][4] = { 6, 7, 8, 8 ,19, 21, 23, 26 ,6, 7, 8, 8 ,19, 21, 23, 26 };
float32_t b[4][4] = {};
arm_matrix_instance_f32 A_mat = {.numRows = 4,.numCols = 4,.pData = (float32_t *)a
},
B_mat = {.numRows = 4,.numCols = 4,.pData = (float32_t *)b
};
arm_mat_trans_f32(A_mat, B_mat);printf("Matrix A:\n");
for(int i = 0; i 4; i++) {for(int j = 0; j 4; j++) {printf("%4.1f ",a[i][j]);}printf("\n");
}printf("Matrix A's transpose:\n");
for(int i = 0; i 4; i++) {for(int j = 0; j 4; j++) {printf("%4.1f ",b[i][j]);}printf("\n");
}arm_mat_inverse_f32(A_mat, B_mat);
printf("Matrix A's inverse:\n");
for(int i = 0; i 4; i++) {for(int j = 0; j 4; j++) {printf("%.3f ",b[i][j]);}printf("\n");
}输出结果
Matrix A:6.0 7.0 8.0 8.0
19.0 21.0 23.0 26.0 6.0 7.0 8.0 8.0
19.0 21.0 23.0 26.0
Matrix A's transpose:6.0 19.0 6.0 19.0 7.0 21.0 7.0 21.0 8.0 23.0 8.0 23.0 8.0 26.0 8.0 26.0
Matrix A's inverse:
-3.000 1.000 0.000 0.000
2.714 -0.857 0.000 0.000
-1.000 0.000 1.000 0.000
0.000 -1.000 0.000 1.000 参考
https://diversant.sk/CMSIS/Documentation/DSP/html/index.html
https://chat.deepseek.com
.markdown-preview { width: 100%; height: 100%; box-sizing: border-box }
.markdown-preview ul { list-style: disc }
.markdown-preview ul ul { list-style: circle }
.markdown-preview ul ul ul { list-style: square }
.markdown-preview ol { list-style: decimal }
.markdown-preview ol ol, .markdown-preview ul ol { list-style-type: lower-roman }
.markdown-preview ol ol ol, .markdown-preview ol ul ol, .markdown-preview ul ol ol, .markdown-preview ul ul ol { list-style-type: lower-alpha }
.markdown-preview .newpage, .markdown-preview .pagebreak { page-break-before: always }
.markdown-preview pre.line-numbers { position: relative; padding-left: 3.8em; counter-reset: linenumber 0 }
.markdown-preview pre.line-numbers>code { position: relative }
.markdown-preview pre.line-numbers .line-numbers-rows { position: absolute; pointer-events: none; top: 1em; font-size: 100%; left: 0; width: 3em; letter-spacing: -1px; border-right: 1px solid rgba(153, 153, 153, 1); -webkit-user-select: none; -moz-user-select: none; -ms-user-select: none; user-select: none }
.markdown-preview pre.line-numbers .line-numbers-rows>span { pointer-events: none; display: block; counter-increment: linenumber 1 }
.markdown-preview pre.line-numbers .line-numbers-rows>span:before { content: counter(linenumber); color: rgba(153, 153, 153, 1); display: block; padding-right: 0.8em; text-align: right }
.markdown-preview .mathjax-exps .MathJax_Display { text-align: center !important }
.markdown-preview:not([data-for="preview"]) .code-chunk .code-chunk-btn-group { display: none }
.markdown-preview:not([data-for="preview"]) .code-chunk .status { display: none }
.markdown-preview:not([data-for="preview"]) .code-chunk .output-div { margin-bottom: 16px }
.markdown-preview .md-toc { padding: 0 }
.markdown-preview .md-toc .md-toc-details { margin: 0; padding: 0; border: none }
.markdown-preview .md-toc summary.md-toc-link-wrapper { cursor: pointer; list-style-position: outside }
.markdown-preview .md-toc summary.md-toc-link-wrapper::marker { font-size: 0.85em }
.markdown-preview .md-toc .md-toc-link-wrapper .md-toc-link { display: inline; padding: 0.25rem 0 }
.markdown-preview .md-toc .md-toc-link-wrapper .md-toc-link div, .markdown-preview .md-toc .md-toc-link-wrapper .md-toc-link p { display: inline }
.markdown-preview .md-toc .md-toc-link-wrapper.highlighted .md-toc-link { font-weight: 800 }
.markdown-preview .md-toc .md-toc-link-wrapper[data-level="0"] { margin-left: 12px }
.markdown-preview .md-toc .md-toc-link-wrapper[data-level="1"] { margin-left: 28px }
.markdown-preview .md-toc .md-toc-link-wrapper[data-level="2"] { margin-left: 44px }
.markdown-preview .md-toc .md-toc-link-wrapper[data-level="3"] { margin-left: 60px }
.markdown-preview .md-toc .md-toc-link-wrapper[data-level="4"] { margin-left: 76px }
.markdown-preview .md-toc .md-toc-link-wrapper[data-level="5"] { margin-left: 80px }
.markdown-preview a.tag, .markdown-preview a.tag:active, .markdown-preview a.tag:link, .markdown-preview a.tag:visited, .markdown-preview span.tag:not(.token) { display: inline-block; background-color: var(--vscode-badge-background,rgba(127,127,127,.18)); color: var(--vscode-badge-foreground,currentColor); border-radius: 2em; font-size: 0.85em; font-weight: 500; padding: 0.05em 0.55em; line-height: 1.5; white-space: nowrap; text-decoration: none; text-shadow: none; box-shadow: none; transition: background-color 0.12s ease, opacity 0.12s ease }
.markdown-preview a.tag { cursor: pointer }
.markdown-preview a.tag:focus, .markdown-preview a.tag:hover { background-color: var(--vscode-badge-background,rgba(127,127,127,.32)); color: var(--vscode-badge-foreground,currentColor); opacity: 0.85; text-decoration: none; outline: 0 }
.markdown-preview .wikilink-embed-content { border-left: 3px solid var(--vscode-textBlockQuote-border,rgba(128,128,128,.3)); background: var(--vscode-textBlockQuote-background,rgba(128,128,128,.05)); margin: 0.5em 0; padding: 0.5em 0.75em; border-radius: 0 3px 3px 0 }
.markdown-preview .wikilink-embed-content.wikilink-embed-error { border-left-color: var(--vscode-errorForeground,#cf222e); color: var(--vscode-errorForeground,#cf222e); font-style: italic; opacity: 0.8 }
.markdown-preview .wikilink-embed-content h1:first-child, .markdown-preview .wikilink-embed-content h2:first-child, .markdown-preview .wikilink-embed-content h3:first-child, .markdown-preview .wikilink-embed-content h4:first-child, .markdown-preview .wikilink-embed-content h5:first-child, .markdown-preview .wikilink-embed-content h6:first-child { margin-top: 0.25em }
.markdown-preview .wikilink-embed-content .wikilink-embed-content { margin: 0.25em 0; padding: 0.25em 0.5em }
.markdown-preview .block-id { display: none }
.scrollbar-style::-webkit-scrollbar { width: 8px }
.scrollbar-style::-webkit-scrollbar-track { border-radius: 10px; background-color: rgba(0, 0, 0, 0) }
.scrollbar-style::-webkit-scrollbar-thumb { border-radius: 5px; background-color: var(--vscode-scrollbarSlider-background,rgba(150,150,150,.66)); border: 4px solid var(--vscode-scrollbarSlider-background,rgba(150,150,150,.66)); background-clip: content-box }
.scrollbar-style::-webkit-scrollbar-thumb:hover { background-color: var(--vscode-scrollbarSlider-hoverBackground,rgba(120,120,120,.66)) }
html body[for="html-export"]:not([data-presentation-mode]) { position: relative; width: 100%; height: 100%; top: 0; left: 0; margin: 0; padding: 0; overflow: auto }
html body[for="html-export"]:not([data-presentation-mode]) .markdown-preview { position: relative; top: 0; min-height: 100vh }
@media screen and (min-width: 914px) { html body[for="html-export"]:not([data-presentation-mode]) .markdown-preview { padding: 2em calc(50% - 457px + 2em) } }
@media screen and (max-width: 914px) { html body[for="html-export"]:not([data-presentation-mode]) .markdown-preview { padding: 2em } }
@media screen and (max-width: 450px) { html body[for="html-export"]:not([data-presentation-mode]) .markdown-preview { font-size: 14px !important; padding: 1em } }
@media print { html body[for="html-export"]:not([data-presentation-mode]) #sidebar-toc-btn { display: none } }
html body[for="html-export"]:not([data-presentation-mode]) #sidebar-toc-btn { position: fixed; bottom: 8px; left: 8px; font-size: 28px; cursor: pointer; color: inherit; z-index: 99; width: 32px; text-align: center; opacity: 0.4 }
html body[for="html-export"]:not([data-presentation-mode])[html-show-sidebar-toc] #sidebar-toc-btn { opacity: 1 }
html body[for="html-export"]:not([data-presentation-mode])[html-show-sidebar-toc] .md-sidebar-toc { position: fixed; top: 0; left: 0; width: 300px; height: 100%; padding: 32px 0 48px; font-size: 14px; box-shadow: 0 0 4px rgba(150, 150, 150, 0.33); box-sizing: border-box; overflow: auto; background-color: inherit }
html body[for=html-export]:not([data-presentation-mode])[html-show-sidebar-toc] .md-sidebar-toc::-webkit-scrollbar { width: 8px }
html body[for=html-export]:not([data-presentation-mode])[html-show-sidebar-toc] .md-sidebar-toc::-webkit-scrollbar-track { border-radius: 10px; background-color: rgba(0, 0, 0, 0) }
html body[for=html-export]:not([data-presentation-mode])[html-show-sidebar-toc] .md-sidebar-toc::-webkit-scrollbar-thumb { border-radius: 5px; background-color: var(--vscode-scrollbarSlider-background,rgba(150,150,150,.66)); border: 4px solid var(--vscode-scrollbarSlider-background,rgba(150,150,150,.66)); background-clip: content-box }
html body[for=html-export]:not([data-presentation-mode])[html-show-sidebar-toc] .md-sidebar-toc::-webkit-scrollbar-thumb:hover { background-color: var(--vscode-scrollbarSlider-hoverBackground,rgba(120,120,120,.66)) }
html body[for="html-export"]:not([data-presentation-mode])[html-show-sidebar-toc] .md-sidebar-toc a { text-decoration: none }
html body[for="html-export"]:not([data-presentation-mode])[html-show-sidebar-toc] .md-sidebar-toc .md-toc { padding: 0 16px }
html body[for="html-export"]:not([data-presentation-mode])[html-show-sidebar-toc] .md-sidebar-toc .md-toc .md-toc-details { margin: 0; padding: 0; border: none }
html body[for="html-export"]:not([data-presentation-mode])[html-show-sidebar-toc] .md-sidebar-toc .md-toc .md-toc-link-wrapper .md-toc-link { display: inline; padding: 0.25rem 0 }
html body[for="html-export"]:not([data-presentation-mode])[html-show-sidebar-toc] .md-sidebar-toc .md-toc .md-toc-link-wrapper .md-toc-link div, html body[for="html-export"]:not([data-presentation-mode])[html-show-sidebar-toc] .md-sidebar-toc .md-toc .md-toc-link-wrapper .md-toc-link p { display: inline }
html body[for="html-export"]:not([data-presentation-mode])[html-show-sidebar-toc] .md-sidebar-toc .md-toc .md-toc-link-wrapper.highlighted .md-toc-link { font-weight: 800 }
html body[for="html-export"]:not([data-presentation-mode])[html-show-sidebar-toc] .md-sidebar-toc .md-toc .md-toc-link-wrapper[data-level="0"] { padding-left: 0 }
html body[for="html-export"]:not([data-presentation-mode])[html-show-sidebar-toc] .md-sidebar-toc .md-toc .md-toc-link-wrapper[data-level="1"] { padding-left: 16px }
html body[for="html-export"]:not([data-presentation-mode])[html-show-sidebar-toc] .md-sidebar-toc .md-toc .md-toc-link-wrapper[data-level="2"] { padding-left: 32px }
html body[for="html-export"]:not([data-presentation-mode])[html-show-sidebar-toc] .md-sidebar-toc .md-toc .md-toc-link-wrapper[data-level="3"] { padding-left: 48px }
html body[for="html-export"]:not([data-presentation-mode])[html-show-sidebar-toc] .md-sidebar-toc .md-toc .md-toc-link-wrapper[data-level="4"] { padding-left: 64px }
html body[for="html-export"]:not([data-presentation-mode])[html-show-sidebar-toc] .md-sidebar-toc .md-toc .md-toc-link-wrapper[data-level="5"] { padding-left: 80px }
html body[for="html-export"]:not([data-presentation-mode])[html-show-sidebar-toc] .markdown-preview { left: 300px; width: calc(100% - 300px); padding: 2em calc(50% - 457px - 300px / 2); margin: 0; box-sizing: border-box }
@media screen and (max-width: 1274px) { html body[for="html-export"]:not([data-presentation-mode])[html-show-sidebar-toc] .markdown-preview { padding: 2em } }
@media screen and (max-width: 450px) { html body[for="html-export"]:not([data-presentation-mode])[html-show-sidebar-toc] .markdown-preview { width: 100% } }
html body[for="html-export"]:not([data-presentation-mode]):not([html-show-sidebar-toc]) .markdown-preview { left: 50%; transform: translateX(-50%) }
html body[for="html-export"]:not([data-presentation-mode]):not([html-show-sidebar-toc]) .md-sidebar-toc { display: none }
code[class*="language-"], pre[class*="language-"] { color: rgba(51, 51, 51, 1); background: left top; font-family: Consolas, "Liberation Mono", Menlo, Courier, monospace; text-align: left; white-space: pre; word-spacing: normal; word-break: normal; word-wrap: normal; line-height: 1.4; -moz-tab-size: 8; -o-tab-size: 8; tab-size: 8; -webkit-hyphens: none; -moz-hyphens: none; -ms-hyphens: none; hyphens: none }
pre[class*="language-"] { padding: 0.8em; overflow: auto; border-radius: 3px; background: rgba(245, 245, 245, 1) }
:not(pre)>code[class*="language-"] { padding: 0.1em; border-radius: 0.3em; white-space: normal; background: rgba(245, 245, 245, 1) }
.token.blockquote, .token.comment { color: rgba(150, 152, 150, 1) }
.token.cdata { color: rgba(24, 54, 145, 1) }
.token.doctype, .token.macro.property, .token.punctuation, .token.variable { color: rgba(51, 51, 51, 1) }
.token.builtin, .token.important, .token.keyword, .token.operator, .token.rule { color: rgba(167, 29, 93, 1) }
.token.attr-value, .token.regex, .token.string, .token.url { color: rgba(24, 54, 145, 1) }
.token.atrule, .token.boolean, .token.code, .token.command, .token.constant, .token.entity, .token.number, .token.property, .token.symbol { color: rgba(0, 134, 179, 1) }
.token.prolog, .token.selector, .token.tag { color: rgba(99, 163, 92, 1) }
.token.attr-name, .token.class, .token.class-name, .token.function, .token.id, .token.namespace, .token.pseudo-class, .token.pseudo-element, .token.url-reference .token.variable { color: rgba(121, 93, 163, 1) }
.token.entity { cursor: help }
.token.title, .token.title .token.punctuation { font-weight: 700; color: rgba(29, 62, 129, 1) }
.token.list { color: rgba(237, 106, 67, 1) }
.token.inserted { background-color: rgba(234, 255, 234, 1); color: rgba(85, 165, 50, 1) }
.token.deleted { background-color: rgba(255, 236, 236, 1); color: rgba(189, 44, 0, 1) }
.token.bold { font-weight: 700 }
.token.italic { font-style: italic }
.language-json .token.property { color: rgba(24, 54, 145, 1) }
.language-markup .token.tag .token.punctuation { color: rgba(51, 51, 51, 1) }
.language-css .token.function, code.language-css { color: rgba(0, 134, 179, 1) }
.language-yaml .token.atrule { color: rgba(99, 163, 92, 1) }
code.language-yaml { color: rgba(24, 54, 145, 1) }
.language-ruby .token.function { color: rgba(51, 51, 51, 1) }
.language-markdown .token.url { color: rgba(121, 93, 163, 1) }
.language-makefile .token.symbol { color: rgba(121, 93, 163, 1) }
.language-makefile .token.variable { color: rgba(24, 54, 145, 1) }
.language-makefile .token.builtin { color: rgba(0, 134, 179, 1) }
.language-bash .token.keyword { color: rgba(0, 134, 179, 1) }
pre[data-line] { position: relative; padding: 1em 0 1em 3em }
pre[data-line] .line-highlight-wrapper { position: absolute; top: 0; left: 0; background-color: rgba(0, 0, 0, 0); display: block; width: 100% }
pre[data-line] .line-highlight { position: absolute; left: 0; right: 0; padding: inherit 0; margin-top: 1em; background: linear-gradient(90deg, rgba(153, 122, 102, 0.1) 70%, rgba(153, 122, 102, 0)); pointer-events: none; line-height: inherit; white-space: pre }
pre[data-line] .line-highlight:before, pre[data-line] .line-highlight[data-end]:after { content: attr(data-start); position: absolute; top: 0.4em; left: 0.6em; min-width: 1em; padding: 0 0.5em; background-color: rgba(153, 122, 102, 0.4); color: rgba(244, 241, 239, 1); font: bold 65% / 1.5 sans-serif; text-align: center; vertical-align: 0.3em; border-radius: 999px; text-shadow: none; box-shadow: 0 1px rgba(255, 255, 255, 1) }
pre[data-line] .line-highlight[data-end]:after { content: attr(data-end); top: auto; bottom: 0.4em }
.emoji { height: 0.8em }
html body { font-family: "Helvetica Neue", Helvetica, "Segoe UI", Arial, freesans, sans-serif; font-size: 16px; line-height: 1.6; color: rgba(51, 51, 51, 1); background-color: rgba(255, 255, 255, 1); overflow: initial; box-sizing: border-box; word-wrap: break-word }
html body>:first-child { margin-top: 0 }
html body h1, html body h2, html body h3, html body h4, html body h5, html body h6 { line-height: 1.2; margin-top: 1em; margin-bottom: 16px; color: rgba(0, 0, 0, 1) }
html body h1 { font-size: 2.25em; font-weight: 300; padding-bottom: 0.3em }
html body h2 { font-size: 1.75em; font-weight: 400; padding-bottom: 0.3em }
html body h3 { font-size: 1.5em; font-weight: 500 }
html body h4 { font-size: 1.25em; font-weight: 600 }
html body h5 { font-size: 1.1em; font-weight: 600 }
html body h6 { font-size: 1em; font-weight: 600 }
html body h1, html body h2, html body h3, html body h4, html body h5 { font-weight: 600 }
html body h5 { font-size: 1em }
html body h6 { color: rgba(92, 92, 92, 1) }
html body strong { color: rgba(0, 0, 0, 1) }
html body del { color: rgba(92, 92, 92, 1) }
html body a:not([href]) { color: inherit; text-decoration: none }
html body a { color: rgba(0, 136, 204, 1); text-decoration: none }
html body a:hover { color: rgba(0, 163, 245, 1); text-decoration: none }
html body img { max-width: 100% }
html body>p { margin-top: 0; margin-bottom: 16px; word-wrap: break-word }
html body>ol, html body>ul { margin-bottom: 16px }
html body ol, html body ul { padding-left: 2em }
html body ol.no-list, html body ul.no-list { padding: 0; list-style-type: none }
html body ol ol, html body ol ul, html body ul ol, html body ul ul { margin-top: 0; margin-bottom: 0 }
html body li { margin-bottom: 0 }
html body li.task-list-item { list-style: none }
html body li>p { margin-top: 0; margin-bottom: 0 }
html body .task-list-item-checkbox { margin: 0 0.2em 0.25em -1.8em; vertical-align: middle }
html body .task-list-item-checkbox:hover { cursor: pointer }
html body blockquote { margin: 16px 0; font-size: inherit; padding: 0 15px; color: rgba(92, 92, 92, 1); background-color: rgba(240, 240, 240, 1); border-left: 4px solid rgba(214, 214, 214, 1) }
html body blockquote>:first-child { margin-top: 0 }
html body blockquote>:last-child { margin-bottom: 0 }
html body hr { height: 4px; margin: 32px 0; background-color: rgba(214, 214, 214, 1); border: 0 none }
html body table { margin: 10px 0 15px; border-collapse: collapse; border-spacing: 0; display: block; width: 100%; overflow: auto; word-break: keep-all }
html body table th { font-weight: 700; color: rgba(0, 0, 0, 1) }
html body table td, html body table th { border: 1px solid rgba(214, 214, 214, 1); padding: 6px 13px }
html body dl { padding: 0 }
html body dl dt { padding: 0; margin-top: 16px; font-size: 1em; font-style: italic; font-weight: 700 }
html body dl dd { padding: 0 16px; margin-bottom: 16px }
html body code { font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 0.85em; color: rgba(0, 0, 0, 1); background-color: rgba(240, 240, 240, 1); border-radius: 3px; padding: 0.2em 0 }
html body code::after, html body code::before { letter-spacing: -0.2em; content: " " }
html body pre>code { padding: 0; margin: 0; word-break: normal; white-space: pre; background: left top; border: 0 }
html body .highlight { margin-bottom: 16px }
html body .highlight pre, html body pre { padding: 1em; overflow: auto; line-height: 1.45; border: rgba(214, 214, 214, 1); border-radius: 3px }
html body .highlight pre { margin-bottom: 0; word-break: normal }
html body pre code, html body pre tt { display: inline; max-width: initial; padding: 0; margin: 0; overflow: initial; line-height: inherit; word-wrap: normal; background-color: rgba(0, 0, 0, 0); border: 0 }
html body pre code:after, html body pre code:before, html body pre tt:after, html body pre tt:before { content: normal }
html body blockquote, html body dl, html body ol, html body p, html body pre, html body ul { margin-top: 0; margin-bottom: 16px }
html body kbd { color: rgba(0, 0, 0, 1); border-top: 1px solid rgba(214, 214, 214, 1); border-right: 1px solid rgba(214, 214, 214, 1); border-bottom: 2px solid rgba(199, 199, 199, 1); border-left: 1px solid rgba(214, 214, 214, 1); padding: 2px 4px; background-color: rgba(240, 240, 240, 1); border-radius: 3px }
@media print { html body { background-color: rgba(255, 255, 255, 1) } html body h1, html body h2, html body h3, html body h4, html body h5, html body h6 { color: rgba(0, 0, 0, 1); page-break-after: avoid } html body blockquote { color: rgba(92, 92, 92, 1) } html body pre { page-break-inside: avoid } html body table { display: table } html body img { display: block; max-width: 100%; max-height: 100% } html body code, html body pre { word-wrap: break-word; white-space: pre } }