使用 lmfit 将信号拟合为方波形式

编程问答 2022-07-11

问题描述

我想使用 lmfit 模块拟合遵循方波模式的数据。我遇到的问题是，在拟合代码时，我对频率或数据偏移的初始猜测没有变化。我阅读了有关此问题的常见问题解答 (https://lmfit.github.io/lmfit-py/faq.html)，它说当变量的微小变化不会影响拟合的残差时，这些变量根本没有变化，就会发生这种情况。

另外，有没有什么好的方法来估计信号的频率以产生一个好的起始猜测？

from matplotlib import pyplot as plt
import numpy as np
from lmfit import Model
from numpy import random

def square_wave(x,f,a,h,phi):
    return h + a * np.sign(np.cos( 2*np.pi*f*x + phi ))

def analyse(x,y):
    supermodel = Model(square_wave)

    print('Parameter names: {}'.format(supermodel.param_names))
    print('Indipendent variables: {}'.format(supermodel.independent_vars))

    params = supermodel.make_params(f = 5,a = (np.max(y)-np.min(y))/2,h = np.mean(y),phi = 0 )
    result = supermodel.fit(y,params,x=x)

    print(result.fit_report())

    plt.plot(x,y,'b')
    plt.plot(x,result.init_fit,'k--',label='initial fit')
    plt.plot(x,result.best_fit,'r-',label='best fit')
    plt.legend(loc='best')
    plt.show()
                                    
t = np.linspace(0,1,1000)
f = 6

y = square_wave(t,a = 1,h = 0,phi = 0) + random.normal(0,0.1,t.size)

analyse(t,y)
plt.show()

解决方法

真阶函数有问题。它们是不可微的，如果跳跃发生在数据点之间，频率或相位的微小变化可能不会改变卡方值，从而使最小化算法错误地停止。一个简单的迭代方法很容易解决这个问题。选择一个可微的函数，并逐步使其更像正方形。

这是我的方法：

from matplotlib import pyplot as plt
import numpy as np
from numpy.random import normal
from scipy.optimize import least_squares

def square_wave(x,f,a,h,phi):
    return h + a * np.sign(np.cos( 2 * np.pi *f * x + phi ) )

def soft_wave( x,phi,s=1 ):
    return h + a * np.tanh( s * np.cos( 2 * np.pi *f * x + phi ) )

def residuals( params,xdata,ydata,s ):
    f,phi = params
    yt = soft_wave( xdata,s=s )
    return yt - ydata


tl = np.linspace( -2,8.3,67 )
sl = square_wave( tl,0.39,3.6,2.1,0.25 )
sl += normal( size=len(tl),scale=0.2 )
swl = soft_wave( tl,0.25,s=10 )

soldict = dict()
swdict=dict()
previous = ( 0.4,4,2,0 )
for mys in range( 2,26,4 ):
    res = least_squares( residuals,x0=previous,args=( tl,sl,mys ) )
    print( res.x )
    previous = res.x
    soldict[mys] = np.append( res.x,mys )
    swdict[mys] = soft_wave( tl,*( soldict[mys] ) )

print (np.linalg.inv( np.dot( np.transpose( res.jac ),res.jac ) ) )

fig = plt.figure( figsize=(10,6))
ax= fig.add_subplot( 1,1,1 )
ax.plot( tl,ls='',marker='+' )
ax.plot( tl,swl )
for mys in range( 2,4 ):
    ax.plot( tl,swdict[mys] )
plt.show( )
curve-fittinglmfitnumpypython