我们在嵌入式上跑矩阵运算时候,会遇到这样一个问题。
假设将矩阵设置成N*N维的二维数组后,我们想求两个矩阵相乘,那就需要按照矩阵计算规则编写矩阵相乘函数,这样的话4*4矩阵得编一个,5*5矩阵又得编一个,要求逆还得编一个,求行列式还得编。
自己写的函数代码效率容易低,将导致本来要跑在单片机上的算法,难达到想象计算速度。
这篇教程将教会你如何使用arm内核库的矩阵计算函数,让你降低代码编写难度还能提高运算效率。据笔者所知,M4内核自带DSP库。
1.创建新模板工程(这里使用的是keil创建工程)
创建一个对应MCU的新模板工程(添加必要的库引用和创建.C文件,详情略)
2.引入arm内核的DSP库文件(含矩阵计算库)

勾选CMSIS下的DSP,确认后将添加库引用至工程

3.添加矩阵函数库
在main函数添加引用,”arm_math.h”

还需在”option of target”中”C/C++”中的define栏添加“,ARM_MATH_CM4”以声明开启arm_math库

再到”option of target”中”target”勾选一下”Use MicroLIB”

4.矩阵计算
//声明大小为NN的浮点一维数组NN=N*N,此处将矩阵转换成一维数组
float32_t Data1[NN] = {};
float32_t Data2[NN] = {};
float32_t Data3[NN];
//声明ARM矩阵类型指针
arm_matrix_instance_f32 Matrix_data1;//再设立
arm_matrix_instance_f32 Matrix_data2;
arm_matrix_instance_f32 Matrix_data3;
//将一维数组地址赋予给ARM矩阵指针
arm_mat_init_f32(&Matrix_data1,N,N,(float32_t *)Data1);
arm_mat_init_f32(&Matrix_data2,N,N,(float32_t *)Data2);
arm_mat_init_f32(&Matrix_data3,N,N,(float32_t *)Data3);
//使用矩阵相乘函数arm_mat_mult_f32
arm_mat_mult_f32(&Matrix_data1,&Matrix_data2,&Matrix_data3);
总结分析
在ARM内核的单片机中, ARM提供了矩阵运算库”arm_math”。在设置好库引用后,需存入用一维数组表示的矩阵值,再将矩阵的一维数组传递给ARM_matrix的指针,在声明函数中将设置矩阵的行列数,再调用库中的各类运算指令就可。
Ps:在本人stm32单片机上测试,arm自带函数只能计算至16维以内(不含16维)的矩阵计算,大于等于16维矩阵单片机编译没报错,但单片机不正常运行。小编水平有限这个存在疑惑,欢迎大佬们交流。
附图
经测试80M的系统时钟下,两个15维矩阵相乘时间约430us,两个4维矩阵相乘时间仅需14us,这计算速度还是很让人满意的。

