Numpy 的仿制 2

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

上一篇的

《Numpy 仿制 01》

，只实现了简单的功能，离能用还差几个核心函数，这篇文章将继续实现仿制 Numpy 的一些功能。

01 

N维数组的 Slice

Slice 是切片的意思。Numpy 的 Slice 可以切取 N 维数组的任意部分。像在切蛋糕一样。Numpy 的 Slice 格式：[ star:end:step, start:end:step]。想必常用Numpy 的小伙伴应该会比较熟悉。这个格式，Numpy 的内部代码叫做 indicator。

仔细观察，这个种格式，有两种模式：一是截取 A 点到 B 点之间的元素、二是就选择某一个元素。例如：[0:3, 3]，只是要提取第一维元素中 0 到 3 的元素，而第零维提取第三个元素。

那么如果在我们用 C 仿制的 ultra_array 中，这种 slice 该如何实现呢？下面讲解我理解的 Numpy 的 Slice 算法。

• 第一步，Indicator 的设计

为了简化起见，Slice 格式中的 step 就先忽略不管。Slice 最重要的功能还是截取功能。那么根据前面的分析，我们也先来设计一种 Indicator，也是包含两种模式：截断和提取，且还要包含作用的维度信息：

struct _ua_indicator {
 // 作用的维度，例如是作用在第二维。
 int __axis;
 // 选择模式下有用。
 int __picked;
 // 截取模式下起始下标。
 int __start;
 // 截取模式下结束的下标 。
 int __tail;
 // 链接下一个 router
 struct _ua_indicator* next;
};

 

• 第二步，估计截取后的体积

毫无疑问，通过 Slice 截取后，它仍然是一个N维数组，除非它的 Slice 的格式错误了。于是我们即可通过 slice 的 Indicator 来计算出截取后的N维数组的维度信息，然后计算出N维数组的体积。例如 3*3*3 的N维数组，通过 [2,0:1] indicator 截取后便是一个 [1, 3] 的 N 维数组。

• 第三步，计算被选中的元素的地址

这一步无疑是最难啊的。原因是，存储数据的是一个一维数组。通过 Indicator 的分析，然后计算截取范围，最后再去定位。

这里设计两个结构体：

• 1 ua_data_chuck_t，这个对象非常重要，用于记录在计算过程中遇到符合 slice 条件的元素块的地址：

typedef struct _ua_data_chunk ua_data_chunk_t;

struct _ua_data_chunk {
 char* chunk_addr;
 size_t chunk_size;
 struct _ua_data_chunk* next;
};

其中：

chunk_addr 就是符合 slice 条件元素块的首地址。

chunk_size 则记录后续连续的符合 slice 条件的元素块的大小。

• 2  ua_chunk_note_t，是保存 ua_data_chunk_t 的链表：

struct _ua_chunk_note {
 size_t* shape;
 int axis_n;
 ua_data_chunk_t* chunk_map;
};

这个等于是一个 新的 N 维数组的蓝图，在计算完 indicator 后，ua_chunck_note_t 会包含新的 N 维数组的维度信息，以及所有被选中的元素的地址。

02 

计算 slice 的元素地址

在计算 slice 过程，主要的函数是 UArray_indicator_analysis。此函数用于分析 Indicator，计算出 slice 后 N 维数组的维度信息，在分析完维度信息后，便开始使用UArray_survey_chuck_address 计算 slice 的元素地址 ：

int UArray_indicator_analysis(ua_indicator_t* indicator_list, 
u_array_t* arr, 
 ua_chunk_note_t* chunk_note) 
{
 chunk_note->shape = NULL;
 chunk_note->axis_n = 0;
 chunk_note->chunk_map = NULL;
 
 ua_indicator_t* ptr = indicator_list;
 int last_axis = -1;
 // 计算总的维数
 while(ptr != NULL) {
 if (ptr->__picked == -1) (chunk_note->axis_n)++;
 last_axis = ptr->__axis;
 ptr = ptr->next;
 }

 if (last_axis >= arr->axis_n) {
 return -1;
 }

 chunk_note->axis_n = chunk_note->axis_n + UA_axisn(arr) - (last_axis+1);

 if (chunk_note->axis_n > 0) 
 chunk_note->shape = malloc( chunk_note->axis_n * sizeof(size_t) );
 else 
 return -1;

 ptr = indicator_list;
 int axis_index = 0;
 while(ptr != NULL) {
 if (ptr->__picked == -1) {
 int tail = (ptr->__tail<=0?UA_shape_axis(arr, ptr->__axis) + ptr->__tail:ptr->__tail);
 (chunk_note->shape)[axis_index++] = tail - ptr->__start;
 } 
 ptr = ptr->next;
 }

 for (int i = (last_axis+1); i<arr->axis_n; ++i){
 (chunk_note->shape)[axis_index++] = UA_shape_axis(arr, i);
 }

 // -------------------------
 ptr = indicator_list; 
 if (ptr != NULL) {
 UArray_survey_chuck_address(arr, UA_data_ptr(arr), ptr, chunk_note);
 } else {
 chunk_note->chunk_map = UArray_datachunk_create(UA_data_ptr(arr), UA_size(arr));
 }
 return 0; 
}

