Серия Фурье Fit в Python

Question

У меня есть некоторые данные, которые я хочу уместить, используя ряды Фурье 2-й, 3-й или 4-й степени.

В то время как этот вопрос и ответ о переполнении стека приближается к тому, что я хочучтобы использовать scipy, они уже заранее определяют свои коэффициенты как tau = 0.045 всегда.Я хочу, чтобы моя подборка нашла возможные коэффициенты (a0, w1, w2, w3 и т. Д.) С 95% -ным доверительным интервалом, точно так же, как MATLAB fit fit , эквивалентный для ряда Фурье.Другой вариант, который я видел, использовал fourier_series от sympy , однако эта функция работает только с символическими параметрами, соответствующими определенной функции, а не необработанными данными.

1) Есть ли способ для симпози fourier_series получать необработанные данные, а не функцию или другой способ использования этой библиотеки?

2) Или подгонка кривой для данныхчто есть несколько неизвестных (коэффициенты)

Answer 1

Я думаю, что, вероятно, будет проще использовать этот API для вызова функций MATLAB в скриптах Python для выполнения всей тяжелой работы вместо того, чтобы разбираться со всеми деталями sympy и scipy.

Мое решение заключалось в следующем:

1. установить набор инструментов функций Matlab и примерки
2. Установить ядро ​​Matlab для python, запустив python setup.py install в корневой папке matlab в extern / двигателей /питон.ПРИМЕЧАНИЕ: он работает только для Python 2.7, 3.5 и 3.6

3. Используйте следующий код:

   import matlab.engine
   import numpy as np
   
   eng = matlab.engine.start_matlab()
   eng.evalc("idx = transpose(1:{})".format(len(diff)))
   eng.workspace['y'] = eng.transpose(eng.cell2mat(diff.tolist()))
   eng.evalc("f = fit(idx, y, 'fourier3')")
   y_f = eng.evalc("f(idx)").replace('ans =', '')
   
   y_f = np.fromstring(y_f, dtype=float, sep='\n')
   

Несколько замечаний:

1. eng.workspace['myVariable'] - объявлять переменные Matlab с использованием результатов Python, которые впоследствии могут быть вызваны с помощью evalc
2. eng.evalc возвращает строку в форме 'ans = ...'
3. diff в этом коде была просто разница между некоторыми данными и их линией наименьших квадратов, и она имела тип серии.
4. Списки Python эквивалентны типу ячейки в MATLAB.

Answer 2

Вы можете остаться очень близко к sympy коду для подбора данных, если хотите, используя пакет, который я написал для этой цели и который называется symfit.Он в основном обёртывает scipy с использованием интерфейса sympy.Используя symfit, вы можете сделать что-то вроде следующего:

from symfit import parameters, variables, sin, cos, Fit
import numpy as np
import matplotlib.pyplot as plt

def fourier_series(x, f, n=0):
    """
    Returns a symbolic fourier series of order `n`.

    :param n: Order of the fourier series.
    :param x: Independent variable
    :param f: Frequency of the fourier series
    """
    # Make the parameter objects for all the terms
    a0, *cos_a = parameters(','.join(['a{}'.format(i) for i in range(0, n + 1)]))
    sin_b = parameters(','.join(['b{}'.format(i) for i in range(1, n + 1)]))
    # Construct the series
    series = a0 + sum(ai * cos(i * f * x) + bi * sin(i * f * x)
                     for i, (ai, bi) in enumerate(zip(cos_a, sin_b), start=1))
    return series

x, y = variables('x, y')
w, = parameters('w')
model_dict = {y: fourier_series(x, f=w, n=3)}
print(model_dict)

Это напечатает желаемую символическую модель:

{y: a0 + a1*cos(w*x) + a2*cos(2*w*x) + a3*cos(3*w*x) + b1*sin(w*x) + b2*sin(2*w*x) + b3*sin(3*w*x)}

Далее я подгоню это к простому шагуфункция, чтобы показать вам, как это работает:

# Make step function data
xdata = np.linspace(-np.pi, np.pi)
ydata = np.zeros_like(xdata)
ydata[xdata > 0] = 1
# Define a Fit object for this model and data
fit = Fit(model_dict, x=xdata, y=ydata)
fit_result = fit.execute()
print(fit_result)

# Plot the result
plt.plot(xdata, ydata)
plt.plot(xdata, fit.model(x=xdata, **fit_result.params).y, color='green', ls=':')

Это напечатает:

Parameter Value        Standard Deviation
a0        5.000000e-01 2.075395e-02
a1        -4.903805e-12 3.277426e-02
a2        5.325068e-12 3.197889e-02
a3        -4.857033e-12 3.080979e-02
b1        6.267589e-01 2.546980e-02
b2        1.986491e-02 2.637273e-02
b3        1.846406e-01 2.725019e-02
w         8.671471e-01 3.132108e-02
Fitting status message: Optimization terminated successfully.
Number of iterations:   44
Regression Coefficient: 0.9401712713086535

И выдает следующий график:

Это так просто!Остальное оставляю на ваше воображение.Для получения дополнительной информации вы можете найти документацию здесь .

Answer 3

Это простой код, который я недавно использовал.Я знал частоту, поэтому она немного изменена, чтобы соответствовать ей.Начальное предположение должно быть довольно хорошим, хотя (я думаю).С другой стороны, есть несколько хороших способов получить базовую частоту.

import numpy as np
from scipy.optimize import leastsq


def make_sine_graph( params, xData):
    """
    take amplitudes A and phases P in form [ A0, A1, A2, ..., An, P0, P1,..., Pn ]
    and construct function f = A0 sin( w t + P0) + A1 sin( 2 w t + Pn ) + ... + An sin( n w t + Pn )
    and return f( x )
    """
    fr = params[0]
    npara = params[1:]
    lp =len( npara )
    amps = npara[ : lp // 2 ]
    phases = npara[ lp // 2 : ]
    fact = range(1, lp // 2 + 1 )
    return [ sum( [ a * np.sin( 2 * np.pi * x * f * fr + p ) for a, p, f in zip( amps, phases, fact ) ] ) for x in xData ]

def sine_residuals( params , xData, yData):
    yTh = make_sine_graph( params, xData )
    diff = [ y -  yt for y, yt in zip( yData, yTh ) ]
    return diff

def sine_fit_graph( xData, yData, freqGuess=100., sineorder = 3 ):
    aStart = sineorder * [ 0 ]
    aStart[0] = max( yData )
    pStart = sineorder * [ 0 ]
    result, _ = leastsq( sine_residuals, [ freqGuess ] + aStart + pStart, args=( xData, yData ) )
    return result


if __name__ == '__main__':
    import matplotlib.pyplot as plt
    timeList = np.linspace( 0, .1, 777 )
    signalList = make_sine_graph( [  113.7 ] + [1,.5,-.3,0,.1, 0,.01,.02,.03,.04], timeList )

    result = sine_fit_graph( timeList, signalList, freqGuess=110., sineorder = 3 )
    print result
    fitList =  make_sine_graph( result, timeList )
    fig = plt.figure()
    ax = fig.add_subplot( 1, 1 ,1 )
    ax.plot( timeList, signalList )
    ax.plot( timeList, fitList, '--' )
    plt.show()

Обеспечение

<< [ 1.13699742e+02  9.99722859e-01 -5.00511764e-01  3.00772260e-01
     1.04248878e-03 -3.13050074e+00 -3.12358208e+00 ]