Numpy 的研究仿制 1

作者：祖维

2022 年 6 月 29 日
本文字数：2654 字
阅读完需：约 9 分钟

曾经接触过一个项目，需要提取一段音频的声纹特征向量，用于语音识别。这需要一些数学计算。主要包括快速傅立叶变换、高频滤波、倒梅尔系数的计算。网上参考了一些资料，有用的都是使用 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;}