UArray_survey_chuck_address，此函数用于计算符合 slice 的元素的地址:

int UArray_survey_chuck_address(u_array_t* arr, 
char* chunk_start_from, 
 ua_indicator_t* indicator, 
 ua_chunk_note_t* chunk_note) 
{

 size_t sub_chunk_size = (indicator->__axis < UA_axisn(arr) -1 ? UArray_axis_mulitply(arr, indicator->__axis + 1) : 1) * sizeof(vfloat_t);
 int sub_chunk_number = 1;

 if (indicator->next == NULL) {
 // 最后一个 route nod 
 size_t offset = 0;

 // 计算下一个维度每一个块的大小。

 if (indicator->__picked == -1) {
 sub_chunk_number = (indicator->__tail <= 0 ? UA_shape_axis(arr, indicator->__axis) + indicator->__tail : indicator->__tail) - indicator->__start;
 offset = indicator->__start * sub_chunk_size;
 } else {
 offset = indicator->__picked * sub_chunk_size;
 }
 
 ua_data_chunk_t* new_chunk = UArray_datachunk_create(chunk_start_from + offset, sub_chunk_size * sub_chunk_number);
 UArray_datachunk_addto(&(chunk_note->chunk_map), new_chunk);
 
 } else {

 if (indicator->__picked == -1) {
 // : 的情况
 int tail = (indicator->__tail <= 0 ? UA_shape_axis(arr, indicator->__axis) + indicator->__tail : indicator->__tail);
 for (int i=indicator->__start; i<tail; ++i) {
 char* sub_chunk_start_from = chunk_start_from + i * sub_chunk_size;
 UArray_survey_chuck_address(arr, sub_chunk_start_from, indicator->next, chunk_note);
 }
 } else {
 // picked 的情况
 char* sub_chunk_start_from = chunk_start_from + indicator->__picked * sub_chunk_size;
 UArray_survey_chuck_address(arr, sub_chunk_start_from, indicator->next, chunk_note);
 }
 }
 return 0;
}

在做完 Indicator 分析后，ua_chunck_note_t 会包含新的 N 维数组的维度信息，以及所有被选中的元素的地址。最后将根据计算结果，将数据从原来的 N 维数组中复制到新的 N 维数组即可：

u_array_t UArray_fission_with_indicators(u_array_t* a, ua_indicator_t* indicators) 
{
 // 分析计算 indicator
 ua_chunk_note_t chunk_note;
 UArray_indicator_analysis(indicators, a, &chunk_note);
 u_array_t fission = UArray_create_with_axes_array(chunk_note.axis_n, chunk_note.shape);

 ua_data_chunk_t* ptr = chunk_note.chunk_map;
 char* data_ptr = UA_data_ptr(&fission);
 
 // 
 while( ptr != NULL) {
 // 拷贝每一块选中的元素到新的 N 维数组中去。
 memcpy(data_ptr, ptr->chunk_addr, ptr->chunk_size);
 data_ptr += ptr->chunk_size;
 ptr = ptr->next;

 }
 
 UArray_chunk_note_finalize(&chunk_note);
 // 返回新的 N 维数组
 return fission;
}

03 

测试与对比

python 代码：

import numpy as np

if __name__ == "__main__":
 arr = np.arange(3*3*3).reshape((3,3,3))
 print(arr)
 print("arr[0:2,0:1]")
 print(arr[0:2,0:1])
 print(arr[0:2,0:1].shape)

运行结果：

C 版的代码：

 u_array_t u1 = _UArray3d(3,3,3);
 UA_arange(&u1, 3*3*3);

 u_array_t u2 = UA_slice(&u1, "0:2,0:1");
 
 printf("u1: \n");
 UA_display(&u1);
 printf("\n u2: \n");
 UA_display(&u2);

 UArray_(&u1);
 UArray_(&u2);


运行结果：

结果显示一致。

04

 

 

总结

这个 N 维数组的 slice 是比较复杂的，可惜文笔功力有限，不能流畅第展示整个 slice 算法过程，只能缝上代码，供有心的读者研究。