搜索
bottom↓
回复: 222

最小二乘法C算法终极整理版本,绝对原创!

  [复制链接]

出10入23汤圆

发表于 2013-5-22 11:40:35 | 显示全部楼层 |阅读模式
      前两天发了一个关于最小二乘法C语言算法的帖子(http://www.amobbs.com/thread-5535084-1-1.html),这次重新发一个,但是绝对不是属于重复发帖,如果你对算法感兴趣,这个绝对适合你的口味,这次综合了最小二乘法高阶函数拟合(对之前存在的溢出问题有了比较好的解决办法,没有彻底解决),直线拟合(经过高度优化,初始化代码6行,运算代码7行!),抛物线拟合(经过高度优化,初始化代码12行,运算代码20行!)先上代码,工程在附件里面,后续会有进一步的分析(目前还有一个最小二乘法的算法正在研究,思路跟这个版本相差很大,尚不知道是否能实现)





小二

乘法

拟合直线(包含测试函数):
#include "stdafx.h"
#include "process.h"

/*******************************************************************************
阶乘函数宏定义
*******************************************************************************/
#define     my_pow2(x)      ((x) * (x))

/*******************************************************************************
定义一个2行3列的系数矩阵
*******************************************************************************/
static double Para[2][3] = {0.0};

/*******************************************************************************
系数矩阵的初始化,没有撒谎,真的只有6行 + 一个循环
*******************************************************************************/
static int ParaInit(double* X, double* Y, int Amount)
{
        Para[1][1] = Amount;
        for ( ; Amount; Amount--, X++, Y++)
        {
                Para[0][0] += my_pow2(*X);
                Para[0][1] += (*X);
                Para[0][2] += (*X) * (*Y);
                Para[1][2] += (*Y);
        }
        Para[1][0] = Para[0][1];
        return 0;
}

/*******************************************************************************
系数矩阵的运算,没有撒谎,真的只有7行
*******************************************************************************/
static int ParaDeal(void)
{
        Para[0][0] -= Para[1][0] * (Para[0][1] / Para[1][1]);
        Para[0][2] -= Para[1][2] * (Para[0][1] / Para[1][1]);
        Para[0][1] = 0;
        Para[1][2] -= Para[0][2] * (Para[1][0] / Para[0][0]);
        Para[1][0] = 0;
        Para[0][2] /= Para[0][0];
        //Para[0][0] = 1.0;
        Para[1][2] /= Para[1][1];
        //Para[1][1] = 1.0;
        return 0;
}

/***********************************************************************************
从txt文件里读取double型的X,Y数据
***********************************************************************************/
static int GetXY(const char* FileName, double* X, double* Y, int* Amount)
{
        FILE* File = fopen(FileName, "r");
        if (!File)
                return -1;
        for (*Amount = 0; !feof(File); X++, Y++, (*Amount)++)
                if (2 != fscanf(File, (const char*)"%lf %lf", X, Y))
                        break;
        fclose(File);
        return 0;
}

/*******************************************************************************
*******************************************************************************/
int Test(const char* FileName)
{
        int Amount;
        double BufferX[1024], BufferY[1024];
        if (GetXY(FileName, (double*)BufferX, (double*)BufferY, &Amount))
        {
                printf("读取数据文件失败!\r\n");
                return -1;
        }
        printf("读取数据文件成功,共有%d个点!\r\n", Amount);
        ParaInit((double*)BufferX, (double*)BufferY, Amount);
        ParaDeal();
        printf("拟合数据成功,拟合直线为:\r\ny = (%lf) * x + (%lf);\r\n", Para[0][2], Para[1][2]);
        system("pause");
        return 0;
}

int main(int argc, char* argv[])
{
        Test((const char*)"test.txt");
        return 0;
}

拟合抛物线(包含测试函数):
#include "stdafx.h"
#include "process.h"

/*******************************************************************************
阶乘函数宏定义
*******************************************************************************/
#define     my_pow2(x)      ((x) * (x))
#define     my_pow3(x)      ((x) * my_pow2(x))
#define     my_pow4(x)      (my_pow2(x) * my_pow2(x))

/*******************************************************************************
定义一个2行3列的系数矩阵
*******************************************************************************/
static double Para[3][4] = {0.0};

/*******************************************************************************
系数矩阵的初始化
*******************************************************************************/
static int ParaInit(const double* X, const double* Y, int Amount)
{
        Para[2][2] = Amount;
        for ( ; Amount; Amount--, X++, Y++)
        {
                Para[0][0]  +=  my_pow4(*X);
                Para[0][1]  += my_pow3(*X);
                Para[0][2]  += my_pow2(*X);
                Para[1][2]  += (*X);
                Para[0][3]  +=  my_pow2(*X) * (*Y);
                Para[1][3]  +=  (*X) * (*Y);
                Para[2][3]  +=  (*Y);
        }
        Para[1][0] = Para[0][1];
        Para[1][1] = Para[0][2];
        Para[2][0] = Para[1][1];
        Para[2][1] = Para[1][2];
        return 0;
}

/*******************************************************************************
系数矩阵的运算
*******************************************************************************/
static int ParaDeal(void)
{
        //用Para[2][2]把Para[0][2]消成0
        Para[0][0] -= Para[2][0] * (Para[0][2] / Para[2][2]);
        Para[0][1] -= Para[2][1] * (Para[0][2] / Para[2][2]);
        Para[0][3] -= Para[2][3] * (Para[0][2] / Para[2][2]);
        Para[0][2] = 0;
        //用Para[2][2]把Para[1][2] 消成0
        Para[1][0] -= Para[2][0] * (Para[1][2] / Para[2][2]);
        Para[1][1] -= Para[2][1] * (Para[1][2] / Para[2][2]);
        Para[1][3] -= Para[2][3] * (Para[1][2] / Para[2][2]);
        Para[1][2] = 0;
        //用Para[1][1]把Para[0][1]消成0
        Para[0][0] -= Para[1][0] * (Para[0][1] / Para[1][1]);
        Para[0][3] -= Para[1][3] * (Para[0][1] / Para[1][1] );
        Para[0][1] = 0;
        //至此,已经成成为三角矩阵
        //用Para[0][0]把Para[1][0]消成0
        Para[1][3] -= Para[0][3] * (Para[1][0] / Para[0][0]);
        Para[1][0] = 0;
        //用Para[0][0]把Para[2][0]消成0
        Para[2][3] -= Para[0][3] * (Para[2][0] / Para[0][0]);
        Para[2][0] = 0;
        //用Para[1][1]把Para[2][1]消成0
        Para[2][3] -= Para[1][3] * (Para[2][1] / Para[1][1]);
        Para[2][1] = 0;
        //至此,已经成成为对角矩阵
        //把Para[0][0],Para[1][1],Para[2][2]消成1
        Para[0][3] /= Para[0][0];//这个就是最后a的值
        //Para[0][0] = 1.0;
        Para[1][3] /= Para[1][1]; //这个就是最后b的值
        //Para[1][1] = 1.0;
        Para[2][3] /= Para[2][2]; //这个就是最后c的值
        //Para[2][2] = 1.0;
        return 0;
}

/***********************************************************************************
从txt文件里读取double型的X,Y数据
***********************************************************************************/
static int GetXY(const char* FileName, double* X, double* Y, int* Amount)
{
        FILE* File = fopen(FileName, "r");
        if (!File)
                return -1;
        for (*Amount = 0; !feof(File); X++, Y++, (*Amount)++)
                if (2 != fscanf(File, (const char*)"%lf %lf", X, Y))
                        break;
        fclose(File);
        return 0;
}

/*******************************************************************************
*******************************************************************************/
int Test(const char* FileName)
{
        int Amount;
        double BufferX[1024], BufferY[1024];
        if (GetXY(FileName, (double*)BufferX, (double*)BufferY, &Amount))
        {
                printf("读取数据文件失败!\r\n");
                return -1;
        }
        printf("读取数据文件成功,共有%d个点!\r\n", Amount);
        ParaInit((double*)BufferX, (double*)BufferY, Amount);
        ParaDeal();
        printf("拟合数据成功,拟合直线为:\r\ny = (%lf) * x^2 + (%lf) * x + (%lf);\r\n", Para[0][3], Para[1][3], Para[2][3]);
        system("pause");
        return 0;
}

int main(int argc, char* argv[])
{
        Test((const char*)"test.txt");
        return 0;
}

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?注册

x

阿莫论坛20周年了!感谢大家的支持与爱护!!

一只鸟敢站在脆弱的枝条上歇脚,它依仗的不是枝条不会断,而是自己有翅膀,会飞。

出10入23汤圆

 楼主| 发表于 2013-5-22 11:41:10 | 显示全部楼层
自己顶一下

出10入23汤圆

 楼主| 发表于 2013-5-22 11:43:12 | 显示全部楼层
上传比较完整的分析文档
低阶拟合算法的优化:
当阶数比较低时,如直线或者抛物线,解行列式可以优化,不需要做成递推模式,直接强行解就是,而且一般不会出现溢出问题。
还有pow(double x, double y)函数应该自己写替换,原因是在我们的应用中y肯定是整数,而这个函数按照double运算的,效率应该会差很多,所以来个阉割版的pow(double x, int y).在结束低的时候,直接写两个宏定义。
#define     my_pow0(x)      (1.0)
#define     my_pow1(x)      (x)
#define     my_pow2(x)      (my_pow1(x) * my_pow1(x))
#define     my_pow2(x)      (my_pow1(x) * my_pow2(x))
#define     my_pow4(x)      (my_pow2(x) * my_pow2(x))
#define     my_pow5(x)      (my_pow2(x) * my_pow3(x))
#define     my_pow6(x)      (my_pow3(x) * my_pow3(x))
直线拟合的简化:
假设最后拟合直线:
y = a * x + b
残差方程:
(a * x + b – y)^2 = x^2 * a^2 + b^2 + y^2 + 2 * (x * a * b – x * y * a – y * b)
残差方程对a偏导:
x^2 * a + x^1 * b – y * x^1
残差方程对a偏导:
x^1 * a + x^0 * b – y * x^0
当残差方程之和取极值时,偏导为零。
a * ∑x^2 + b * ∑x^1 = ∑(y * x^1)
a * ∑x^1 + b * ∑x^0 = ∑(y * x^0)
(1)首先是定义一个2行3列的系数矩阵,并初始化
double Para[2] [3] = {0.0};
int ParaInit(const double* X, const double* Y, int Amount)
{
        Para[1][1] = Amount;
        for ( ; Amount; Amount--, X++, Y++)
{
        Para[0][0]  +=  my_pow2(*X);
Para[0][1]  +=  my_pow1(*X);
//Para[1][1]  +=  my_pow0(*Y);
Para[0][2]  +=  my_pow1(*X) * my_pow1(*Y);
//Para[1][2]  +=  my_pow0(*X) * my_pow1(*Y);
Para[1][2]  +=  my_pow1(*Y);
}
Para[1][0] = Para[0][1];
return 0;
}

(2)解行列式:
由方程组
a * ∑x^2 + b * ∑x^1 = ∑(y * x^1)
a * ∑x^1 + b * ∑x^0 = ∑(y * x^0)
首先可以看到:
A,Para[0][0]和Para[1][1]肯定大于0
B,Para[0][1] = Para[1][0]
C,        可以证明[2][2]行列式的值
Para[0][0] * Para[1][1] - Para[0][1] * Para[1][0]
(∑x^2) * (∑x^0) - (∑x^1) * (∑x^1)> 0
也即是说,方阵正定,偏导为零的点就是最优值的点
可以看出,第1行的系数比第0行的系数要小很多,所以先采用第1行消第0行,再采用第0行消第1行。
int ParaDeal(void)
{
/*先把Para[0][1]消成0,由于有Para[1][1]肯定大于0,所以用乘法(除法),除法更好一些,数据溢出的问题几乎不存在,适应性会更好*/
Para[0][0] = Para[0][0] - Para[1][0] * (Para[0][1] / Para[1][1];
Para[0][2] = Para[0][2] - Para[1][2] * (Para[0][1] / Para[1][1];
Para[0][1] = 0;
//再把Para[1][0] 消成0
Para[1][2] = Para[1][2] - Para[0][2] * (Para[1][0] / Para[0][0]);
Para[1][0] = 0;
//把Para[0][0],Para[1][1] 消成1
Para[0][2] /= Para[0][0];//这个就是最后a的值
Para[0][0] = 1.0;
Para[1][2] /= Para[1][1]; //这个就是最后b的值
Para[1][1] = 1.0;
return 0;
}
至此,一阶直线拟合的推导过程和C代码全部完毕,可以看出运算量可以说很小,两个函数均只需调用一次即可得到结果,在源数据的个数和大小不是特别大的时候,完全可以移植到MCU中进行运算,实用价值比较高。
再来个二阶抛物线拟合的优化:
前期推导和一阶直线差不多,略去,系数矩阵方程组:
a * ∑x^4 + b * ∑x^3 + c * ∑x^2 = ∑(y * x^2)
a * ∑x^3 + b * ∑x^2 + c * ∑x^1 = ∑(y * x^1)
a * ∑x^2 + b * ∑x^1 + c * ∑x^0 = ∑(y * x^0)
(1)首先是定义一个3行4列的系数矩阵,并初始化
double Para[3] [4] = {0.0};
int ParaInit(const double* X, const double* Y, int Amount)
{
        Para[1][1] = Amount;
        for ( ; Amount; Amount--, X++, Y++)
{
Para[0][0]  +=  my_pow4(*X);
Para[0][1]  += my_pow3(*X);
Para[0][2]  += my_pow2(*X);
Para[1][2]  += my_pow1(*X);
//Para[2][2] +=  my_pow0(*X);
Para[0][3]  +=  my_pow2(*X) * my_pow1(*Y);
Para[1][3]  +=  my_pow1(*X) * my_pow1(*Y);
Para[2][3]  +=  my_pow0(*X) * my_pow1(*Y);
}
Para[1][0] = Para[0][1];
Para[1][1] = Para[0][2];
Para[2][0] = Para[1][1];
Para[2][1] = Para[1][2];
return 0;
}

a * ∑x^4 + b * ∑x^3 + c * ∑x^2 = ∑(y * x^2)
a * ∑x^3 + b * ∑x^2 + c * ∑x^1 = ∑(y * x^1)
a * ∑x^2 + b * ∑x^1 + c * ∑x^0 = ∑(y * x^0)
首先可以看到:
A,Para[0][0],Para[1][1], Para[2][2]肯定大于0
B,左边3*3的方阵是对称矩阵
D,        可以证明[3][3]行列式的值 > 0,也即是说,方阵正定,偏导为零的点就是最优值的点
可以看出,第k + 1行的系数比第k行的系数要小很多,所以先采用第k +1行消第k行,再采用第k行消第k + 1行。

int ParaDeal(void)
{
//用Para[2][2]把Para[0][2]消成0
Para[0][0] = Para[0][0] - Para[2][0] * (Para[0][2] / Para[2][2]);
Para[0][1] = Para[0][1] - Para[2][1] * (Para[0][2] / Para[2][2]);
Para[0][3] = Para[0][3] - Para[2][3] * (Para[0][2] / Para[2][2]);
Para[0][2] = 0;
//用Para[2][2]把Para[1][2] 消成0
Para[1][0] = Para[1][0] - Para[2][0] * (Para[1][2] / Para[2][2]);
Para[1][1] = Para[1][1] - Para[2][1] * (Para[1][2] / Para[2][2]);
Para[1][3] = Para[1][3] - Para[2][3] * (Para[1][2] / Para[2][2]);
Para[1][2] = 0;
//用Para[1][1]把Para[0][1]消成0
Para[0][0] = Para[0][0] - Para[1][0] * (Para[0][1] / Para[1][1]);
Para[0][3] = Para[0][3] - Para[1][3] * (Para[0][1] / Para[1][1] );
Para[0][1] = 0;
//至此,已经成成为三角矩阵
//用Para[0][0]把Para[1][0]消成0
Para[1][3] = Para[1][3] - Para[0][3] * (Para[1][0] / Para[0][0]);
Para[1][0] = 0;
//用Para[0][0]把Para[2][0]消成0
Para[2][3] = Para[2][3] - Para[0][3] * (Para[2][0] / Para[0][0]);
Para[2][0] = 0;
//用Para[1][1]把Para[2][1]消成0
Para[2][3] = Para[2][3] - Para[1][3] * (Para[2][1] / Para[1][1]);
Para[2][1] = 0;
//至此,已经成成为对角矩阵
//把Para[0][0],Para[1][1],Para[2][2]消成1
Para[0][3] /= Para[0][0];//这个就是最后a的值
Para[0][0] = 1.0;
Para[1][3] /= Para[1][1]; //这个就是最后b的值
Para[1][1] = 1.0;
Para[2][3] /= Para[2][2]; //这个就是最后c的值
Para[2][2] = 1.0;
return 0;
}

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?注册

x

出10入23汤圆

 楼主| 发表于 2013-5-22 11:54:03 | 显示全部楼层

上传拟合高阶曲线的VC6.0工程,经测试在20个数据,6阶的拟合中与matlab拟合结果基本吻合,没有出现数据溢出,由于取消了限幅函数,数据量大(线性增加),阶数高(觅增长),数据大(指数增长)的时候估计还是会溢出,欢迎测试

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?注册

x

出0入0汤圆

发表于 2013-5-22 11:57:04 | 显示全部楼层
LZ有没有指数曲线拟合的C程序?

出10入23汤圆

 楼主| 发表于 2013-5-22 11:59:41 | 显示全部楼层
无级电工 发表于 2013-5-22 11:57
LZ有没有指数曲线拟合的C程序?

目前没有,你要的话可以试着帮你推一个,这个多项式拟合也就前天才开始推导,但指数比这个难一点,需要比较严谨的数学方法证明才行

出10入23汤圆

 楼主| 发表于 2013-5-22 12:09:49 | 显示全部楼层
无级电工 发表于 2013-5-22 11:57
LZ有没有指数曲线拟合的C程序?

如果你仔细看那个一阶或二阶的推导过程,应该你也能把指数的推出来,只不过把多项式换成指数函数,(最好经过严格的数学证明,或者有参考工具可以校验,比如matlab)

出10入23汤圆

 楼主| 发表于 2013-5-22 12:13:32 | 显示全部楼层
zouzhichao 发表于 2013-5-22 12:09
如果你仔细看那个一阶或二阶的推导过程,应该你也能把指数的推出来,只不过把多项式换成指数函数,(最好 ...

说错了,指数跟这个还是很不相同的,第一步无法得到一个确定的系数矩阵,所以得换思路才行

出0入0汤圆

发表于 2013-5-22 12:20:19 | 显示全部楼层
太感谢楼主了,雪中送炭。
不过楼主你能否弄个 递推抛物拟合。就是不断的有新数据进入,老数据推出;或者对老数据遗忘处理?

出10入23汤圆

 楼主| 发表于 2013-5-22 12:54:18 | 显示全部楼层
kmani 发表于 2013-5-22 12:20
太感谢楼主了,雪中送炭。
不过楼主你能否弄个 递推抛物拟合。就是不断的有新数据进入,老数据推出;或者对 ...

不太明白你的意思,是动态的拟合吗?每次的XY数据有类似FIFO的更新?

出10入23汤圆

 楼主| 发表于 2013-5-22 12:55:45 | 显示全部楼层
kmani 发表于 2013-5-22 12:20
太感谢楼主了,雪中送炭。
不过楼主你能否弄个 递推抛物拟合。就是不断的有新数据进入,老数据推出;或者对 ...

抛物线拟合的代码和数学推导都比较简单,你试着改一改吧

出0入0汤圆

发表于 2013-5-22 12:56:31 | 显示全部楼层
zouzhichao 发表于 2013-5-22 12:54
不太明白你的意思,是动态的拟合吗?每次的XY数据有类似FIFO的更新?

嗯,是的。
因为我想递推拟合的话,计算量可能会小一些。
直接拟合每次都重新计算了。

出0入0汤圆

发表于 2013-5-22 12:58:06 | 显示全部楼层
zouzhichao 发表于 2013-5-22 12:55
抛物线拟合的代码和数学推导都比较简单,你试着改一改吧

直接拟合的我会,就是递推的不行,不太会推公式...

出0入0汤圆

发表于 2013-5-22 13:01:16 | 显示全部楼层
顶一个先

出10入23汤圆

 楼主| 发表于 2013-5-22 13:02:35 | 显示全部楼层
kmani 发表于 2013-5-22 12:56
嗯,是的。
因为我想递推拟合的话,计算量可能会小一些。
直接拟合每次都重新计算了。 ...

你的意思是,假设这里有1000个数据,但是这1000个不是立马就有的,而是一边运行,一边产生的?
拟合1——100个数
拟合2——101个数
拟合3——102个数
拟合4——103个数

是这样吗?那么当前的拟合结果和上一次拟合结果需要关联吗?

出10入23汤圆

 楼主| 发表于 2013-5-22 13:05:36 | 显示全部楼层
kmani 发表于 2013-5-22 12:58
直接拟合的我会,就是递推的不行,不太会推公式...

明白了,你要拟合1000个数据,但一次只拟合100个,但是之前的拟合结果会影响这次的拟合结果,最后的效果是可以做到实时拟合,并且不需要计算1000个数的拟合,只需要计算100个数的拟合,是这样吧?

出10入23汤圆

 楼主| 发表于 2013-5-22 13:06:40 | 显示全部楼层
那这个的逻辑关系就复杂了,暂时没有好的思路

出10入23汤圆

 楼主| 发表于 2013-5-22 13:07:12 | 显示全部楼层
至少先数学推导才敢下笔

出0入0汤圆

发表于 2013-5-22 17:45:08 | 显示全部楼层
本帖最后由 kmani 于 2013-5-22 18:02 编辑
zouzhichao 发表于 2013-5-22 13:05
明白了,你要拟合1000个数据,但一次只拟合100个,但是之前的拟合结果会影响这次的拟合结果,最后的效果 ...


嗯,就是这样的,每次我只用100个数据拟合。下一次计算时进来一个新的数据,同时去掉一个最旧的数据,再次拟合。
最终的目的就是实时拟合,新数据不断地进来,旧数据不断地出去,保证每次拟合都是100个数据。
我为啥想用递推算法呢,因为假如使用楼主的算法,每次都计算了100个点的拟合,下一次拟合还是计算100个点,就没有使用上一次计算的部分结果,这样增大了计算量。
就跟拉格朗日插值一样,只要新增加点之后,所有的项都要重新计算。但是牛顿插值就只用计算新的项就可以了,不用全部从头算。但是我的想法是有了新的数据后,还有自动去掉旧的数据。
假如用递推方式计算,我估计只能采用加遗忘因子的方法了,不能准确的使用100个数据。

出0入0汤圆

发表于 2013-5-22 18:21:29 | 显示全部楼层
加精吧。很好

出10入23汤圆

 楼主| 发表于 2013-5-22 18:25:15 | 显示全部楼层
kmani 发表于 2013-5-22 17:45
嗯,就是这样的,每次我只用100个数据拟合。下一次计算时进来一个新的数据,同时去掉一个最旧的数据,再 ...

要实现精确的数学模型来实现有点难度,可能模糊地实现(不保证数学上的最优点,但是更接近最优点)可能好实现一些

出10入23汤圆

 楼主| 发表于 2013-5-22 18:25:38 | 显示全部楼层
jacobson 发表于 2013-5-22 18:21
加精吧。很好

多谢夸奖

出10入23汤圆

 楼主| 发表于 2013-5-22 18:34:15 | 显示全部楼层
kmani 发表于 2013-5-22 17:45
嗯,就是这样的,每次我只用100个数据拟合。下一次计算时进来一个新的数据,同时去掉一个最旧的数据,再 ...

多交流哈,递推的拟合你要是有什么好的思路教教我,突然对这个很有兴趣

出10入23汤圆

 楼主| 发表于 2013-5-22 18:38:30 | 显示全部楼层
目前临近毕业了,杂七杂八的事情多了,沉不下心看书,等忙完了再思考思考上面提到的指数曲线拟合和递推拟合

出0入0汤圆

发表于 2013-5-22 23:22:26 | 显示全部楼层
厉害,向楼主学习

出10入23汤圆

 楼主| 发表于 2013-5-22 23:26:48 | 显示全部楼层
xiaowenshao 发表于 2013-5-22 23:22
厉害,向楼主学习

出0入0汤圆

发表于 2013-5-23 07:41:35 | 显示全部楼层
厉害,向楼主学习最佩服有知识的人了

出0入0汤圆

发表于 2013-5-23 08:29:35 | 显示全部楼层
收下      

出0入0汤圆

发表于 2013-5-23 08:31:57 | 显示全部楼层
收藏 谢谢。

出0入0汤圆

发表于 2013-5-23 09:40:16 | 显示全部楼层
楼主很佩服你啊!学问做的这么好!还这么谦虚!

出10入23汤圆

 楼主| 发表于 2013-5-23 09:46:03 | 显示全部楼层
jsszdfdn 发表于 2013-5-23 09:40
楼主很佩服你啊!学问做的这么好!还这么谦虚!

学问真不好,这点实在惭愧,现在快毕业了还差十来个学分没修完 ,这个也不是用了什么高深的算法,推导过程不难,只是要花点时间写草稿演算

出0入0汤圆

发表于 2013-5-23 09:53:25 | 显示全部楼层
做个标记,有时间研究

出0入0汤圆

发表于 2013-5-23 09:59:08 | 显示全部楼层
徐士良 常用算法程序集
里面都有

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?注册

x

出10入23汤圆

 楼主| 发表于 2013-5-23 10:10:29 | 显示全部楼层
biansf2001 发表于 2013-5-23 09:59
徐士良 常用算法程序集
里面都有

非常感谢

出10入23汤圆

 楼主| 发表于 2013-5-23 10:13:08 | 显示全部楼层
biansf2001 发表于 2013-5-23 09:59
徐士良 常用算法程序集
里面都有

浏览了一下目录,好书!

出0入0汤圆

发表于 2013-5-23 11:17:17 | 显示全部楼层
ma入库 好,有空看看

出0入0汤圆

发表于 2013-5-29 10:14:07 | 显示全部楼层
收下,不容错过!学习

出10入23汤圆

 楼主| 发表于 2013-5-29 10:22:03 | 显示全部楼层
wangfei1956 发表于 2013-5-29 10:14
收下,不容错过!学习

谢谢支持

出0入0汤圆

发表于 2013-5-29 11:01:46 | 显示全部楼层
不错的帖子,好好学习下

出10入23汤圆

 楼主| 发表于 2013-6-3 08:29:29 | 显示全部楼层
无级电工 发表于 2013-5-22 11:57
LZ有没有指数曲线拟合的C程序?


拟合指数函数的示例,是这样子的吗?

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?注册

x

出0入0汤圆

发表于 2013-6-3 08:33:27 | 显示全部楼层
zouzhichao 发表于 2013-6-3 08:29
拟合指数函数的示例,是这样子的吗?

看起来效果不错啊。

出10入23汤圆

 楼主| 发表于 2013-6-3 08:46:27 | 显示全部楼层
无级电工 发表于 2013-6-3 08:33
看起来效果不错啊。

用的迭代,有个问题,有时候不收敛,数据也不能太大

出10入23汤圆

 楼主| 发表于 2013-6-3 08:53:24 | 显示全部楼层
无级电工 发表于 2013-6-3 08:33
看起来效果不错啊。

方法还是一阶偏导为零,但是没有一个确定的系数矩阵,不能直接求解行列式
解决办法是:增加一个二阶偏导(一个3*3的矩阵)
(1)给a,b,c一组初始值,求出一阶偏导的值,二阶偏导矩阵
(2)一阶偏导 / 二阶偏导矩阵,得到da,db,dc。(求解行列式)
(3)a = a - da; b = b - db; c = c - dc; 更新a,b,c的值,回到(1),进入迭代

出0入0汤圆

发表于 2013-6-3 08:54:07 | 显示全部楼层
zouzhichao 发表于 2013-6-3 08:46
用的迭代,有个问题,有时候不收敛,数据也不能太大

这倒是个问题。不急,将来搞成熟了再发布也行啊。 网上指数函数拟合的C程序没见过,你这应该算头一份了。

出10入23汤圆

 楼主| 发表于 2013-6-3 09:09:18 | 显示全部楼层
本帖最后由 zouzhichao 于 2013-6-3 09:13 编辑
无级电工 发表于 2013-6-3 08:33
看起来效果不错啊。


m代码分四个:
(1)求解一阶导数GetYijieDaoshu.m
%得到一阶导数的值,公式如下
%2 * exp(b * x) * (c - y + a * exp(b * x))
%2 * a * x * exp(b * x) * (c - y + a * exp(b * x))
%2 * c - 2 * y + 2 * a * exp(b * x)
function Res = GetYijieDaoshu(Xy, Para)
Res = zeros(3, 1);
for i = 1 : size(Xy, 1)
    Res(1, 1) = Res(1, 1) + 2 * exp(Para(2) * Xy(i, 1)) * (Para(3) - Xy(i, 2) + Para(1) * exp(Para(2) * Xy(i, 1)));
    Res(2, 1) = Res(2, 1) + 2 * Para(1) * Xy(i, 1) * exp(Para(2) * Xy(i, 1)) * (Para(3) - Xy(i, 2) + Para(1) * exp(Para(2) * Xy(i, 1)));
    Res(3, 1) = Res(3, 1) + 2 * Para(3) - 2 * Xy(i, 2) + 2 * Para(1) * exp(Para(2) * Xy(i, 1));
end
end

(2)求解二阶导数GetErjieDaoshu.m
%得到二阶导数的值,公式如下
%2 * exp(2 * b * x)
%2 * x * exp(b * x) * (c - y + a * exp(b * x)) + 2 * a * x * exp(2 * b * x)
%2 * exp(b * x)
%2 * x * exp(b * x) * (c - y + a * exp(b * x)) + 2 * a * x * exp(2 * b * x)
%2 * a^2 * x^2 * exp(2 * b * x) + 2 * a * x^2 * exp(b * x) * (c - y + a * exp(b * x))
%2 * a * x * exp(b * x)
%2 * exp(b * x)
%2 * a * x * exp(b * x)
%2
function Res = GetErjieDaoshu(Xy, Para)
Res = zeros(3, 3);
for i = 1 : size(Xy, 1)
    Res(1, 1) = Res(1, 1) + 2 * exp(2 * Para(2) * Xy(i, 1));
    Res(1, 2) = Res(1, 2) + 2 * Xy(i, 1) * exp(Para(2) * Xy(i, 1)) * (Para(3) - Xy(i, 2) + Para(1) * exp(Para(2) * Xy(i, 1))) + 2 * Para(1) * Xy(i, 1) * exp(2 * Para(2) * Xy(i, 1));
    Res(1, 3) = Res(1, 3) + 2 * exp(Para(2) * Xy(i, 1));
    Res(2, 1) = Res(2, 1) + 2 * Xy(i, 1) * exp(Para(2) * Xy(i, 1)) * (Para(3) - Xy(i, 2) + Para(1) * exp(Para(2) * Xy(i, 1))) + 2 * Para(1) * Xy(i, 1) * exp(2 * Para(2) * Xy(i, 1));
    Res(2, 2) = Res(2, 2) + 2 * Para(1)^2 * Xy(i, 1)^2 * exp(2 * Para(2) * Xy(i, 1)) + 2 * Para(1) * Xy(i, 1)^2 * exp(Para(2) * Xy(i, 1)) * (Para(3) - Xy(i, 2) + Para(1) * exp(Para(2) * Xy(i, 1)));
    Res(2, 3) = Res(2, 3) + 2 * Para(1) * Xy(i, 1) * exp(Para(2) * Xy(i, 1));
    Res(3, 1) = Res(3, 1) + 2 * exp(Para(2) * Xy(i, 1));
    Res(3, 2) = Res(3, 2) + 2 * Para(1) * Xy(i, 1) * exp(Para(2) * Xy(i, 1));
    Res(3, 3) = Res(3, 3) + 2;
end
end

(3)迭代计算CalPara.m
%计算系数的值
function Res = CalPara(Xy, ParaInit, N)
Res = ParaInit;
for i = 1 : N
    ParaN = GetYijieDaoshu(Xy, Res);
    ParaNN = GetErjieDaoshu(Xy, Res);
    Res = Res - 1 * (ParaNN \ ParaN);
end
end

(4)测试Test.m
clear;
clc;
%目标函数y = a * exp(b * x) + c
%预定参数
a = 1.7;
b = 0.2;
c = -4.6;
%待拟合的点,带上(rand(size(X)) - 0.5) * 0.2的随机偏差
X = -3 : 0.1 : 2.9;
Y = a * exp(b * X) + c + (rand(size(X)) - 0.5) * 0.2;
%拟合,注意需要给出a, b, c的初始值,不要离得太远,否则会发散导致迭代失败
%这里我们取[1.0, 1.0, -5.2]',跟预定参数[1.7, 0.2, -4.6]'距离不是很远,迭代200次
ParaABC = CalPara([X', Y'], [1.0, 1.0, -5.2]', 200);
%绘图比较
X1 = min(X) : (max(X)- min(X))/9999 : max(X);
Z = ParaABC(1) * exp(ParaABC(2) * X1) + ParaABC(3);
plot(X, Y, '*', X1, Z);

测试图:


本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?注册

x

出10入23汤圆

 楼主| 发表于 2013-6-3 09:16:59 | 显示全部楼层
无级电工 发表于 2013-6-3 08:54
这倒是个问题。不急,将来搞成熟了再发布也行啊。 网上指数函数拟合的C程序没见过,你这应该算头一份了 ...

目前还没有写C,但是没有用matlab现成的拟合函数,最高级的也只是用了一个矩阵除法(解行列式用),结合之前高阶多项式拟合的时候写的那个纯C的行列式求解代码,很容易纯C实现

出0入0汤圆

发表于 2013-6-3 09:23:37 | 显示全部楼层
好帖子,顶一下!

出10入23汤圆

 楼主| 发表于 2013-6-3 09:24:26 | 显示全部楼层
gfy200866 发表于 2013-6-3 09:23
好帖子,顶一下!

多谢帮顶

出0入0汤圆

发表于 2013-6-3 09:25:05 | 显示全部楼层
zouzhichao 发表于 2013-6-3 08:29
拟合指数函数的示例,是这样子的吗?

这样的,曲线该怎么实现算法。不知楼主有没有好的原创案例。

出10入23汤圆

 楼主| 发表于 2013-6-3 09:27:23 | 显示全部楼层
gfy200866 发表于 2013-6-3 09:25
这样的,曲线该怎么实现算法。不知楼主有没有好的原创案例。

?没听明白

出0入0汤圆

发表于 2013-6-3 09:27:33 | 显示全部楼层
如果不是一个标准的指数函数。但也接近。LZ有什么好的办法减少失真?

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?注册

x

出10入23汤圆

 楼主| 发表于 2013-6-3 09:32:13 | 显示全部楼层
gfy200866 发表于 2013-6-3 09:27
如果不是一个标准的指数函数。但也接近。LZ有什么好的办法减少失真?


*的点不是指数函数了,加了随机偏差,拟合的话就是找一条标准的指数函数曲线,使之跟那些点尽可能的近,那条实线就是拟合出来的

出0入0汤圆

发表于 2013-6-3 09:34:27 | 显示全部楼层
zouzhichao 发表于 2013-6-3 09:32
*的点不是指数函数了,加了随机偏差,拟合的话就是找一条标准的指数函数曲线,使之跟那些点尽可能的近, ...

LZ的意思是需要先找一个接近实际需要的标准函数,然后在想办法进行拟合吗?

出0入0汤圆

发表于 2013-6-3 09:35:42 | 显示全部楼层
这个就是所说的最小二乘法吗?以前没接触过这个算法。有没有相关的资料分享一下吗?谢谢!

出0入0汤圆

发表于 2013-6-3 09:36:05 | 显示全部楼层
牛....

出10入23汤圆

 楼主| 发表于 2013-6-3 09:39:21 | 显示全部楼层
本帖最后由 zouzhichao 于 2013-6-3 09:41 编辑
gfy200866 发表于 2013-6-3 09:35
这个就是所说的最小二乘法吗?以前没接触过这个算法。有没有相关的资料分享一下吗?谢谢! ...


资料全部在上面,最后那一部分还是今天早上刚出炉的,别的资料我只能说大学课本+自己的思维了

出0入0汤圆

发表于 2013-6-3 09:46:01 | 显示全部楼层
支持一下~

出10入23汤圆

 楼主| 发表于 2013-6-3 09:47:39 | 显示全部楼层
x9fish 发表于 2013-6-3 09:46
支持一下~

感谢支持

出10入23汤圆

 楼主| 发表于 2013-6-3 11:07:02 | 显示全部楼层
本帖最后由 zouzhichao 于 2013-6-3 11:08 编辑
gfy200866 发表于 2013-6-3 09:35
这个就是所说的最小二乘法吗?以前没接触过这个算法。有没有相关的资料分享一下吗?谢谢! ...


楼上有网友推荐的那本算法书不错
33楼,徐士良 《常用算法程序集》

出0入0汤圆

发表于 2013-6-3 11:45:59 | 显示全部楼层
非常感谢楼主!入库收藏。

出10入23汤圆

 楼主| 发表于 2013-6-3 11:53:48 | 显示全部楼层
woshimajia222 发表于 2013-6-3 11:45
非常感谢楼主!入库收藏。

还不太完善,有待改进

出0入0汤圆

发表于 2013-6-3 15:06:11 | 显示全部楼层
zouzhichao 发表于 2013-6-3 11:07
楼上有网友推荐的那本算法书不错
33楼,徐士良 《常用算法程序集》

非常感谢!

出0入0汤圆

发表于 2013-6-3 15:18:27 来自手机 | 显示全部楼层
学习 谢谢

出70入0汤圆

发表于 2013-6-13 18:08:34 | 显示全部楼层
楼主太厉害了,虽然看不懂,先mark!

出0入0汤圆

发表于 2013-6-14 15:58:10 | 显示全部楼层
非常感谢。

出0入0汤圆

发表于 2013-8-29 08:56:05 | 显示全部楼层
谢谢分享,下载看看

出0入0汤圆

发表于 2013-8-31 13:39:55 | 显示全部楼层
MARK一下
头像被屏蔽

出0入0汤圆

发表于 2013-9-19 15:09:45 | 显示全部楼层
COOL !

出0入0汤圆

发表于 2013-9-19 15:13:15 | 显示全部楼层
收藏!~                                 

出0入0汤圆

发表于 2013-9-19 15:25:56 | 显示全部楼层
每次看算法都头疼。。。

出0入0汤圆

发表于 2013-9-19 15:40:45 | 显示全部楼层
这个帖子不错哦

出0入0汤圆

发表于 2013-9-20 04:54:01 | 显示全部楼层
本帖最后由 blue.fox 于 2013-9-20 04:56 编辑

LZ代码好多啊。
LZ这个是在整非线性最小二乘法的叠代吧?有在DSP上实际测试过执行时间没?
另外,LZ矩阵求逆在DSP里怎么实现的?没看到相关的部分
谢谢

出0入0汤圆

发表于 2013-9-20 08:03:23 | 显示全部楼层
mark

出10入23汤圆

 楼主| 发表于 2013-9-23 12:47:14 | 显示全部楼层
blue.fox 发表于 2013-9-20 04:54
LZ代码好多啊。
LZ这个是在整非线性最小二乘法的叠代吧?有在DSP上实际测试过执行时间没?
另外,LZ矩阵求 ...

是的,没有在DSP中跑,都是在GCC和VC6.0,WinAVR GCC和GCC中编译测试的,前面的发帖是多项式拟合,里面有矩阵求逆的代码(解一个方程组),后面的发帖是为了解决其中有一个坛友提出的指数曲线拟合的问题,用的是迭代的方法,测试了指数函数拟合以及正弦函数拟合,都还行,只是对迭代初始化位置比较敏感,尤其是正弦函数的频率参数非常敏感,改进方法暂时没有新的思路

出0入0汤圆

发表于 2013-9-23 12:55:49 | 显示全部楼层
zouzhichao 发表于 2013-9-23 12:47
是的,没有在DSP中跑,都是在GCC和VC6.0,WinAVR GCC和GCC中编译测试的,前面的发帖是多项式拟合,里面有 ...

如果是多项式线性拟合的话,(A^T * A)^-1 * A^T部分是个常数矩阵,可以先用matlab算好,然后C语言就只需要做一次矩阵乘法了

指数曲线的迭代我也分析过。里面的求逆,判断奇异的运算量太大了。而且的确对初始位置比较敏感,DSP里面难以实现。最后换的别的方法。

出0入0汤圆

发表于 2013-10-2 15:47:31 | 显示全部楼层
好贴!!!

出0入0汤圆

发表于 2013-10-16 19:54:01 | 显示全部楼层
matlab不行吗!                              

出10入23汤圆

 楼主| 发表于 2013-10-16 22:18:56 | 显示全部楼层
笑傲江湖 发表于 2013-10-16 19:54
matlab不行吗!

从来就没说matlab不行啊

出0入0汤圆

发表于 2013-10-17 07:38:36 来自手机 | 显示全部楼层
很厉害的感觉

出0入0汤圆

发表于 2014-1-13 12:07:12 | 显示全部楼层
强帖留名,刚好有需要,谢谢楼主分享.

出0入0汤圆

发表于 2014-1-13 14:32:28 | 显示全部楼层
楼主厉害,赞一个!

出0入0汤圆

发表于 2014-1-21 00:23:05 | 显示全部楼层
学习学习。。。

出0入0汤圆

发表于 2014-1-21 08:57:35 | 显示全部楼层
mark 一下

出0入0汤圆

发表于 2014-1-21 10:43:40 | 显示全部楼层
收藏成功

出0入0汤圆

发表于 2014-1-21 12:10:29 | 显示全部楼层
先赞一个,有需要时再来下载。

出0入0汤圆

发表于 2014-1-21 12:50:37 | 显示全部楼层
顶顶。。。。收藏了

出0入4汤圆

发表于 2014-1-21 12:53:25 | 显示全部楼层
不错,看哈

出0入0汤圆

发表于 2014-1-21 13:52:07 | 显示全部楼层
mark下先,以后可能用到

出0入0汤圆

发表于 2014-1-21 14:34:30 | 显示全部楼层
小白学习中……

出0入0汤圆

发表于 2014-1-21 15:55:04 | 显示全部楼层
mark

出0入0汤圆

发表于 2014-1-21 21:38:23 | 显示全部楼层
晕死,matlab这些都还给老师了。

出0入0汤圆

发表于 2014-1-21 22:23:23 | 显示全部楼层
膜拜大神、、

出0入0汤圆

发表于 2014-1-22 08:37:58 | 显示全部楼层
mark,有空再看

出10入23汤圆

 楼主| 发表于 2014-1-22 12:36:01 | 显示全部楼层

惭愧,非大神

出675入8汤圆

发表于 2014-1-22 12:43:33 | 显示全部楼层
绝对的好贴

出0入0汤圆

发表于 2014-1-22 14:14:37 | 显示全部楼层
make下,回去学习

出0入0汤圆

发表于 2014-1-23 09:31:55 | 显示全部楼层
写得好,值得保留

出0入0汤圆

发表于 2014-1-24 13:14:57 | 显示全部楼层
make下,回去学习

出0入0汤圆

发表于 2014-1-24 15:06:54 | 显示全部楼层
mark,good

出0入0汤圆

发表于 2014-2-18 10:05:18 | 显示全部楼层
标记mark
回帖提示: 反政府言论将被立即封锁ID 在按“提交”前,请自问一下:我这样表达会给举报吗,会给自己惹麻烦吗? 另外:尽量不要使用Mark、顶等没有意义的回复。不得大量使用大字体和彩色字。【本论坛不允许直接上传手机拍摄图片,浪费大家下载带宽和论坛服务器空间,请压缩后(图片小于1兆)才上传。压缩方法可以在微信里面发给自己(不要勾选“原图),然后下载,就能得到压缩后的图片】。另外,手机版只能上传图片,要上传附件需要切换到电脑版(不需要使用电脑,手机上切换到电脑版就行,页面底部)。
您需要登录后才可以回帖 登录 | 注册

本版积分规则

手机版|Archiver|amobbs.com 阿莫电子技术论坛 ( 粤ICP备2022115958号, 版权所有:东莞阿莫电子贸易商行 创办于2004年 (公安交互式论坛备案:44190002001997 ) )

GMT+8, 2024-3-29 12:58

© Since 2004 www.amobbs.com, 原www.ourdev.cn, 原www.ouravr.com

快速回复 返回顶部 返回列表