以下的测试程序均基于windows平台用visual studio 2019调试。
拉格朗日拟合测试文件
程序清单1拉格朗日二次分段插值法拟合测试Main.h
________________________________________
1. #ifndef MAIN_H
2. #define MAIN_H
3.
4. /*该工程作为单片机分段二次拉格朗日插值法的测试工程*/
5.
6. //变量声明放在c文件，函数声明放在h文件
7.
8. //为了配合单片机，变量尽量用float
9.
10. //参考资料：计算机数值方法 65-66 72-73
11.
12. /*********************************************************************************************************
13. * 模块名称: 拉格朗日二次插值测试工程
14. * 摘 要:
15. * 当前版本: 1.00
16. * 作 者: guangjie2333
17. * 完成日期: 2021-4-15
18. *
19. **********************************************************************************************************
20.
21.
________________________________________
程序清单2拉格朗日二次分段插值法拟合测试ALG.h
________________________________________
1. #ifndef ALG_H
2. #define ALG_H
3.
4. /*该文件是算法函数的头文件*/
5.
6. #include "stdio.h"
7. #include "math.h"
8. #include "stdlib.h"
9.
10. /*tips :
11. malloc 之后记得free
12. 传入*p改变值，传入**p改变地址
13. 单片机堆区512字节 float 128个
14. */
15.
16. //这个工程多用一下指针，看看自己水平如何
17.
18. /************************************
19. * 宏定义
20. ************************************/
21.
22. /************************************
23. * 结构体声明
24. ************************************/
25.
26. typedef struct {
27. float value;
28. unsigned short index;
29. }Node;
30.
31. //0到250的拉格朗日分段二次拟合
32. void LagrangeFit(float* cuff, unsigned short* paluse, int s_iPluseStepCnt, float* y);
33.
34. //数组逆排续
35. void u16_reserverl(unsigned short a[], int n);
36. void float_reserverl(float a[], int n);
37.
38. //找最大值和最大值所在的数组
39. Node FindMax(float* y);
40. #endif
41.
________________________________________
程序清单3拉格朗日二次分段插值法拟合测试Main.c
________________________________________
1. #include "ALG.h"
2.
3. /*内部变量定义*/
4.
5. //拟合曲线的袖带压AD
6. static float cuffPressAD[14] = { 2863,2717,2586,2461,2346,2220,2117,
7. 2053,1950,1882,1807,1718,1643,1557 };
8. //拟合曲线的脉搏波AD
9. static unsigned short paluseAD[14] = { 5622,4007,3004,4229,8464,9874,12125,
10. 10680,8704,6426,4256,4114,3555,3243 };
11. //放气阶数
12. int s_iPluseStepCnt = 14;
13.
14. void main()
15. {
16. float y[250] = {0};
17. Node maxNode;
18. float dap, sap, map;
19.
20. for (int i = 0; i < s_iPluseStepCnt; i++)
21. {
22. cuffPressAD[i] = (float)(cuffPressAD[i] * 0.097 - 108); //AD-真实值 线性转换公式
23. paluseAD[i] = (unsigned short)(paluseAD[i] / 16); //脉搏波AD放大了16倍
24. }
25. u16_reserverl(paluseAD,14);
26. float_reserverl(cuffPressAD,14);
27. //拉格朗日分段二次插值
28. LagrangeFit(cuffPressAD,paluseAD,s_iPluseStepCnt,y);
29.
30. for (int i = 0; i < 250; i++)
31. {
32. printf("%f \n",y[i]);
33. }
34.
35. maxNode = FindMax(y);
36.
37. //找舒张压
38. for (int i = maxNode.index; i > 0; i--)
39. {
40. if (y[i] <= maxNode.value * 0.8)
41. {
42. printf("dap = %d\n",i);
43. dap = i;
44. break;
45. }
46. }
47.
48. //找收缩压
49. for (int i = maxNode.index; i < 250; i++)
50. {
51. if (y[i] <= maxNode.value * 0.6)
52. {
53. printf("sap = %d\n", i);
54. sap = i;
55. break;
56. }
57. }
58.
59. printf("map = %f \n",dap*2/3+sap/3);
60.
61. }
________________________________________
程序清单4拉格朗日二次分段插值法拟合测试ALG.c
________________________________________
1. #include "ALG.h"
2.
3. /*********************************************************************************************************
4. * 内部函数定义
5. *********************************************************************************************************/
6.
7. //返回拉格朗日插值的l0
8. static float l_0(Node node0, Node node1, Node node2,float x)
9. {
10. float r;
11.
12. r = (x - node1.index) * (x - node2.index) / (node0.index - node1.index) / (node0.index - node2.index);
13.
14. return r;
15. }
16.
17. //返回拉格朗日插值的l1
18. static float l_1(Node node0, Node node1, Node node2, float x)
19. {
20. float r;
21.
22. r = (x - node0.index) * (x - node2.index) / (node1.index - node0.index) / (node1.index - node2.index);
23.
24. return r;
25. }
26.
27. //返回拉格朗日插值的l2
28. static float l_2(Node node0, Node node1, Node node2, float x)
29. {
30. float r;
31.
32. r = (x - node0.index) * (x - node1.index) / (node2.index - node0.index) / (node2.index - node1.index);
33.
34. return r;
35.
36. }
37.
38. //返回分段拉格朗日插值的数组 value-脉搏波 index-袖带压
39. static float L2(Node node0, Node node1, Node node2,int x)
40. {
41. return (float)( node0.value * l_0(node0,node1,node2,x) + node1.value * l_1(node0, node1, node2, x) + node2.value * l_2(node0, node1, node2, x));
42. }
43.
44. /*********************************************************************************************************
45. * 外部函数定义
46. *********************************************************************************************************/
47. void LagrangeFit(float* cuff,unsigned short * paluse,int s_iPluseStepCnt,float* y)
48. {
49. int i,j;
50. int diff_Left, diff_Right;
51. //二次多项式拟合需要三个节点
52. Node node0;
53. Node node1;
54. Node node2;
55.
56. for (i = 0;i < 250; i++)
57. {
58. if (i<cuff[0] || i >= cuff[s_iPluseStepCnt-1]) //在采样数据之外的数据先赋值为0，也可以用直线连一下
59. {
60.
61. y[i] = 0;
62. if (i >= cuff[0] - 1)
63. {
64. printf("你过来啊 ");
65. }
66. }
67. else
68. {
69. printf("你过来啊 ");
70. if (i>=cuff[0] && i<cuff[1])
71. {
72. node0.index = cuff[0];
73. node0.value = paluse[0];
74. node1.index = cuff[1];
75. node1.value = paluse[1];
76. node2.index = cuff[2];
77. node2.value = paluse[2];
78. y[i] = L2(node0,node1,node2,i);
79. }
80. else if (i>=cuff[s_iPluseStepCnt-2] && i <cuff[s_iPluseStepCnt-1])
81. {
82. node0.index = cuff[s_iPluseStepCnt-3];
83. node0.value = paluse[s_iPluseStepCnt-3];
84. node1.index = cuff[s_iPluseStepCnt-2];
85. node1.value = paluse[s_iPluseStepCnt-2];
86. node2.index = cuff[s_iPluseStepCnt-1];
87. node2.value = paluse[s_iPluseStepCnt-1];
88. y[i] = L2(node0, node1, node2, i);
89. }
90. else
91. {
92. for ( j = 0; j < s_iPluseStepCnt; j++)
93. {
94. if (i >= cuff[j] && i < cuff[j+1])
95. {
96. diff_Left = abs(i - cuff[j-1]);
97. diff_Right = abs(i - cuff[j+2]);
98. if (diff_Left > diff_Right) //这个点离左边比较远
99. {
100. node0.index = cuff[j];
101. node0.value = paluse[j];
102. node1.index = cuff[j+1];
103. node1.value = paluse[j+1];
104. node2.index = cuff[j+2];
105. node2.value = paluse[j+2];
106. y[i] = L2(node0, node1, node2, i);
107. }
108. else
109. {
110. node0.index = cuff[j-1];
111. node0.value = paluse[j-1];
112. node1.index = cuff[j];
113. node1.value = paluse[j];
114. node2.index = cuff[j + 1];
115. node2.value = paluse[j + 1];
116. y[i] = L2(node0, node1, node2, i);
117. }
118. }
119. }
120. }
121. }
122. }
123. }
124.
125. //数组逆排
126. void float_reserverl(float a[], int n)
127. {
128. int m = (n + 1) / 2;
129. int i;
130. int j;
131. float temp;
132. for (i = 0; i < m; i++)
133. {
134. j = n - 1 - i;
135. temp = a[i];
136. a[i] = a[j];
137. a[j] = temp;
138. }
139. }
140.
141. //数组逆排
142. void u16_reserverl(unsigned short a[], int n)
143. {
144. int m = (n + 1) / 2;
145. int i;
146. int j;
147. unsigned short temp;
148. for (i = 0; i < m; i++)
149. {
150. j = n - 1 - i;
151. temp = a[i];
152. a[i] = a[j];
153. a[j] = temp;
154. }
155. }
156.
157. //找最大值和最大值所在的数组
158. Node FindMax(float *y)
159. {
160. int i;
161. Node maxNode;
162. maxNode.value = y[0]; maxNode.index = 0;
163. for (i = 0; i< 250;i++)
164. {
165. if (y[i] > maxNode.value)
166. {
167. maxNode.value = y[i];
168. maxNode.index = i;
169. }
170. }
171. return maxNode;
172. }
173.
________________________________________
程序清单5最小二乘法高斯拟合Main.h
________________________________________
1. #ifndef GUASSFITTEST_H
2. #define GUASSFITTEST_H
3. #include "ALG.h"
4. /*该工程作为单片机最小二乘法高斯拟合的测试工程*/
5.
6. //变量声明放在c文件，函数声明放在h文件
7.
8. //高斯公式原理可参考:
9. //https://blog.csdn.net/dingzj2000/article/details/103719368
10. //数值计算方法教材 103-108页
11. // fai0 = 1 fai1 = x fai2 = x^2
12.
13. #endif
________________________________________
程序清单6最小二乘法高斯拟合ALG.h
________________________________________
1. #ifndef ALG_H
2. #define ALG_H
3.
4. /*该文件是算法函数的头文件*/
5.
6. #include "stdio.h"
7. #include "math.h"
8.
9.
10. /*************************************
11. * 宏定义区
12. **************************************/
13.
14. #define ROW 3 //行
15. #define LINE 4 //列
16.
17. typedef struct
18. {
19. float b0;
20. float b1;
21. float b2;
22. }B_parameter;
23.
24. typedef struct
25. {
26. float a;
27. float b;
28. float c;
29. }Gaussian_parameter;
30.
31. typedef struct
32. {
33. float value;
34. int index;
35. }matrixData;
36.
37. //高斯列主元素消元法，返回b0，b1, b2
38. B_parameter GaussianElimination(float* y, int* x, int s_iPluseStepCnt);
39. //计算高斯曲线参数
40. Gaussian_parameter CalcGaussianParameter(B_parameter B);
41. //0~250高斯拟合
42. void GuassFit(Gaussian_parameter gaussian);
43. #endif // !ALG
44.
________________________________________
程序清单7最小二乘法高斯拟合Main.c
________________________________________
1. \/**************************************
2. 2021 04 10
3. 作者：SZU guangjie2333
4. 功能：c语言实现最小二乘法高斯拟合
5. *************************************/
6.
7.
8. /*头文件引用*/
9. #include "GuassFitTest.h"
10.
11.
12. /*内部变量定义*/
13.
14. //拟合曲线的袖带压AD
15. static int cuffPressAD[14] = { 2863,2717,2586,2461,2346,2220,2117,
16. 2053,1950,1882,1807,1718,1643,1557 };
17. //拟合曲线的脉搏波AD
18. static int paluseAD[14] = { 5622,4007,3004,4229,8464,9874,12125,
19. 10680,8704,6426,4256,4114,3555,3243 };
20. //放气阶数
21. int s_iPluseStepCnt = 14;
22.
23.
24. /*主函数*/
25. void main()
26. {
27. B_parameter B;
28. Gaussian_parameter guassian;
29. float Z[14]; //高斯函数两边取对数，变形，得到一个二次函数。 Z[i] = ln y[i] = ln paluseAD[i]
30. for (int i = 0; i < s_iPluseStepCnt; i++)
31. {
32. cuffPressAD[i] =(int)(cuffPressAD[i] * 0.097 - 108); //AD-真实值 线性转换公式
33. paluseAD[i] =(int)(paluseAD[i] / 16); //脉搏波AD放大了16倍
34. Z[i] = log(paluseAD[i]);
35. printf("cuffPressAD[%d] = %d \n paluseAD[%d] = %d \n Z[%d] = %f \n",i, cuffPressAD[i], i,paluseAD[i],i,Z[i]);
36. }
37. //计算最小二乘法的B向量
38. B = GaussianElimination(Z, cuffPressAD,s_iPluseStepCnt);
39. //计算高斯曲线的参数
40. guassian = CalcGaussianParameter(B);
41. //0~250点进行高斯拟合
42. GuassFit(guassian);
43.
44. printf(" \n\n b0 = %f b1 = %f b2 = %f \n\n",B.b0,B.b1,B.b2);
45. printf(" \n\n a = %f b = %f c = %f \n\n",guassian.a, guassian.b, guassian.c);
46. }
47.
48.
________________________________________
程序清单8最小二乘法高斯拟合测试ALG.c
________________________________________
1. #include"ALG.h"
2.
3. /*该文件用于高斯拟合最小二乘法算法的c语言实现*/
4.
5.
6.
7. /*************************************
8. * 内部变量定义区
9. **************************************/
10.
11.
12. //权重（由于前面的点会由于刚刚关阀的影响，脉搏波波动大，故权重为0）
13. static int w[14] = {0,0,0,1,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 };
14.
15. /*************************************
16. * 内部函数定义区
17. **************************************/
18.
19. //(fai0,fai0)的内积，返回1 求和
20. static float xfun0(int s_iPluseStepCnt)
21. {
22. float sum = 0;
23. for (int i = 0; i < s_iPluseStepCnt; i++)
24. {
25. sum += 1 * w[i];
26. }
27.
28. return sum;
29. }
30.
31. //(fai0,fai1)的内积，返回x 求和
32. static float xfun1(int* x, int s_iPluseStepCnt)
33. {
34. float sum = 0;
35. for (int i = 0;i < s_iPluseStepCnt ; i++)
36. {
37. sum += x[i] * w[i];
38. }
39.
40. return sum;
41. }
42.
43. //(fai0,fai2) || (fai1,fai1)的内积，返回x^2 求和
44. static float xfun2(int* x, int s_iPluseStepCnt)
45. {
46. float sum = 0;
47. for (int i = 0; i < s_iPluseStepCnt; i++)
48. {
49. sum += x[i]*x[i] * w[i];
50. }
51.
52. return sum;
53. }
54.
55.
56. //(fai1,fai2)的内积，返回x^3 求和
57. static float xfun3(int* x, int s_iPluseStepCnt)
58. {
59. float sum = 0;
60. for (int i = 0; i < s_iPluseStepCnt; i++)
61. {
62. sum += x[i] * x[i] * x[i] * w[i];
63. }
64.
65. return sum;
66. }
67.
68. //(fai2,fai2)的内积，返回x^4 求和
69. static float xfun4(int* x, int s_iPluseStepCnt)
70. {
71. float sum = 0;
72. for (int i = 0; i < s_iPluseStepCnt; i++)
73. {
74. sum += x[i] * x[i] * x[i] * x[i] * w[i];
75. }
76.
77. return sum;
78. }
79.
80.
81. //(y,fai0)的内积，返回y*1 求和
82. static float yfun0(float* y, int* x ,int s_iPluseStepCnt)
83. {
84. float sum = 0;
85. for (int i = 0; i < s_iPluseStepCnt; i++)
86. {
87. sum += y[i] * 1 * w[i];
88. }
89.
90. return sum;
91. }
92.
93. //(y,fai1)的内积，返回y*x 求和
94. static float yfun1(float* y, int* x, int s_iPluseStepCnt)
95. {
96. float sum = 0;
97. for (int i = 0; i < s_iPluseStepCnt; i++)
98. {
99. sum += y[i] * x[i] * w[i];
100. }
101.
102. return sum;
103. }
104.
105. //(y,fai2)的内积，返回y*1 求和
106. static float yfun2(float* y, int* x, int s_iPluseStepCnt)
107. {
108. float sum = 0;
109. for (int i = 0; i < s_iPluseStepCnt; i++)
110. {
111. sum += y[i] * x[i] * x[i] * w[i];
112. }
113.
114. return sum;
115. }
116.
117. //选第k次交换的列最大值当主元，返回所在行
118. static matrixData row_k_max(float matrix[ROW][LINE],int k)
119. {
120. matrixData line[ROW] ;
121. matrixData Max;
122.
123. //赋初值
124. for (int i = k; i < ROW; i++)
125. {
126. line[i].value = matrix[i][k];
127. line[i].index = i;
128. }
129. //比较大小
130. Max.value = 0;
131. for (int i = k; i < ROW; i++)
132. {
133. if (fabs(line[i].value) > fabs(Max.value))
134. {
135. Max = line[i];
136. }
137. }
138. return Max;
139. }
140.
141. //交换两行
142. static void swapRow(float matrix[ROW][LINE],int k, int maxRowE)
143. {
144. float temp = 0;
145. for (int i = 0; i < LINE; i++)
146. {
147. temp = matrix[k][i];
148. matrix[k][i] = matrix[maxRowE][i];
149. matrix[maxRowE][i] = temp;
150. }
151. }
152.
153. static void printfMatrix(float matrix[ROW][LINE])
154. {
155. for (int i = 0; i < ROW; i++)
156. {
157. for (int j = 0 ; j < LINE; j++)
158. {
159. printf(" %f ", matrix[i][j]);
160. }
161. printf("\n");
162. }
163. }
164.
165.
166. /*
167. 外部函数调用
168. */
169.
170. //高斯列主元素消元法，返回b0，b1, b2
171. B_parameter GaussianElimination(float* y, int* x, int s_iPluseStepCnt)
172. {
173. B_parameter B; //最小二乘法的待求解向量
174. matrixData Max; //第k次交换，列的最大值和索引
175. float m;
176.
177. //计算增广矩阵中的每个值
178. float matrix[ROW][LINE];
179. for (int i = 0; i < ROW; i++)
180. {
181. for (int j = 0; j < LINE; j++)
182. {
183. if (LINE-1 == j)
184. {
185. //计算最后一列即Z向量的值
186. switch (i)
187. {
188. case 0: matrix[i][j] = yfun0(y, x, s_iPluseStepCnt); break;
189. case 1: matrix[i][j] = yfun1(y, x, s_iPluseStepCnt); break;
190. case 2: matrix[i][j] = yfun2(y, x, s_iPluseStepCnt); break;
191. default: break;
192. }
193. }
194. else
195. {
196. //计算x矩阵的每个值
197. switch (i+j)
198. {
199. case 0: matrix[i][j] = xfun0(s_iPluseStepCnt); break;
200. case 1: matrix[i][j] = xfun1(x, s_iPluseStepCnt); break;
201. case 2: matrix[i][j] = xfun2(x, s_iPluseStepCnt); break;
202. case 3: matrix[i][j] = xfun3(x, s_iPluseStepCnt); break;
203. case 4: matrix[i][j] = xfun4(x, s_iPluseStepCnt); break;
204. default: break;
205. }
206. }
207.
208. }
209. }
210.
211. printf("打印初始的矩阵\n\n");
212. printfMatrix(matrix);
213.
214. //进行列主元素消元
215. for (int i = 0; i < ROW-1; i++) //i代表第i次交换
216. {
217. //选主元
218. Max = row_k_max(matrix,i);
219. //交换行
220. if (Max.index != i)
221. {
222. swapRow(matrix,i,Max.index);
223. }
224. //消主元
225. for (int j = i+1; j < ROW; j++)
226. {
227. m = matrix[j][i] / matrix[i][i]; //小除大，避免小值当分母
228. for (int k = i; k < LINE ; k++)
229. {
230. matrix[j][k] = matrix[j][k] - matrix[i][k] * m; //消元
231. }
232. }
233.
234. printf("\n\n打印第%d次变化后的矩阵\n\n",i);
235. printfMatrix(matrix);
236.
237. }
238.
239. //计算最小二乘法的B列向量
240. B.b2 = matrix[ROW - 1][LINE - 1] / matrix[ROW - 1][LINE - 2];
241. B.b1 = (matrix[ROW - 2][LINE - 1] - B.b2 * matrix[ROW - 2][LINE - 2])/ matrix[ROW - 2][LINE - 3];
242. B.b0 = (matrix[ROW - 3][LINE - 1] - B.b2 * matrix[ROW - 3][LINE - 2] - B.b1 * matrix[ROW - 3][LINE - 3])/matrix[ROW - 3][LINE - 4];
243.
244. return B;
245.
246. }
247.
248. //计算高斯拟合参数a，b，c
249. Gaussian_parameter CalcGaussianParameter(B_parameter B)
250. {
251. Gaussian_parameter gaussian;
252. gaussian.a = exp((B.b0 - (B.b1*B.b1/(4*B.b2))));
253. gaussian.b = -1 / B.b2;
254. gaussian.c = -B.b1 / (2 * B.b2);
255. return gaussian;
256. }
257.
258. //从0到250 进行高斯拟合
259. void GuassFit(Gaussian_parameter gaussian)
260. {
261. float y[250];
262. for (int i = 0; i<250;i++)
263. {
264. y[i] = gaussian.a * exp(-(i-gaussian.c)* (i - gaussian.c) /gaussian.b);
265. printf("%f \n",y[i]);
266. }
267. }
268.
________________________________________
程序清单9三次样条插值法拟合测试Main.h
________________________________________
1. #ifndef MAIN_H
2. #define MAIN_H
3.
4. /*该工程作为单片机三次样条插值法的测试工程*/
5.
6. //变量声明放在c文件，函数声明放在h文件
7.
8. //为了配合单片机，变量尽量用float
9.
10. //参考资料：计算机数值方法 65-66 72-73
11.
12.
13. /*********************************************************************************************************
14. * 模块名称: 三次样条插值法测试工程
15. * 摘 要:
16. * 当前版本: 1.00
17. * 作 者: guangjie2333
18. * 完成日期: 2021-4-23
19. *
20. **********************************************************************************************************
21.
22.
23. /*********************************************************************************************************
24. * 包含头文件
25. *********************************************************************************************************/
26.
27.
28. /*********************************************************************************************************
29. * 宏定义
30. *********************************************************************************************************/
31.
32.
33. /*********************************************************************************************************
34. * 内部变量
35. *********************************************************************************************************/
36.
37.
38. /*********************************************************************************************************
39. * 内部函数声明
40. *********************************************************************************************************/
41.
42. /*********************************************************************************************************
43. * 内部函数定义
44. *********************************************************************************************************/
45.
46.
47.
48. #endif // ! MAIN_H
49.
50.
________________________________________
程序清单10三次样条插值法拟合测试Main.c
________________________________________
1. #include "ALG.h"
2.
3. /*内部变量定义*/
4.
5. //拟合曲线的袖带压AD
6. static float cuffPressAD[14] = { 2863,2717,2586,2461,2346,2220,2117,
7. 2053,1950,1882,1807,1718,1643,1557 };
8. //拟合曲线的脉搏波AD
9. static unsigned short paluseAD[14] = { 5622,4007,3004,4229,8464,9874,12125,
10. 10680,8704,6426,4256,4114,3555,3243 };
11. //放气阶数
12. int s_iPluseStepCnt = 14;
13.
14.
15. void main()
16. {
17. float y[250] = { 0 };
18. Node maxNode;
19. float dap, sap, map;
20.
21. for (int i = 0; i < s_iPluseStepCnt; i++)
22. {
23. cuffPressAD[i] = (float)(cuffPressAD[i] * 0.097 - 108); //AD-真实值 线性转换公式
24. paluseAD[i] = (unsigned short)(paluseAD[i] / 16); //脉搏波AD放大了16倍
25. }
26. u16_reserverl(paluseAD, 14);
27. float_reserverl(cuffPressAD, 14);
28.
29. cubicSpline(cuffPressAD, paluseAD, s_iPluseStepCnt, y);
30.
31.
32. //结果显示
33. for (int i = 0; i < 250; i++)
34. {
35. printf("%f \n", y[i]);
36. }
37.
38. maxNode = FindMax(y);
39.
40. //找舒张压
41. for (int i = maxNode.index; i > 0; i--)
42. {
43. if (y[i] <= maxNode.value * 0.8)
44. {
45. printf("dap = %d\n", i);
46. dap = i;
47. break;
48. }
49. }
50.
51. //找收缩压
52. for (int i = maxNode.index; i < 250; i++)
53. {
54. if (y[i] <= maxNode.value * 0.6)
55. {
56. printf("sap = %d\n", i);
57. sap = i;
58. break;
59. }
60. }
61.
62. printf("map = %f \n", dap * 2 / 3 + sap / 3);
63. }
________________________________________
程序清单11三次样条插值法拟合测试ALG.h
________________________________________
1. #ifndef ALG_H
2. #define ALG_H
3.
4. /*该文件是算法函数的头文件*/
5.
6. #include "stdio.h"
7. #include "math.h"
8. #include "stdlib.h"
9.
10. /************************************
11. * 宏定义
12. ************************************/
13.
14.
15. /************************************
16. * 结构体声明
17. ************************************/
18.
19. typedef struct {
20. float value;
21. unsigned short index;
22. }Node;
23.
24.
25.
26.
27. //数组逆排续
28. void u16_reserverl(unsigned short a[], int n);
29. void float_reserverl(float a[], int n);
30.
31. //三次样条插值
32. void cubicSpline(float* cuff, unsigned short* paluse, int s_iPluseStepCnt, float* yFit);
33.
34. //找最大值和最大值所在的数组
35. Node FindMax(float* y);
36. #endif
37.
38.
________________________________________
程序清单12三次样条插值法拟合测试ALG.c
________________________________________
1. #include "ALG.h"
2.
3.
4.
5.
6. /*********************************************************************************************************
7. * 内部函数定义
8. *********************************************************************************************************/
9.
10. /*********************************************************************************************************
11. * 外部函数定义
12. *********************************************************************************************************/
13. void cubicSpline(float* cuff,unsigned short* paluse,int s_iPluseStepCnt, float* yFit)
14. {
15. //参数含义难以注释，详情参考计算机数值方法教材
16. //循环参数
17. int i,j;
18. //找到区间标志
19. char findIndexRange = 0;
20.
21. //矩阵参数 97页
22. float h[20] = { 0 };
23. float lambda[21] = { 0 }; //阶数最多不超过20,为了和教材统一下标
24. float Mu[21] = { 0 };
25. float g[21] = { 0 };
26. //57页
27. float alpha[21] = { 0 };
28. float beta[21] = { 0 };
29. float gamma[21] = { 0 };
30. //58页
31. float u[21] = { 0 };
32. float s[21] = { 0 };
33. float y[21] = { 0 };
34. float x[21] = { 0 };
35.
36. //计算h系数
37. for (i = 0; i < s_iPluseStepCnt;i++)
38. {
39. h[i] = cuff[i+1] - cuff[i];
40. }
41. //计算Mu
42. for (i = 1; i <= s_iPluseStepCnt;i++)
43. {
44. Mu[i] = h[i - 1] / (h[0] + h[i - 1]);
45. }
46. //计算lambda
47. for (i = 1; i <= s_iPluseStepCnt; i++)
48. {
49. lambda[i] = h[0] / (h[0] + h[i - 1]);
50. }
51. //计算g
52. for (i = 1; i <= s_iPluseStepCnt; i++)
53. {
54. g[i] = 3 * (Mu[i]*(paluse[1]- paluse[0])/h[0] + lambda[i]*(paluse[i]- paluse[i-1])/h[i-1]);
55. }
56. //计算alpha & beta & u & s
57. for (i = 1; i <= s_iPluseStepCnt; i++)
58. {
59. if (1 == i)
60. {
61. alpha[i] = 2;
62. beta[i] = Mu[i] / alpha[i];
63. u[i] = lambda[i] / alpha[i];
64. s[i] = Mu[s_iPluseStepCnt];
65. }
66. else if (s_iPluseStepCnt - 1 == i)
67. {
68. alpha[i] = 2 - lambda[i] * beta[i - 1];//alpha要求从1到n
69. gamma[i] = lambda[i];//gamma只要求到n-1
70. u[i] = (Mu[i] - lambda[i] * u[i - 1]) / alpha[i];
71. s[i] = lambda[s_iPluseStepCnt] - s[s_iPluseStepCnt - 1] * beta[s_iPluseStepCnt - 1];
72. }
73. else if (s_iPluseStepCnt == i)
74. {
75. for (j = 1; j < s_iPluseStepCnt - 1; j++)
76. {
77. alpha[i] += s[j] * u[j];
78. }
79. alpha[i] = 2 - alpha[i];
80. }
81. else
82. {
83. alpha[i] = 2 - lambda[i] * beta[i - 1];//alpha要求从1到n
84. beta[i] = Mu[i] / alpha[i];//beta 只要求到n-2
85. gamma[i] = lambda[i];//gamma只要求到n-1
86. u[i] = -lambda[i] * Mu[i - 1] / alpha[i];
87. s[i] = -s[i - 1] * beta[i - 1];
88. }
89.
90. }
91.
92. //计算y
93. for ( i = 1; i <= s_iPluseStepCnt; i++)
94. {
95. if (1 == i)
96. {
97. y[i] = g[i] / alpha[i];
98. }
99. else if (s_iPluseStepCnt == i)
100. {
101. for ( j = 1; j <= s_iPluseStepCnt-1; j++)
102. {
103. y[i] += s[j] * y[j];
104. }
105. y[i] = g[i] - y[i];
106. y[i] /= alpha[i];
107. }
108. else
109. {
110. y[i] = (g[i] - lambda[i] * y[i - 1]) / alpha[i];
111. }
112. }
113.
114.
115. //计算x
116. x[s_iPluseStepCnt] = y[s_iPluseStepCnt];
117. x[0] = x[s_iPluseStepCnt];
118. for (i = s_iPluseStepCnt; i >= 1; i--)
119. {
120. if (s_iPluseStepCnt == i)
121. {
122. x[i] = y[i];
123. }
124. else if (s_iPluseStepCnt - 1 == i)
125. {
126. x[i] = y[i] - u[i]*x[s_iPluseStepCnt];
127. }
128. else
129. {
130. x[i] = y[i] - beta[i] * x[i + 1] - u[i]*x[s_iPluseStepCnt];
131. }
132. }
133.
134. //拟合曲线
135. for ( i = 0; i < 250; i++)
136. {
137. for ( j = 0; j < s_iPluseStepCnt - 1; j++)
138. {
139. if (i >= cuff[j] && i <cuff[j+1])
140. {
141. findIndexRange = 1;
142. break;
143. }
144. }
145.
146. if (1 == findIndexRange)
147. {
148. yFit[i] = (1 - 2 * (i - cuff[j]) / (-h[j])) * ((i - cuff[j + 1]) / (-h[j])) * ((i - cuff[j + 1]) / (-h[j])) * paluse[j] +
149. (1 - 2 * (i - cuff[j + 1]) / h[j]) * ((i - cuff[j]) / (h[j])) * ((i - cuff[j]) / h[j]) * paluse[j + 1] +
150. (i - cuff[j]) * ((i - cuff[j + 1]) / (-h[j])) * ((i - cuff[j + 1]) / (-h[j])) * x[j] +
151. (i - cuff[j + 1]) * ((i - cuff[j]) / (h[j])) * ((i - cuff[j]) / h[j]) * x[j + 1];
152. findIndexRange = 0;
153. }
154. else
155. {
156. yFit[i] = 0; //在采样点范围之外
157. }
158. }
159.
160.
161.
162.
163. }
164.
165.
166.
167. //数组逆排
168. void float_reserverl(float a[], int n)
169. {
170. int m = (n + 1) / 2;
171. int i;
172. int j;
173. float temp;
174. for (i = 0; i < m; i++)
175. {
176. j = n - 1 - i;
177. temp = a[i];
178. a[i] = a[j];
179. a[j] = temp;
180. }
181. }
182.
183. //数组逆排
184. void u16_reserverl(unsigned short a[], int n)
185. {
186. int m = (n + 1) / 2;
187. int i;
188. int j;
189. unsigned short temp;
190. for (i = 0; i < m; i++)
191. {
192. j = n - 1 - i;
193. temp = a[i];
194. a[i] = a[j];
195. a[j] = temp;
196. }
197. }
198.
199. //找最大值和最大值所在的数组
200. Node FindMax(float* y)
201. {
202. int i;
203. Node maxNode;
204. maxNode.value = y[0]; maxNode.index = 0;
205. for (i = 0; i < 250; i++)
206. {
207. if (y[i] > maxNode.value)
208. {
209. maxNode.value = y[i];
210. maxNode.index = i;
211. }
212. }
213. return maxNode;
214. }
215.
________________________________________
扫一扫
08-13
234
234
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
