footPDR.c 23 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930
  1. #include <stdio.h>
  2. #include <string.h>
  3. #include <stdlib.h>
  4. #include <math.h>
  5. #include <stdint.h>
  6. #define ZUPT_threshold 0.81
  7. #define SIGMA 0.01
  8. #define SIGMA_V 0.01
  9. #define PI 3.1416
  10. #define RUN 1
  11. #define STAND 0
  12. //当地的重力加速度
  13. float g = 9.788f;
  14. float dt = 0.01f;
  15. float P[81], acc_n[3];
  16. float Temporary_array1[9], Temporary_array2[9];
  17. float K[27], P_prev[81], delta_x[9];
  18. float C[9], C_prev[9];
  19. float vel_n[3], pos_n[3];
  20. float last_pos_n[3];
  21. float pos_offset[3];
  22. int frame_index = 0;
  23. int stand_num = 0;
  24. float gyr_norm_window[10];
  25. float gyr_z_window[10];
  26. float gyr_extreme[6];
  27. float gyr_mean[3];
  28. float num_peak;
  29. float acc_mean[3];
  30. float gyrBias[3];
  31. float last_gyr[3];
  32. float last_acc[3];
  33. float last_vel_n[3];
  34. float accSum;
  35. float accSize;
  36. int press_data[10];
  37. int ZUPT_STATUS;
  38. int HAS_RESULT;
  39. int IS_DOWN;
  40. int press_index;
  41. float acc_shoes[3];
  42. float last_acc_shoes[3];
  43. int last_move_index;
  44. int down_pass1;
  45. int down_pass2;
  46. float theta;
  47. //last_stage:0 为 走路状态;
  48. //last_stage:1 为 静止状态
  49. int last_stage;
  50. unsigned int step_count;
  51. unsigned int step_distance;
  52. enum _CMD_MOTION{
  53. MOTION_STOP = 0,
  54. MOTION_RUN = 1,
  55. MOTION_JUMP = 2,
  56. MOTION_DOWN = 3,
  57. MOTION_LEFT = 4,
  58. MOTION_RIGHT = 5,
  59. MOTION_FRONT = 6,
  60. MOTION_BACK = 7,
  61. };
  62. void cal_step_data(void)
  63. {
  64. step_count = step_count + 1;
  65. float step_length = sqrt((last_pos_n[0] - pos_n[0])*(last_pos_n[0] - pos_n[0]) + (last_pos_n[1] - pos_n[1])*(last_pos_n[1] - pos_n[1]));
  66. if(step_length > 3.0f)
  67. {
  68. step_length = 1.5f;
  69. }
  70. step_distance = step_distance + (uint32_t)(step_length * 100.0f);
  71. }
  72. uint32_t get_step_count(void)
  73. {
  74. return step_count;
  75. }
  76. uint32_t get_step_length(void)
  77. {
  78. return step_distance;
  79. }
  80. void invert3x3(float * src, float * dst)
  81. {
  82. float det;
  83. /* Compute adjoint: */
  84. dst[0] = +src[4] * src[8] - src[5] * src[7];
  85. dst[1] = -src[1] * src[8] + src[2] * src[7];
  86. dst[2] = +src[1] * src[5] - src[2] * src[4];
  87. dst[3] = -src[3] * src[8] + src[5] * src[6];
  88. dst[4] = +src[0] * src[8] - src[2] * src[6];
  89. dst[5] = -src[0] * src[5] + src[2] * src[3];
  90. dst[6] = +src[3] * src[7] - src[4] * src[6];
  91. dst[7] = -src[0] * src[7] + src[1] * src[6];
  92. dst[8] = +src[0] * src[4] - src[1] * src[3];
  93. /* Compute determinant: */
  94. det = src[0] * dst[0] + src[1] * dst[3] + src[2] * dst[6];
  95. /* Multiply adjoint with reciprocal of determinant: */
  96. det = 1.0f / det;
  97. dst[0] *= det;
  98. dst[1] *= det;
  99. dst[2] *= det;
  100. dst[3] *= det;
  101. dst[4] *= det;
  102. dst[5] *= det;
  103. dst[6] *= det;
  104. dst[7] *= det;
  105. dst[8] *= det;
  106. }
  107. void multiply3x3(float *a, float *b, float *dst)
  108. {
  109. dst[0] = a[0] * b[0] + a[1] * b[3] + a[2] * b[6];
  110. dst[1] = a[0] * b[1] + a[1] * b[4] + a[2] * b[7];
  111. dst[2] = a[0] * b[2] + a[1] * b[5] + a[2] * b[8];
  112. dst[3] = a[3] * b[0] + a[4] * b[3] + a[5] * b[6];
  113. dst[4] = a[3] * b[1] + a[4] * b[4] + a[5] * b[7];
  114. dst[5] = a[3] * b[2] + a[4] * b[5] + a[5] * b[8];
  115. dst[6] = a[6] * b[0] + a[7] * b[3] + a[8] * b[6];
  116. dst[7] = a[6] * b[1] + a[7] * b[4] + a[8] * b[7];
  117. dst[8] = a[6] * b[2] + a[7] * b[5] + a[8] * b[8];
  118. }
  119. void multiply3x1(float *a, float *b, float *dst)
  120. {
  121. dst[0] = a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
  122. dst[1] = a[3] * b[0] + a[4] * b[1] + a[5] * b[2];
  123. dst[2] = a[6] * b[0] + a[7] * b[1] + a[8] * b[2];
  124. }
  125. void init_attitude_matrix(float *C, float *acc)
  126. {
  127. float pitch = asin(-acc[0] / g);
  128. float roll = atan2(acc[1], acc[2]);
  129. float pitch_sin = sin(pitch);
  130. float pitch_cos = cos(pitch);
  131. float roll_sin = sin(roll);
  132. float roll_cos = cos(roll);
  133. C[0] = cos(pitch);
  134. C[1] = pitch_sin * roll_sin;
  135. C[2] = pitch_sin * roll_cos;
  136. C[4] = roll_cos;
  137. C[5] = -roll_sin;
  138. C[6] = -pitch_sin;
  139. C[7] = pitch_cos * roll_sin;
  140. C[8] = pitch_cos * roll_cos;
  141. }
  142. void reset_yaw_C(float *C)
  143. {
  144. float pitch = asin(-C[6]);
  145. float roll = atan2(C[7], C[8]);
  146. float pitch_sin = sin(pitch);
  147. float pitch_cos = cos(pitch);
  148. float roll_sin = sin(roll);
  149. float roll_cos = cos(roll);
  150. C[0] = pitch_cos;
  151. C[1] = pitch_sin * roll_sin;
  152. C[2] = pitch_sin * roll_cos;
  153. C[3] = 0.0;
  154. C[4] = roll_cos;
  155. C[5] = -roll_sin;
  156. }
  157. // F * P * F'
  158. void State_covariance_matrix_update(float *P, float *acc_n)
  159. {
  160. // P2 + P3*dt,
  161. P[3] = P[3] + P[6] * dt;
  162. P[4] = P[4] + P[7] * dt;
  163. P[5] = P[5] + P[8] * dt;
  164. P[12] = P[12] + P[15] * dt;
  165. P[13] = P[13] + P[16] * dt;
  166. P[14] = P[14] + P[17] * dt;
  167. P[21] = P[21] + P[24] * dt;
  168. P[22] = P[22] + P[25] * dt;
  169. P[23] = P[23] + P[26] * dt;
  170. //P4 + P7*dt,
  171. P[27] = P[27] + P[54] * dt;
  172. P[28] = P[28] + P[55] * dt;
  173. P[29] = P[29] + P[56] * dt;
  174. P[36] = P[36] + P[63] * dt;
  175. P[37] = P[37] + P[64] * dt;
  176. P[38] = P[38] + P[65] * dt;
  177. P[45] = P[45] + P[72] * dt;
  178. P[46] = P[46] + P[73] * dt;
  179. P[47] = P[47] + P[74] * dt;
  180. // P5 + P8*dt + dt*(P6 + P9*dt)
  181. P[30] = P[30] + P[57] * dt + dt * (P[33] + P[60] * dt);
  182. P[31] = P[31] + P[58] * dt + dt * (P[34] + P[61] * dt);
  183. P[32] = P[32] + P[59] * dt + dt * (P[35] + P[62] * dt);
  184. P[39] = P[39] + P[66] * dt + dt * (P[42] + P[69] * dt);
  185. P[40] = P[40] + P[67] * dt + dt * (P[43] + P[70] * dt);
  186. P[41] = P[41] + P[68] * dt + dt * (P[44] + P[71] * dt);
  187. P[48] = P[48] + P[75] * dt + dt * (P[51] + P[78] * dt);
  188. P[49] = P[49] + P[76] * dt + dt * (P[52] + P[79] * dt);
  189. P[50] = P[50] + P[77] * dt + dt * (P[53] + P[80] * dt);
  190. //P6 + P9*dt + matr*(P4 + P7*dt)
  191. P[33] = P[33] + P[60] * dt + acc_n[2] * dt * P[28] - acc_n[1] * dt * P[29];
  192. P[34] = P[34] + P[61] * dt - acc_n[2] * dt * P[27] + acc_n[0] * dt * P[29];
  193. P[35] = P[35] + P[62] * dt + acc_n[1] * dt * P[27] - acc_n[0] * dt * P[28];
  194. P[42] = P[42] + P[69] * dt + acc_n[2] * dt * P[37] - acc_n[1] * dt * P[38];
  195. P[43] = P[43] + P[70] * dt - acc_n[2] * dt * P[36] + acc_n[0] * dt * P[38];
  196. P[44] = P[44] + P[71] * dt + acc_n[1] * dt * P[36] - acc_n[0] * dt * P[37];
  197. P[51] = P[51] + P[78] * dt + acc_n[2] * dt * P[46] - acc_n[1] * dt * P[47];
  198. P[52] = P[52] + P[79] * dt - acc_n[2] * dt * P[45] + acc_n[0] * dt * P[47];
  199. P[53] = P[53] + P[80] * dt + acc_n[1] * dt * P[45] - acc_n[0] * dt * P[46];
  200. //P7 + P1*matr
  201. P[54] = P[54] + P[9] * acc_n[2] * dt - P[18] * acc_n[1] * dt;
  202. P[55] = P[55] + P[10] * acc_n[2] * dt - P[19] * acc_n[1] * dt;
  203. P[56] = P[56] + P[11] * acc_n[2] * dt - P[20] * acc_n[1] * dt;
  204. P[63] = P[63] - P[0] * acc_n[2] * dt + P[18] * acc_n[0] * dt;
  205. P[64] = P[64] - P[1] * acc_n[2] * dt + P[19] * acc_n[0] * dt;
  206. P[65] = P[65] - P[2] * acc_n[2] * dt + P[20] * acc_n[0] * dt;
  207. P[72] = P[72] + P[0] * acc_n[1] * dt - P[9] * acc_n[0] * dt;
  208. P[73] = P[73] + P[1] * acc_n[1] * dt - P[10] * acc_n[0] * dt;
  209. P[74] = P[74] + P[2] * acc_n[1] * dt - P[11] * acc_n[0] * dt;
  210. //P8 + P2*matr + dt*(P9 + P3*matr),
  211. P[57] = P[57] + dt * P[60] + P[12] * acc_n[2] * dt - P[21] * acc_n[1] * dt;
  212. P[58] = P[58] + dt * P[61] + P[13] * acc_n[2] * dt - P[22] * acc_n[1] * dt;
  213. P[59] = P[59] + dt * P[62] + P[14] * acc_n[2] * dt - P[23] * acc_n[1] * dt;
  214. P[66] = P[66] + dt * P[69] - P[3] * acc_n[2] * dt + P[21] * acc_n[0] * dt;
  215. P[67] = P[67] + dt * P[70] - P[4] * acc_n[2] * dt + P[22] * acc_n[0] * dt;
  216. P[68] = P[68] + dt * P[71] - P[5] * acc_n[2] * dt + P[23] * acc_n[0] * dt;
  217. P[75] = P[75] + dt * P[78] + P[3] * acc_n[1] * dt - P[12] * acc_n[0] * dt;
  218. P[76] = P[76] + dt * P[79] + P[4] * acc_n[1] * dt - P[13] * acc_n[0] * dt;
  219. P[77] = P[77] + dt * P[80] + P[5] * acc_n[1] * dt - P[14] * acc_n[0] * dt;
  220. // P9 + P3*matr + matr*(P7 + P1*matr)
  221. P[60] = P[60] + P[15] * acc_n[2] * dt - P[24] * acc_n[1] * dt + acc_n[2] * dt * P[55] - acc_n[1] * dt * P[56];
  222. P[61] = P[61] + P[16] * acc_n[2] * dt - P[25] * acc_n[1] * dt - acc_n[2] * dt * P[54] + acc_n[0] * dt * P[56];
  223. P[62] = P[62] + P[17] * acc_n[2] * dt - P[26] * acc_n[1] * dt + acc_n[1] * dt * P[54] - acc_n[0] * dt * P[55];
  224. P[69] = P[69] - P[6] * acc_n[2] * dt + P[24] * acc_n[0] * dt + acc_n[2] * dt * P[64] - acc_n[1] * dt * P[65];
  225. P[70] = P[70] - P[7] * acc_n[2] * dt + P[25] * acc_n[0] * dt - acc_n[2] * dt * P[63] + acc_n[0] * dt * P[65];
  226. P[71] = P[71] - P[8] * acc_n[2] * dt + P[26] * acc_n[0] * dt + acc_n[1] * dt * P[63] - acc_n[0] * dt * P[64];
  227. P[78] = P[78] + P[6] * acc_n[1] * dt - P[15] * acc_n[0] * dt + acc_n[2] * dt * P[73] - acc_n[1] * dt * P[74];
  228. P[79] = P[79] + P[7] * acc_n[1] * dt - P[16] * acc_n[0] * dt - acc_n[2] * dt * P[72] + acc_n[0] * dt * P[74];
  229. P[80] = P[80] + P[8] * acc_n[1] * dt - P[17] * acc_n[0] * dt + acc_n[1] * dt * P[72] - acc_n[0] * dt * P[73];
  230. //P3 + P1 * matr
  231. P[6] = P[6] + P[1] * acc_n[2] * dt - P[2] * acc_n[1] * dt;
  232. P[7] = P[7] - P[0] * acc_n[2] * dt + P[2] * acc_n[0] * dt;
  233. P[8] = P[8] + P[0] * acc_n[1] * dt - P[1] * acc_n[0] * dt;
  234. P[15] = P[15] + P[10] * acc_n[2] * dt - P[11] * acc_n[1] * dt;
  235. P[16] = P[16] - P[9] * acc_n[2] * dt + P[11] * acc_n[0] * dt;
  236. P[17] = P[17] + P[9] * acc_n[1] * dt - P[10] * acc_n[0] * dt;
  237. P[24] = P[24] + P[19] * acc_n[2] * dt - P[20] * acc_n[1] * dt;
  238. P[25] = P[25] - P[18] * acc_n[2] * dt + P[20] * acc_n[0] * dt;
  239. P[26] = P[26] + P[18] * acc_n[1] * dt - P[19] * acc_n[0] * dt;
  240. float noise = SIGMA * SIGMA * dt *dt;
  241. for (int i = 0; i < 9; i++)
  242. {
  243. P[i * 9 + i] += noise;
  244. }
  245. }
  246. void Kalfman_gain(float *P, float *Temporary_array, float *Temporary_array1, float *K)
  247. {
  248. Temporary_array[0] = P[60] + SIGMA_V * SIGMA_V * 0.01f;
  249. Temporary_array[1] = P[61];
  250. Temporary_array[2] = P[62];
  251. Temporary_array[3] = P[69];
  252. Temporary_array[4] = P[70] + SIGMA_V * SIGMA_V * 0.01f;
  253. Temporary_array[5] = P[71];
  254. Temporary_array[6] = P[78];
  255. Temporary_array[7] = P[79];
  256. Temporary_array[8] = P[80] + SIGMA_V * SIGMA_V * 0.01f;
  257. invert3x3(Temporary_array, Temporary_array1);
  258. memcpy(Temporary_array, Temporary_array1, 9 * sizeof(float));
  259. K[0] = P[6] * Temporary_array[0] + P[7] * Temporary_array[3] + P[8] * Temporary_array[6];
  260. K[1] = P[6] * Temporary_array[1] + P[7] * Temporary_array[4] + P[8] * Temporary_array[7];
  261. K[2] = P[6] * Temporary_array[2] + P[7] * Temporary_array[5] + P[8] * Temporary_array[8];
  262. K[3] = P[15] * Temporary_array[0] + P[16] * Temporary_array[3] + P[17] * Temporary_array[6];
  263. K[4] = P[15] * Temporary_array[1] + P[16] * Temporary_array[4] + P[17] * Temporary_array[7];
  264. K[5] = P[15] * Temporary_array[2] + P[16] * Temporary_array[5] + P[17] * Temporary_array[8];
  265. K[6] = P[24] * Temporary_array[0] + P[25] * Temporary_array[3] + P[26] * Temporary_array[6];
  266. K[7] = P[24] * Temporary_array[1] + P[25] * Temporary_array[4] + P[26] * Temporary_array[7];
  267. K[8] = P[24] * Temporary_array[2] + P[25] * Temporary_array[5] + P[26] * Temporary_array[8];
  268. K[9] = P[33] * Temporary_array[0] + P[34] * Temporary_array[3] + P[35] * Temporary_array[6];
  269. K[10] = P[33] * Temporary_array[1] + P[34] * Temporary_array[4] + P[35] * Temporary_array[7];
  270. K[11] = P[33] * Temporary_array[2] + P[34] * Temporary_array[5] + P[35] * Temporary_array[8];
  271. K[12] = P[42] * Temporary_array[0] + P[43] * Temporary_array[3] + P[44] * Temporary_array[6];
  272. K[13] = P[42] * Temporary_array[1] + P[43] * Temporary_array[4] + P[44] * Temporary_array[7];
  273. K[14] = P[42] * Temporary_array[2] + P[43] * Temporary_array[5] + P[44] * Temporary_array[8];
  274. K[15] = P[51] * Temporary_array[0] + P[52] * Temporary_array[3] + P[53] * Temporary_array[6];
  275. K[16] = P[51] * Temporary_array[1] + P[52] * Temporary_array[4] + P[53] * Temporary_array[7];
  276. K[17] = P[51] * Temporary_array[2] + P[52] * Temporary_array[5] + P[53] * Temporary_array[8];
  277. K[18] = P[60] * Temporary_array[0] + P[61] * Temporary_array[3] + P[62] * Temporary_array[6];
  278. K[19] = P[60] * Temporary_array[1] + P[61] * Temporary_array[4] + P[62] * Temporary_array[7];
  279. K[20] = P[60] * Temporary_array[2] + P[61] * Temporary_array[5] + P[62] * Temporary_array[8];
  280. K[21] = P[69] * Temporary_array[0] + P[70] * Temporary_array[3] + P[71] * Temporary_array[6];
  281. K[22] = P[69] * Temporary_array[1] + P[70] * Temporary_array[4] + P[71] * Temporary_array[7];
  282. K[23] = P[69] * Temporary_array[2] + P[70] * Temporary_array[5] + P[71] * Temporary_array[8];
  283. K[24] = P[78] * Temporary_array[0] + P[79] * Temporary_array[3] + P[80] * Temporary_array[6];
  284. K[25] = P[78] * Temporary_array[1] + P[79] * Temporary_array[4] + P[80] * Temporary_array[7];
  285. K[26] = P[78] * Temporary_array[2] + P[79] * Temporary_array[5] + P[80] * Temporary_array[8];
  286. }
  287. void multiply9x3(float *K, float *vel_n, float* delta_x)
  288. {
  289. int i = 0;
  290. for (i = 0; i < 9; i++)
  291. {
  292. delta_x[i] = K[i * 3] * vel_n[0] + K[i * 3 + 1] * vel_n[1] + K[i * 3 + 2] * vel_n[2];
  293. }
  294. }
  295. void State_covariance_matrix_corr(float *P, float *P_tmp, float *K)
  296. {
  297. int i = 0;
  298. int j = 0;
  299. for (i = 0; i < 9; i++) {
  300. for (j = 0; j < 9; j++) {
  301. P_tmp[i * 9 + j] = P[i * 9 + j] - K[3 * i] * P[54 + j] - K[3 * i + 1] * P[63 + j] - K[3 * i + 2] * P[72 + j];
  302. }
  303. }
  304. memcpy(P, P_tmp, 81 * sizeof(float));
  305. }
  306. void Att_matrix_corr(float *C, float *C_prev, float *Temporary_array, float *Temporary_array1, float *delta_x)
  307. {
  308. Temporary_array[0] = 2.0;
  309. Temporary_array[1] = -delta_x[2];
  310. Temporary_array[2] = delta_x[1];
  311. Temporary_array[3] = delta_x[2];
  312. Temporary_array[4] = 2.0;
  313. Temporary_array[5] = -delta_x[0];
  314. Temporary_array[6] = -delta_x[1];
  315. Temporary_array[7] = delta_x[0];
  316. Temporary_array[8] = 2.0;
  317. invert3x3(Temporary_array, Temporary_array1);
  318. Temporary_array[0] = 2.0;
  319. Temporary_array[1] = delta_x[2];
  320. Temporary_array[2] = -delta_x[1];
  321. Temporary_array[3] = -delta_x[2];
  322. Temporary_array[4] = 2.0;
  323. Temporary_array[5] = delta_x[0];
  324. Temporary_array[6] = delta_x[1];
  325. Temporary_array[7] = -delta_x[0];
  326. Temporary_array[8] = 2.0;
  327. multiply3x3(Temporary_array, Temporary_array1, C_prev);
  328. memcpy(Temporary_array, C_prev, 9 * sizeof(float));
  329. multiply3x3(Temporary_array, C, C_prev);
  330. memcpy(C, C_prev, 9 * sizeof(float));
  331. }
  332. void pos_n_corr(float *pos_n, float *delta_x)
  333. {
  334. pos_n[0] -= delta_x[3];
  335. pos_n[1] -= delta_x[4];
  336. pos_n[2] -= delta_x[5];
  337. }
  338. void vel_n_corr(float *vel_n, float *delta_x)
  339. {
  340. vel_n[0] -= delta_x[6];
  341. vel_n[1] -= delta_x[7];
  342. vel_n[2] -= delta_x[8];
  343. }
  344. void State_covariance_matrix_orthogonalization(float *P)
  345. {
  346. int i = 0;
  347. int j = 0;
  348. float temp;
  349. for (i = 0; i < 9; i++)
  350. for (j = i + 1; j < 9; j++)
  351. {
  352. temp = 0.5f*(P[i * 9 + j] + P[j * 9 + i]);
  353. P[i * 9 + j] = temp;
  354. P[j * 9 + i] = temp;
  355. }
  356. }
  357. void Initialize(float *gyr, float *acc)
  358. {
  359. frame_index = 1;
  360. stand_num = 0;
  361. accSize = 1.0f;
  362. accSum = 0.0f;
  363. ZUPT_STATUS = 0;
  364. HAS_RESULT = 0;
  365. memset(last_pos_n, 0, 3 * sizeof(float));
  366. memset(pos_offset, 0, 3 * sizeof(float));
  367. memset(gyr_norm_window, 0, 10 * sizeof(float));
  368. memset(gyr_z_window, 0, 10 * sizeof(float));
  369. memset(P, 0, 81 * sizeof(float));
  370. memset(acc_n, 0, 3 * sizeof(float));
  371. memset(vel_n, 0, 3 * sizeof(float));
  372. memset(pos_n, 0, 3 * sizeof(float));
  373. memset(Temporary_array1, 0, 9 * sizeof(float));
  374. memset(Temporary_array2, 0, 9 * sizeof(float));
  375. memset(K, 0, 27 * sizeof(float));
  376. memset(P_prev, 0, 81 * sizeof(float));
  377. memset(delta_x, 0, 9 * sizeof(float));
  378. memset(C, 0, 9 * sizeof(float));
  379. memset(Temporary_array1, 0, 9 * sizeof(float));
  380. memset(Temporary_array2, 0, 9 * sizeof(float));
  381. memset(press_data, 0, 10 * sizeof(int));
  382. init_attitude_matrix(C, acc);
  383. memcpy(C_prev, C, 9 * sizeof(float));
  384. }
  385. void attitude_matrix_update(float *C, float *Temporary_array1, float *Temporary_array2, float *gyr, float dt)
  386. {
  387. Temporary_array1[0] = 2.0f;
  388. Temporary_array1[1] = dt * gyr[2];
  389. Temporary_array1[2] = -dt * gyr[1];
  390. Temporary_array1[3] = -dt * gyr[2];
  391. Temporary_array1[4] = 2.0f;
  392. Temporary_array1[5] = dt * gyr[0];
  393. Temporary_array1[6] = dt * gyr[1];
  394. Temporary_array1[7] = -dt * gyr[0];
  395. Temporary_array1[8] = 2.0f;
  396. invert3x3(Temporary_array1, Temporary_array2);
  397. memset(Temporary_array1, 0, 9 * sizeof(float));
  398. Temporary_array1[0] = 2 * C[0] + C[1] * dt * gyr[2] - C[2] * dt * gyr[1];
  399. Temporary_array1[1] = 2 * C[1] - C[0] * dt * gyr[2] + C[2] * dt * gyr[0];
  400. Temporary_array1[2] = 2 * C[2] + C[0] * dt * gyr[1] - C[1] * dt * gyr[0];
  401. Temporary_array1[3] = 2 * C[3] + C[4] * dt * gyr[2] - C[5] * dt * gyr[1];
  402. Temporary_array1[4] = 2 * C[4] - C[3] * dt * gyr[2] + C[5] * dt * gyr[0];
  403. Temporary_array1[5] = 2 * C[5] + C[3] * dt * gyr[1] - C[4] * dt * gyr[0];
  404. Temporary_array1[6] = 2 * C[6] + C[7] * dt * gyr[2] - C[8] * dt * gyr[1];
  405. Temporary_array1[7] = 2 * C[7] - C[6] * dt * gyr[2] + C[8] * dt * gyr[0];
  406. Temporary_array1[8] = 2 * C[8] + C[6] * dt * gyr[1] - C[7] * dt * gyr[0];
  407. multiply3x3(Temporary_array1, Temporary_array2, C);
  408. }
  409. float max_window_val(float *window, int window_size)
  410. {
  411. float val = window[0];
  412. for (int i = 0; i < window_size; i++)
  413. {
  414. if (window[i] > val)
  415. val = window[i];
  416. }
  417. return val;
  418. }
  419. int max_window_int(int *window, int window_size)
  420. {
  421. int val = window[0];
  422. for (int i = 0; i < window_size; i++)
  423. {
  424. if (window[i] > val)
  425. val = window[i];
  426. }
  427. return val;
  428. }
  429. float min_window_val(float *window, int window_size)
  430. {
  431. float val = window[0];
  432. for (int i = 0; i < window_size; i++)
  433. {
  434. if (window[i] < val)
  435. val = window[i];
  436. }
  437. return val;
  438. }
  439. int min_window_int(int *window, int window_size)
  440. {
  441. int val = window[0];
  442. for (int i = 0; i < window_size; i++)
  443. {
  444. if (window[i] < val)
  445. val = window[i];
  446. }
  447. return val;
  448. }
  449. //press_tren 函数功能:提供走路过程中上升沿,下降沿
  450. //1 为上升 2 为 下降 0为不需要得状态
  451. int press_trend(int index ,int *window, int window_size)
  452. {
  453. int i;
  454. int max_val = window[(index - 1) % window_size];
  455. int max_index = index;
  456. int min_val = max_val;
  457. int min_index = max_index;
  458. for(i = 1; i < window_size + 1; i++)
  459. {
  460. if(max_val < window[(index - i) % window_size])
  461. {
  462. max_index = index - i + 1;
  463. max_val = window[(index - i) % window_size];
  464. }
  465. if(min_val > window[(index - i) % window_size])
  466. {
  467. min_index = index - i + 1;
  468. min_val = window[(index - i) % window_size];
  469. }
  470. }
  471. if(max_index > min_index && max_val > min_val + 50000)
  472. {
  473. return 1;
  474. }
  475. if(max_index < min_index && max_val > min_val + 50000)
  476. {
  477. return 2;
  478. }
  479. return 0;
  480. }
  481. unsigned char footPDR(int num, float *gyr, float *acc, int press, short* pos_res, short* angle_data)
  482. {
  483. unsigned char movement_e = 0;
  484. for (int i = 0; i < 3; i++)
  485. {
  486. gyr[i] *= (PI / 180);
  487. acc[i] *= g;
  488. }
  489. if(num_peak == 0)
  490. {
  491. for(int i = 0; i < 3; i++)
  492. {
  493. gyr_extreme[2 * i] = gyr[i];
  494. gyr_extreme[2 * i + 1] = gyr[i];
  495. }
  496. }
  497. for (int i = 0; i < 3; i++)
  498. {
  499. if (gyr[i] < gyr_extreme[2 * i])
  500. {
  501. gyr_extreme[2 * i] = gyr[i];
  502. }
  503. if (gyr[i] > gyr_extreme[2 * i + 1])
  504. {
  505. gyr_extreme[2 * i + 1] = gyr[i];
  506. }
  507. }
  508. accSum += sqrt(acc[0] * acc[0] + acc[1] * acc[1] + acc[2] * acc[2]);
  509. for (int i = 0; i < 3; i++)
  510. {
  511. gyr_mean[i] += gyr[i];
  512. }
  513. num_peak++;
  514. if (num_peak > 500)
  515. {
  516. if (gyr_extreme[1] - gyr_extreme[0] < 0.005f && gyr_extreme[3] - gyr_extreme[2] < 0.005f && gyr_extreme[5] - gyr_extreme[4] < 0.005f)
  517. {
  518. gyrBias[0] = gyr_mean[0] / num_peak;
  519. gyrBias[1] = gyr_mean[1] / num_peak;
  520. gyrBias[2] = gyr_mean[2] / num_peak;
  521. accSize = g * num_peak /accSum;
  522. }
  523. num_peak = 0;
  524. accSum = 0.0f;
  525. memset(gyr_mean, 0, 3 * sizeof(float));
  526. }
  527. gyr[0] -= gyrBias[0];
  528. gyr[1] -= gyrBias[1];
  529. gyr[2] -= (gyrBias[2]);
  530. acc[0] *= accSize;
  531. acc[1] *= accSize;
  532. acc[2] *= accSize;
  533. float gyr_norm_xyz = sqrt(gyr[0]*gyr[0] + gyr[1]*gyr[1] + gyr[2]*gyr[2]);
  534. //需要一个滑动窗口来判断脚步是否在地上
  535. frame_index++;
  536. //下面为惯导解算
  537. if (num == 1)
  538. {
  539. Initialize(gyr, acc);
  540. return movement_e;
  541. }
  542. //惯导解算拆分为5次迭代,缓解高量程下由于采样率过导致惯导解算有错
  543. //插值法
  544. float gyr_temp[3];
  545. float acc_temp[3];
  546. gyr_temp[0] = gyr[0] - last_gyr[0];
  547. gyr_temp[1] = gyr[1] - last_gyr[1];
  548. gyr_temp[2] = gyr[2] - last_gyr[2];
  549. acc_temp[0] = acc[0] - last_acc[0];
  550. acc_temp[1] = acc[1] - last_acc[1];
  551. acc_temp[2] = acc[2] - last_acc[2];
  552. for(int i = 1; i < 6; i++)
  553. {
  554. last_gyr[0] += 0.2f * gyr_temp[0];
  555. last_gyr[1] += 0.2f * gyr_temp[1];
  556. last_gyr[2] += 0.2f * gyr_temp[2];
  557. last_acc[0] += 0.2f * acc_temp[0];
  558. last_acc[1] += 0.2f * acc_temp[1];
  559. last_acc[2] += 0.2f * acc_temp[2];
  560. attitude_matrix_update(C, Temporary_array1, Temporary_array2, last_gyr, 0.2f*dt);
  561. multiply3x1(C, last_acc, acc_n);
  562. vel_n[0] = last_vel_n[0] + acc_n[0] * dt * 0.2f;
  563. vel_n[1] = last_vel_n[1] + acc_n[1] * dt * 0.2f;
  564. vel_n[2] = last_vel_n[2] + (acc_n[2] - g) * dt * 0.2f;
  565. pos_n[0] = pos_n[0] + (vel_n[0] + last_vel_n[0]) * dt * 0.1f;
  566. pos_n[1] = pos_n[1] + (vel_n[1] + last_vel_n[1]) * dt * 0.1f;
  567. pos_n[2] = pos_n[2] + (vel_n[2] + last_vel_n[2]) * dt * 0.1f;
  568. memcpy(last_vel_n, vel_n, 3 * sizeof(float));
  569. }
  570. // attitude_matrix_update(C, Temporary_array1, Temporary_array2, gyr, dt);
  571. // multiply3x1(C, acc, acc_n);
  572. // vel_n[0] = vel_n[0] + acc_n[0] * dt;
  573. // vel_n[1] = vel_n[1] + acc_n[1] * dt;
  574. // vel_n[2] = vel_n[2] + (acc_n[2] - g) * dt;
  575. // pos_n[0] = pos_n[0] + vel_n[0] * dt;
  576. // pos_n[1] = pos_n[1] + vel_n[1] * dt;
  577. // pos_n[2] = pos_n[2] + vel_n[2] * dt;
  578. memcpy(last_gyr, gyr, 3 * sizeof(float));
  579. memcpy(last_acc, acc, 3 * sizeof(float));
  580. //P = F*P*F' + Q;
  581. State_covariance_matrix_update(P, acc_n);
  582. int window_index = (frame_index - 1) % 10;
  583. float gyr_norm_xz = sqrt(gyr[0] * gyr[0] + gyr[1] * gyr[1] +gyr[2] * gyr[2]);
  584. gyr_norm_window[window_index] = gyr_norm_xz;
  585. press_data[window_index] = press;
  586. //当press_trend函数返回是1,判断为踩地上
  587. // 返回2 的时候,判断为离地
  588. // 返回0 的时候,需要保持状态
  589. int press_trend_val = press_trend(frame_index, press_data, 10);
  590. if(press_trend_val == 1)
  591. {
  592. ZUPT_STATUS = 1;
  593. }
  594. else if(press_trend_val == 2)
  595. {
  596. ZUPT_STATUS = 2;
  597. }
  598. //RUN_ZUPT mean detect on floor when running
  599. int RUN_ZUPT = 0;
  600. if((frame_index > 10 && ZUPT_STATUS == 1))
  601. RUN_ZUPT = 1;
  602. //STAND_ZUPT mean detect on floor when no any moving
  603. int STAND_ZUPT = 0;
  604. if((frame_index > 15 && gyr_norm_window[window_index] < 0.35f && fabs(min_window_val(gyr_norm_window, 10) - max_window_val(gyr_norm_window, 10)) < 0.1f))
  605. STAND_ZUPT = 1;
  606. //zupt
  607. if (STAND_ZUPT ||(RUN_ZUPT))
  608. {
  609. //计算一步的距离及步数+1
  610. if(last_stage == 0)
  611. {
  612. cal_step_data();
  613. }
  614. stand_num = stand_num + 1;
  615. //K = P*H'/(H*P*H' + R);
  616. Kalfman_gain(P, Temporary_array1, Temporary_array2, K);
  617. //delta_x = K * [vel_n(:,i);];
  618. multiply9x3(K, vel_n, delta_x);
  619. State_covariance_matrix_corr(P, P_prev, K);
  620. //这里先从设置 delta_x(3) = atan2(C(2,1),C(1,1));
  621. //意味着每一步的参考方向是IMU X轴方向
  622. delta_x[2] = atan2(C[3], C[0]);
  623. // theta = -1.7801 + atan2(C[3], C[0]);
  624. theta = 0.0f;
  625. Att_matrix_corr(C, C_prev, Temporary_array1, Temporary_array2, delta_x);
  626. pos_n_corr(pos_n, delta_x);
  627. vel_n_corr(vel_n, delta_x);
  628. memset(vel_n, 0, 3*sizeof(float));
  629. last_pos_n[0] = pos_n[0];
  630. last_pos_n[1] = pos_n[1];
  631. last_pos_n[2] = pos_n[2];
  632. HAS_RESULT = 1;
  633. last_stage = 1;
  634. }
  635. else
  636. {
  637. stand_num = 0;
  638. last_stage = 0;
  639. }
  640. State_covariance_matrix_orthogonalization(P);
  641. memcpy(last_vel_n, vel_n, 3 * sizeof(float));
  642. /*theta = -0.61;
  643. temp = [cos(theta), sin(theta); -sin(theta), cos(theta)];
  644. pos_temp = temp * [pos_result(1); pos_result(2)];
  645. pos_result(1) = pos_temp(1);
  646. pos_result(2) = pos_temp(2);*/
  647. float pos_offset_temp0 = pos_n[0] - last_pos_n[0];
  648. float pos_offset_temp1 = pos_n[1] - last_pos_n[1];
  649. float pos_offset_temp2 = pos_n[2] - last_pos_n[2];
  650. // pos_offset[0] = cos(theta) * pos_offset_temp0 + sin(theta) * pos_offset_temp1;
  651. // pos_offset[1] = -sin(theta) * pos_offset_temp0+ cos(theta) * pos_offset_temp1;
  652. pos_offset[0] = pos_offset_temp0;
  653. pos_offset[1] = pos_offset_temp1;
  654. pos_offset[2] = pos_offset_temp2;
  655. // memcpy(pos_res, acc_n, 3 * sizeof(float));
  656. if(HAS_RESULT == 1 && pos_offset[0] > 0.2f)
  657. {
  658. movement_e = 1;
  659. HAS_RESULT = 0;
  660. }
  661. // angle_data[0] = (short) (asin(-C[6]) * 10000.f); //pitch
  662. // angle_data[1] = (short) (atan2(C[7], C[8]) * 10000.f); //roll
  663. // angle_data[2] = (short) (atan2(C[3], C[0]) * 10000.f); //yaw
  664. pos_res[0] = (short) (pos_offset[0] * 100.0f);
  665. pos_res[1] = (short) (pos_offset[1] * 100.0f);
  666. pos_res[2] = (short) (pos_offset[2] * 100.0f);
  667. return movement_e;
  668. }