Numpy 的研究仿制 1

曾经接触过一个项目，需要提取一段音频的声纹特征向量，用于语音识别。这需要一些数学计算。主要包括快速傅立叶变换、高频滤波、倒梅尔系数的计算。网上参考了一些资料，有用的都是使用 Python实现的，其中又以 Numpy 的使用频率最高。但我们的项目需要在手机端提取声纹，只能用 C 来实现音频的声纹提取。这其中的难点就是要 C 写一套 Numpy 的功能（一部分），然后仿效 Python 代码写出 C 版本的声纹提取。

本文主要记录笔者为了在实现声纹提取，而不得已使用 C 仿制 Numpy 的一部分功能。记录在研究 Numpy 过程中遇到一些有趣的过程。

01 

 

N维数组的分析

NumPy 是使用 Python 进行科学计算的基础软件包。NumPy 最重要的一个特点是其 N 维数组对象 ndarray，它是一系列同类型数据的集合，以 0 下标为开始进行集合中元素的索引。ndarray 对象是用于存放同类型元素的多维数组。ndarray 中的每个元素在内存中都有相同存储大小的区域。用 C 来仿制一个 Numpy 该如何设计呢？

关键还是这个 N 维数组。N 是不确定，那么设计对象的时候，内部数组就不是固定的。这对于一个静态语言是一个挑战。但有一点是确定的，那是 N 维数组的元素数量是确定的，元素数量确定，就意味内存空间大小确定。例如 3*5*5 的三维数组，如果是 float 型，那就是 3*5*5*sizeof(float) = 75 * 4 = 300 byte。于是对内我们设计一个大小一致的连续内存块。对外通过接口，将其还原成三维数组。

先给这个用 C 写的冒牌 Numpy 起一个霸气的名字 --- ultra_array，原型如下：

struct _u_array {
 char *start[2];
 int axis_n;
};

对就是这么简洁大气。start 是一个指针数组。总共就只有两位。这两个指针数组，分别指向两个数组，start[0] 是指向存储维度信息的数组指针。例如 3*5*5 则是 [3, 5, 5]，start[1] 则就是指向各个存储数据的内存块指针。而 axis_n 则表明这个这个数组有维度。例如 3 维就是 3。4 维就是 4。axis_n 确认了 start[0] 的边界。而 start[0] 又确认了 start[1] 的边界。

02

 

多维数组元素访问

如何访问这个 N 维数组中的元素呢？N 维数组中数据存放在一维数组中。那么我们在接口处输入元素的坐标，则需要通一个矩阵来转换得到一维数组的坐标，例如有一个 3*4*5*6 数组，我们要访问（2，2，2，3）这个坐标的数据，那么我们先要的计算出这个这个 3*3*5*6 这个转换矩阵。

大家是不是以为转换矩阵就是（3，4，5，6）。错！如果以为是这个，那就犯了经验主义误。它的转换矩阵是要从倒数第二维开始，每一维乘后一维到最后一维为止，得出的矩阵就是将 N 维数组转一维数组的转换矩阵，而最后一维用1来代替。即（3*4*5*6）数组的转换矩阵是 [ 4*5*6，5*6， 6，1 ] => [ 120, 30, 6, 1 ]。

那么（2，2，2，3）这个坐标转换成一维的坐标就是 （2，2，2，3）dot (120，30，6，1)

T 

=> 2*120 + 2*30 + 2*6 + 3*1 => 315。代码实现如下：

static size_t
__xd_coord_to_1d_offset(size_t coord[], size_t axes[], int axis_n) {

 size_t offset = 0, axis_mulitply;
 for (int i=0; i<axis_n; ++i) {
 size_t co = coord[i];
 axis_mulitply = __axis_mulitply(axes, axis_n, i+1);
 offset += co * axis_mulitply;
 }
 return offset;
}

那么一维坐标又如何成为 N 维坐标呢？倒数第一维开始，我们需要用一维的坐标值，除以倒数第二维乘到最后一维的值，得出的商作为当前维数的坐标，得出的余数作下一个维度的总值，用总值再去除以下一个维度到最后一个维度的乘积，一直到最后一维。例如刚刚我们算出来的一维坐标是 315，那么根据以上推算是：

315 / （4*5*6）= 2 余 75

75   / （5*6）   = 2 余 15 

15   / （6）       = 2 余 3   

 3    /    1          = 3     

于是得出坐标是 [2，2，2，3]。代码实现是：

static void
__1d_offset_to_xd_coord( size_t offset, size_t axes[], int axis_n, size_t coord[])
{
 size_t div, mod, i, axis_mulitply, middle_value;
 middle_value = offset;
 for(i=0; i<axis_n-1; ++i) {
 axis_mulitply = __axis_mulitply(axes, axis_n, i+1);
 div = middle_value / axis_mulitply;
 mod = middle_value % axis_mulitply;
 coord[i] = div;
 middle_value = mod;
 }
 coord[i] = mod;
 return;
}

03 

代码实现

• 初始化

/**
 * 输入维度数量，例如 3 维
 * 输入每一个维度，例如 [3, 3, 3]
 */
u_array_t UArray_create(int axis_n, size_t shape[]) 
{
 if (axis_n >= 0) {
 u_array_t n_array;
 n_array.axis_n = axis_n;
 start[0] = __alloc_shape(axis_n, shape);
 start[1] = __alloc_data(__axis_mulitply(shape, axis_n, 0));

 return n_array;
 }
 return ua_unable;
}

• 加载数据

u_array_t* UArray_load(u_array_t* arr, vfloat_t data[])
{
 size_t size_arr = UA_size(arr);
 vfloat_t* ptr = UA_data_ptr(arr);
 memcpy(ptr, data, size_arr); 
 return arr;
}

• 访问数据

float UArray_get(u_array_t* arr, ...) 
{
 va_list valist;
 va_start(valist, arr);
 size_t coord[UA_axisn(arr)];
 for (int i=0; i<UA_axisn(arr); ++i) {
 coord[i] = va_arg(valist, size_t);
 }
 va_end(valist);
 size_t offset = UA_cover_coordinate(arr, coord);
 return ((float*)(UA_data_ptr(arr)))[offset];
}

void UArray_set(u_array_t* arr, ...)
{
 va_list valist;
 va_start(valist, arr);
 size_t coord[UA_axisn(arr)];
 vfloat_t value;
 for (int i=0; i<UA_axisn(arr); ++i) {
 coord[i] = va_arg(valist, size_t);
 }
 value = va_arg(valist, double);
 va_end(valist);
 size_t offset = UA_cover_coordinate(arr, coord);
 ((float*)(UA_data_ptr(arr)))[offset] = value;
 return;
}

04 

测试

int main()
{
 // 定义一个 3 维的 ultra_array
 u_array_t arr3 = UArray3d(2, 3, 4);
 // 填入从 0 到 23 的数字。
 UA_arange(&arr3, 2*3*4);
 // 获取
 float v = UA_get(&arr3, 1, 2, 3);
 // v == 23
 UA_set(&arr3, 1, 2, 3, 5.5);
 v = UA_get(&arr3, 1, 2, 3);
 // v == 5.5
 return 0;
}

至此一个简单的 C 版的多维数组就实现了。以上代码均源自：

https://github.com/zuweie/boring-code/tree/main/src/ultra_array

完！