曾经接触过一个项目，需要提取一段音频的声纹特征向量，用于语音识别。这需要一些数学计算。主要包括快速傅立叶变换、高频滤波、倒梅尔系数的计算。网上参考了一些资料，有用的都是使用 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 算法。
为了简化起见，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 的分析，然后计算截取范围，最后再去定位。
这里设计两个结构体：
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 条件的元素块的大小。
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 算法过程，只能缝上代码，供有心的读者研究。
评论