Imitation of numpy 2

I have been in contact with a project , It is necessary to extract the voiceprint feature vector of an audio , For speech recognition . This requires some mathematical calculation . It mainly includes fast Fourier transform 、 High frequency filtering 、 Calculation of inverse Mel coefficient . I refer to some materials on the Internet , What is useful is to use  Python Realized , Among them, there are  Numpy  The highest frequency of use . But our project needs to extract voiceprints from mobile phones , Only use  C  To achieve audio voiceprint extraction . The difficulty is to  C  Write a set of  Numpy  The function of （ Part of the ）, And then emulate  Python  Code writing  C  Version of voiceprint extraction .

From the last article

《Numpy  Copy  01》

, Only simple functions are realized , There are still several core functions to be used , This article will continue to realize imitation  Numpy  Some functions of .

01 

N Dimension group  Slice

Slice  It means slicing .Numpy  Of  Slice  You can cut  N  Any part of the dimension group . It's like cutting a cake .Numpy  Of  Slice  Format ：[ star:end:step, start:end:step]. Must be commonly used Numpy  My little friend should be familiar with . This format ,Numpy  The internal code of is called  indicator.

Observe carefully , This format , There are two patterns ： First, intercept  A  Point to  B  Elements between points 、 The second is to choose an element . for example ：[0:3, 3], Just extract the first dimension elements  0  To  3  The elements of , And the zero dimension extracts the third element .

So if we use  C  Imitated  ultra_array  in , such  slice  How to achieve it ？ Here is what I understand  Numpy  Of  Slice  Algorithm .

• First step ,Indicator  The design of the

To simplify ,Slice  In the format  step  Just ignore it first .Slice  The most important function is the interception function . Then according to the previous analysis , Let's also design a  Indicator, It also includes two modes ： Truncation and extraction , And also include the dimension information of the role ：

struct _ua_indicator {
 //  The dimension of action , For example, it acts on the second dimension .
 int __axis;
 //  Useful in select mode .
 int __picked;
 //  Start subscript in interception mode .
 int __start;
 //  Subscript ending in intercept mode  .
 int __tail;
 //  Link to next  router
 struct _ua_indicator* next;
};

 

• The second step , Estimate the volume after interception

without doubt , adopt  Slice  After interception , It's still a N Dimension group , Unless it's  Slice  The format of is wrong . So we can pass  slice  Of  Indicator  To calculate the intercepted N Dimension information of dimension group , And then calculate N Volume of dimension group . for example  3*3*3  Of N Dimension group , adopt  [2,0:1] indicator  After interception, there is a  [1, 3]  Of  N  Dimension group .

• The third step , Calculate the address of the selected element

This step is undoubtedly the most difficult . as a result of , The data is stored in a one-dimensional array . adopt  Indicator  Analysis of , Then calculate the interception range , Finally, locate .

Here we design two structures ：

• 1 ua_data_chuck_t, This object is very important , It is used to record the meeting of  slice  Address of element block of condition ：

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;
};

among ：

chunk_addr  It's in line with  slice  The first address of the conditional element block .

chunk_size  Then record the subsequent continuous conformity  slice  The size of the element block of the condition .

• 2  ua_chunk_note_t, Is save  ua_data_chunk_t  The linked list of ：

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

This is equivalent to a   new  N  Blueprint of dimension group , At the end of the calculation  indicator  after ,ua_chunck_note_t  Will include new  N  Dimension information of dimension group , And the addresses of all selected elements .

02 

Calculation  slice  Element address of

In the calculation  slice  The process , The main function is  UArray_indicator_analysis. This function is used to analyze  Indicator, To calculate the  slice  after  N  Dimension information of dimension group , After analyzing the dimension information , I started using UArray_survey_chuck_address  Calculation  slice  Element address of  ：

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;
 //  Calculate the total dimension
 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, This function is used to calculate compliance  slice  The address of the element of :

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) {
 //  the last one  route nod 
 size_t offset = 0;

 //  Calculate the size of each block in the next dimension .

 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) {
 // :  The situation of
 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  The situation of
 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;
}

Finishing  Indicator  After analysis ,ua_chunck_note_t  Will include new  N  Dimension information of dimension group , And the addresses of all selected elements . Finally, according to the calculation results , Change the data from the original  N  Copy the dimension group to the new  N  Dimension array is enough ：

u_array_t UArray_fission_with_indicators(u_array_t* a, ua_indicator_t* indicators) 
{
 //  Analysis and calculation  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) {
 //  Copy each selected element to a new  N  Go to dimension group .
 memcpy(data_ptr, ptr->chunk_addr, ptr->chunk_size);
 data_ptr += ptr->chunk_size;
 ptr = ptr->next;

 }
 
 UArray_chunk_note_finalize(&chunk_note);
 //  Return to new  N  Dimension group
 return fission;
}

03 

Testing and comparison

python  Code ：

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)

Running results ：

C  Version of the code ：

 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);


Running results ：

The results are consistent .

04

 

 

summary

This  N  Dimension group  slice  It's more complicated , Unfortunately, the writing skills are limited , Can't show the whole  slice  The algorithm process , Only code can be sewn , For interested readers to study .