附上代码(注释部分为arm官方例程,包括相乘求逆求行列式)
1#include "sys.h" 2#include "delay.h" 3#include "usart.h" 4#include "arm_math.h" 5#include<stdio.h> 6#include<stdlib.h> 7 8#define N 4 9#define NN 16 10 11float32_t a = 100; 12float32_t Data1[NN] = {}; 13float32_t Data2[NN] = {}; 14float32_t Data3[NN]; 15 16arm_matrix_instance_f32 Matrix_data1; 17arm_matrix_instance_f32 Matrix_data2; 18arm_matrix_instance_f32 Matrix_data3; 19 20//const float32_t B_f32[4] = 21//{ 22// 782.0, 7577.0, 470.0, 4505.0 23//}; 24 25///* -------------------------------------------------------------------------------- 26//* Formula to fit is C1 + C2 * numTaps + C3 * blockSize + C4 * numTaps * blockSize 27//* -------------------------------------------------------------------------------- */ 28 29//const float32_t A_f32[16] = 30//{ 31// /* Const, numTaps, blockSize, numTaps*blockSize */ 32// 1.0, 32.0, 4.0, 128.0, 33// 1.0, 32.0, 64.0, 2048.0, 34// 1.0, 16.0, 4.0, 64.0, 35// 1.0, 16.0, 64.0, 1024.0, 36//}; 37 38 39///* ---------------------------------------------------------------------- 40//* Temporary buffers for storing intermediate values 41//* ------------------------------------------------------------------- */ 42///* Transpose of A Buffer */ 43//float32_t AT_f32[16]; 44///* (Transpose of A * A) Buffer */ 45//float32_t ATMA_f32[16]; 46///* Inverse(Transpose of A * A) Buffer */ 47//float32_t ATMAI_f32[16]; 48///* Test Output Buffer */ 49//float32_t X_f32[4]; 50 51///* ---------------------------------------------------------------------- 52//* Reference ouput buffer C1, C2, C3 and C4 taken from MATLAB 53//* ------------------------------------------------------------------- */ 54//const float32_t xRef_f32[4] = {73.0, 8.0, 21.25, 2.875}; 55 56//float32_t snr; 57 58int i = 0,t = 0,j = 0; 59 60int main(void) 61{ 62 HAL_Init(); //初始化HAL库 63 SystemClock_Config(); //初始化系统时钟为80M 64 __HAL_RCC_GPIOE_CLK_ENABLE(); //开启GPIOE时钟 65 uart_init(115200); //初始化串口,速率为115200 66 while(1) 67 { 68 for(i=0;i<NN;i++) 69 Data1[i]=rand()/(double)(RAND_MAX/100); 70 for(i=0;i<NN;i++) 71 Data2[i]=rand()/(double)(RAND_MAX/100); 72 arm_mat_init_f32(&Matrix_data1,N,N,(float32_t *)Data1); 73 arm_mat_init_f32(&Matrix_data2,N,N,(float32_t *)Data2); 74 arm_mat_init_f32(&Matrix_data3,N,N,(float32_t *)Data3); 75 arm_mat_mult_f32(&Matrix_data1,&Matrix_data2,&Matrix_data3); 76 77 for(i=0;i<N;i++) 78 { 79 for(j=0;j<N;j++) 80 printf("%f ",Data3[i*N+j]); 81 printf("\r\n"); 82 } 83// arm_matrix_instance_f32 A; /* Matrix A Instance */ 84// arm_matrix_instance_f32 AT; /* Matrix AT(A transpose) instance */ 85// arm_matrix_instance_f32 ATMA; /* Matrix ATMA( AT multiply with A) instance */ 86// arm_matrix_instance_f32 ATMAI; /* Matrix ATMAI(Inverse of ATMA) instance */ 87// arm_matrix_instance_f32 B; /* Matrix B instance */ 88// arm_matrix_instance_f32 X; /* Matrix X(Unknown Matrix) instance */ 89 90// uint32_t srcRows, srcColumns; /* Temporary variables */ 91// arm_status status; 92 93// /* Initialise A Matrix Instance with numRows, numCols and data array(A_f32) */ 94// srcRows = 4; 95// srcColumns = 4; 96// arm_mat_init_f32(&A, srcRows, srcColumns, (float32_t *)A_f32); 97// 98// /* Initialise Matrix Instance AT with numRows, numCols and data array(AT_f32) */ 99// srcRows = 4; 100// srcColumns = 4; 101// arm_mat_init_f32(&AT, srcRows, srcColumns, AT_f32); 102// 103 104// /* calculation of A transpose */ 105// status = arm_mat_trans_f32(&A, &AT); 106 107//// 108// for(i=0;i<16;i++) 109// printf("%f,",AT_f32[i]); 110// printf("\r\n"); 111 112// /* Initialise ATMA Matrix Instance with numRows, numCols and data array(ATMA_f32) */ 113// srcRows = 4; 114// srcColumns = 4; 115// arm_mat_init_f32(&ATMA, srcRows, srcColumns, ATMA_f32); 116 117// /* calculation of AT Multiply with A */ 118// status = arm_mat_mult_f32(&AT, &A, &ATMA); 119 120// /* Initialise ATMAI Matrix Instance with numRows, numCols and data array(ATMAI_f32) */ 121// srcRows = 4; 122// srcColumns = 4; 123// arm_mat_init_f32(&ATMAI, srcRows, srcColumns, ATMAI_f32); 124 125// /* calculation of Inverse((Transpose(A) * A) */ 126// status = arm_mat_inverse_f32(&ATMA, &ATMAI); 127 128// /* calculation of (Inverse((Transpose(A) * A)) * Transpose(A)) */ 129// status = arm_mat_mult_f32(&ATMAI, &AT, &ATMA); 130 131// /* Initialise B Matrix Instance with numRows, numCols and data array(B_f32) */ 132// srcRows = 4; 133// srcColumns = 1; 134// arm_mat_init_f32(&B, srcRows, srcColumns, (float32_t *)B_f32); 135 136// /* Initialise X Matrix Instance with numRows, numCols and data array(X_f32) */ 137// srcRows = 4; 138// srcColumns = 1; 139// arm_mat_init_f32(&X, srcRows, srcColumns, X_f32); 140 141// /* calculation ((Inverse((Transpose(A) * A)) * Transpose(A)) * B) */ 142// status = arm_mat_mult_f32(&ATMA, &B, &X); 143 144//// /* Comparison of reference with test output */ 145//// snr = arm_snr_f32((float32_t *)xRef_f32, X_f32, 4); 146 147 } 148}