IT博客汇
  • 首页
  • 精华
  • 技术
  • 设计
  • 资讯
  • 扯淡
  • 权利声明
  • 登录 注册

    python 曲线拟合方法

    Hacper\'s Blog发表于 2021-07-01 00:33:21
    love 0

    背景

    有些产品可能需要用到旋钮来调节音量,比如对讲机,收音机等。如何实现呢?常见的一种实现方案是使用电位器接一个ADC输入,通过ADC检测电压的大小变化来调节音量大小,软件实现上需要定时检测ADC的电压变化,检测周期长短影响音量调节体验,设置的周期长,反应速度慢,会出现突然声音大或者突然声音小的情况,检测时间太频繁,系统又会频繁唤醒,无法进入低功耗状态。

    另一种方式是是使用旋转编码器接几个GPIO,通过旋转编码器可以检测到正转反转,通过其输出信号可以得到旋转角度。

    电位器实物图:

    旋转编码器实物图:

    在某项目中使用ADC+可调电位器的方式实现音量调节功能,电路如下图:

    在使用ADC时发现一个问题,从ADC读取到的电压值并不是实际电压,从ADC读取的电压值范围是0~184mv,测量得到的实际电压为0~1668mv,所以需要对ADC的原始数据做处理,通过软件拟合得到真实电压。遇到问题,解决问题,学习使用python处理曲线拟合的几种方法。

    python 调试环境介绍

    对我的python开发调试环境做个简单介绍,操作系统是win 10,安装的python 是3.8,同时安装了 jupyter lab。

    平时一些简单的python脚本编写、算法验证通常都是在 jupyter lab环境上编码调试。

    为什么使用jupyter lab呢?

    python对我来说是一门辅助性的工具,通常编写的代码量不大,jupyter lab 上可以运行代码,绘图显示方便,同时编写笔记,所以jupyter lab 非常适合我的使用需求,而不需要重量级的IDE开发环境。

    下面是一个来自jupyter lab官网的运行界面图,除了支持python,还支持c++、R、Julia等其他编程语言。

    python 线性拟合方法

    线性拟合需要用到的python包有 numpy、matplotlib 和 scipy,通过pip安装即可,然后导入需要用到的python包

    1
    2
    3
    4
    5
    6
    
    import numpy as np
    import matplotlib as mpl
    import matplotlib.pyplot as plt
    from numpy import polyfit, poly1d
    from scipy.linalg import lstsq
    from scipy.stats import linregress
    

    导入从adc读取到的原始数据和实际测量得到的电量数据,数据量不大,定义为列表,然后绘制ADC原始数据和实际电压的曲线图。

    1
    2
    3
    4
    5
    6
    7
    8
    
    # ADC 采集的原始数据
    adc = [ 192, 333, 372, 399, 421, 435, 450, 453, 459, 462, 456, 436, 419, 386, 340, 199, 184]
    
    # 实际测量的电压数据
    vo = [702, 1212, 1343, 1442, 1525, 1572, 1623, 1638, 1652, 1665, 1647,  1580, 1512, 1395, 1232, 729, 679]
    
    # 绘图
    p = plt.plot(adc, vo, 'rx')
    

    从上图看,原始数据和实际电压数据基本上保持线性关系,这是比较幸运的情况,从图上看应该是一条直线,通过直线拟合的方法来得到实际电压没有问题。

    怎么获取这条直线的参数呢?

    方法1,numpy 提供的多项式拟合函数polyfit计算。

    线性拟合为一阶多项式,一阶多项式$$y = a_1 x + a_0$$拟合成功返回两个系数 $$[a_1, a_0]$$。

    1
    2
    3
    
    coeff1 = polyfit(adc, vo, 1)
    print(coeff1)
    # [ 3.56483415 20.43063565]
    

    使用poly1d生成多项式函数。

    1
    2
    3
    
    f = poly1d(coeff1)
    print(f)
    # 3.565 x + 20.43
    

    然后画出拟合曲线图。

    1
    2
    
    p = plt.plot(adc, f(adc))
    p = plt.plot(adc, vo, 'rx')
    

    方法2,使用scipy提供的线性回归函数linregress求解

    1
    2
    3
    4
    
    # 求解系数
    slope, intercept, r_value, p_value, stderr = linregress(adc, vo)
    print(slope, intercept)
    # 3.5648341453932266 20.430635650877775
    

    绘制拟合曲线图

    1
    2
    3
    4
    5
    6
    
    f = poly1d(coeff)
    print(f)
    # 3.565 x + 20.43
    
    p = plt.plot(adc, vo, 'rx')
    p = plt.plot(adc, f(adc))
    

    方法3,求最小二乘解

    当使用一个 N-1 阶的多项式拟合这 M 个点时,有这样的关系存在$$XC = Y$$ ,

    $$\left[ \begin{matrix} x_0^{N-1} & \dots & x_0 & 1 \\ x_1^{N-1} & \dots & x_1 & 1 \\ \dots & \dots & \dots & \dots \\ x_M^{N-1} & \dots & x_M & 1 \end{matrix}\right] \left[ \begin{matrix} C_{N-1} \\ \dots \\ C_1 \\ C_0 \end{matrix} \right] = \left[ \begin{matrix} y_0 \\ y_1 \\ \dots \\ y_M \end{matrix} \right]$$

    使用一阶多项式求解,即N=2,先将adc扩展为X,vo转换为Y:

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    
    x = np.array(adc)
    Y= np.array(vo).reshape(-1,1)
    
    X = np.hstack((x[:,np.newaxis], np.ones((x.shape[-1],1))))
    print(X[0:8])
    print(Y[0:8])
    
    '''
    [[192.   1.]
     [333.   1.]
     [372.   1.]
     [399.   1.]
     [421.   1.]
     [435.   1.]
     [450.   1.]
     [453.   1.]]
     
    [[ 702]
     [1212]
     [1343]
     [1442]
     [1525]
     [1572]
     [1623]
     [1638]]
    '''
    

    求解:

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    
    C, resid, rank, s = lstsq(X, Y)
    C, resid, rank, s
    
    """
    (array([[ 3.56483415],
            [20.43063565]]),
     array([137.64565162]),
     2,
     array([1.59853058e+03, 9.95236176e-01]))
    """
    

    绘图:

    1
    2
    3
    4
    5
    6
    
    f = poly1d(C.reshape(1,-1)[0])
    print(f)
    # 3.565 x + 20.43
    
    p = plt.plot(adc, vo, 'rx')
    p = plt.plot(adc, f(adc))
    

    非线性拟合

    1. 多项式拟合正弦函数

       1
       2
       3
       4
       5
       6
       7
       8
       9
      10
      11
      12
      13
      14
      15
      16
      17
      18
      19
      20
      
      # 定义sin函数
      x = np.linspace(-3 * np.pi,3 * np.pi,100)
      y = np.sin(x)
      
      # 多项式拟合,从一阶多项式到9阶多项式
      y1 = poly1d(polyfit(x,y,1))
      y3 = poly1d(polyfit(x,y,3))
      y5 = poly1d(polyfit(x,y,5))
      y7 = poly1d(polyfit(x,y,7))
      y9 = poly1d(polyfit(x,y,9))
      
      # 绘图
      p = plt.plot(x, y, 'k-')
      p = plt.plot(x, y1(x))
      p = plt.plot(x, y3(x))
      p = plt.plot(x, y5(x))
      p = plt.plot(x, y7(x))
      p = plt.plot(x, y9(x),"g--")
      
      a = plt.axis([-3 * np.pi, 3 * np.pi, -1.25, 1.25])
      

    2. 拟合自定义曲线

      $$y = a e^{-b sin( f x + \phi)}$$

      定义非线性函数:

      1
      2
      3
      4
      
      def function(x, a , b, f, phi):
          """a function of x with four parameters"""
          result = a * np.exp(-b * np.sin(f * x + phi))
          return result
      
       1
       2
       3
       4
       5
       6
       7
       8
       9
      10
      11
      
      from scipy.optimize import curve_fit
      from scipy.stats import norm
      
      # 绘制原始曲线和加入噪声后曲线
      x = np.linspace(0, 2 * np.pi, 50)
      actual_parameters = [3, 2, 1.25, np.pi / 4]
      y = function(x, *actual_parameters)
      y_noisy = y+ 1.8 * norm.rvs(size=len(x))
      
      p = plt.plot(x, y, 'k-')
      p = plt.plot(x, y_noisy, 'gx')
      

      使用curve_fit求解

      1
      2
      3
      4
      5
      6
      7
      
      p_est, err_est = curve_fit(function, x, y_noisy)
      # 第一个返回的是函数的参数,第二个返回值为各个参数的协方差矩阵
      
      # 绘制结果曲线
      p = plt.plot(x, y_noisy, "gx")
      p = plt.plot(x, y, 'k')
      p = plt.plot(x, function(x, *p_est), "r")
      

      参考资料

      1. 04.04 curve fitting (http://lijin-thu.github.io/04.%20scipy/04.04%20curve%20fitting.html)


沪ICP备19023445号-2号
友情链接