search

tBuLi / symfit

Symbolic Fitting; fitting as it should be.
http://symfit.readthedocs.org
MIT License
233 stars 17 forks source link

Fourier series fit example #216

Closed virsto closed 5 years ago

virsto commented 5 years ago

I do not understand the use of a0, *cos_a = parameters( . . .) since a0 is defined in the string generated. Would you please explain why this is neccessary.

tBuLi commented 5 years ago

See the docs here for more info on the function parameters.

If your issue is with the left hand side of the expression, in python 3 you can tuple-unpack a sequence and allow all the remainder to be assigned to a single variable. In this case, parameters returns a list of Parameter objects. The first one is assigned to a0, the rest to cos_a.

Let me know if there are more thing unclear.

p.s. this line is just a shorthand for doing 

a_0 = Parameter('a_0')
...
a_N = Parameter('a_n')
tBuLi commented 5 years ago

I will close this for now, if you have further questions you are always welcome to ask them, even if the issue is closed I'll keep following this discussion.

virsto commented 5 years ago

I have attached some code that might be useful for those who are interested in learning more about symfit applications.

Best regards, --V. Stokes

När du har kontakt med oss på Uppsala universitet med e-post så innebär det att vi behandlar dina personuppgifter. För att läsa mer om hur vi gör det kan du läsa här: http://www.uu.se/om-uu/dataskydd-personuppgifter/

E-mailing Uppsala University means that we will process your personal data. For more information on how this is performed, please read here: http://www.uu.se/en/about-uu/data-protection-policy

coding:utf-8 # this allows utf-8 characters

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

For setting size and position of matplotlib figure

import matplotlib matplotlib.use("WXAgg") import sys """ Purpose: Example of using callable symfit Fourier Series model to fit
sinusoidal functions and a step function. Also, shows how to extrapolate (frontend and backend) fits. Note:

  1. FS model is linear in coefficients of sinusoidal functions but nonlinear in base frequency (w)
  2. This code can easily be modified to handlle more test cases and fit types
  3. Code based on example at: https://symfit.readthedocs.io/en/stable/examples/ex_fourier_series.html

  4. Python packages used: python - 3.6.5 (v3.6.5:f59c0932b4, Mar 28 2018, 17:00:18) matplotlib - vers. 3.0.2 numpy - vers. 1.15.2+mk1 symfit - vers. 0.4.6

    Coded by: V. Stokes (vs@it.uu.se) Date: 2019.02.28
    """ def ext_xdata(orgL,dx,extF,extB): """ Purpose: add x-values to frontend and/or to backend of xdata Inputs: orgL - original (unextended) length (number) of x-values dx - delta (spacing) between x-values extF - number of x-values to be added to frontend (>= 0) extB - number of x-values to be added to backend (>= 0) Otputs: extx - extended xdata """

    Simple numpy array comprehension for extended x-values

    extx = np.fromiter([k*dx for k in np.arange(-extF,orgL+extB)],np.float_)
    return extx

def addNoise(noiseSF,ydata): """ Purpose: add scaled noise, sfN(0,1) to signal to be modelled where, N(0,1) is standardized normal distribution (mean=0,sd=1) Inputs: noiseSF - scale factor (sf) for noise ydata - TS values to which noise is added Otputs: ydata - TS with noise added """ np.random.seed(131) # for repeatability (same noise series for each execution) N = len(ydata)
ydata_rms = np.sqrt((1/N)np.sum(ydata2)) noiseTS = noiseSFnp.random.randn(np.size(ydata)) # white Gaussian noise noi_rms = np.sqrt((1/N)np.sum(noiseTS2)) ydata += noiseTS
print(" Noise = %6.2fN(0,1)"%noiseSF) print(" SNR(dB) = %8.3f"%(20np.log10(ydata_rms/noi_rms)))
return ydata

def genStepData(): """ Purpose: generate x,y-values for step function centered about x = 0 """
xdata = np.linspace(-np.pi, np.pi) ydata = np.zeros_like(xdata) ydata[xdata > 0] = 20 # step transition from y = 0 to 20 at x = 0
return xdata,ydata

def genSinusoidalData(N,Np,NSigs): """ Purpose: generate x,y-values for sinusoidal functions Inputs: N - number of samples in signal Np - number of samples to be used for fit (Np <= N)
NSigs - number of sinusoidal functions added to generate composite signal Otputs: n - order of FS model to used for fit arg - 2np.pit (xdata) (length = Np) sig - composite signal (length = Np) to be approximated by FS """

Center frequencies of signal components

sigs  = NSigs*[0]
if NSigs == 1:
    fc = [ 3]
    am = [ 5]
    n = 1
elif NSigs == 2:
    fc = [ 1,  3]
    am = [20,  5]
    n = 3 
elif NSigs == 3:    
    fc = [ 1,  2, 3]
    am = [20, 10, 5]
    n = 4            

dt   = 1/N      # time between samples (sec)
fs   = 1/dt     # sample frequency     (Hz)
fny  = fs/2     # Nyquist frequency    (Hz)

t  = dt*np.arange(0.0,Np)    # Np samples to be used for fit
# print('t = ',t)

print (" dt   (time between samples) = %10.6f"%dt)
print (' fs   (Sample freq)          = %6.1f [Hz]'%fs)
print (' fny  (Nyquist freq)         = %6.1f [Hz]'%fny)

FMax = max(fc)  # max freq present in in signals
print (' Fmax (max freq in signal)   = %4d [Hz]'%FMax)
"""
if FMax >= fny:
    print ('!WARNING -- sample frequency (fs) must be > Fmax for PSD estimation') 
    sys.exit()   
"""
arg   = 2*np.pi*t   # (xdata)
if NSigs   == 1:
    sigs[0]  = am[0]*np.sin(arg*fc[0])
elif NSigs == 2:
    sigs[0]  = am[0]*np.sin(arg*fc[0])
    sigs[1]  = am[1]*np.cos(arg*fc[1])  
elif NSigs == 3:
    sigs[0]  = am[0]*np.cos(arg*fc[0])
    sigs[1]  = am[1]*np.sin(arg*fc[1])
    sigs[2]  = am[2]*np.cos(arg*fc[2])

sig     = sum(sigs)      # composite (sum) of sinusoidal functions
return n,arg,sig

def fourier_series(x, f, n=1): """ Purpose: generate symbolic Fourier series model of order-n Inputs: x - independent variable f - base (fundamental) frequency of the Fourier series (radians) n - order of the Fourier series. Otputs: series - symbolic equation for FS of order-n """

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 FS
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

def plotResults(TESTCASE,xdata,ydata,extx,y1_eval,y2_eval,title): """ Purpose: plot the results of Fourier series fit and any extensions to fit Inputs: TESTCASE - number of test case xdata - x-values for original (unextended) data ydata - y-values for original (unextended) data extx - extended x-values y1_eval - y-values for fit to unextended data y2eval - y-values for extended data from unextended fit title - title of plot figure """ fig, = plt.subplots(1,1,constrained_layout=True,figsize=(9.0,5.0)) fig.canvas.set_window_title('SymFit - Fourier Series Approximation, TESTCASE = %d'%TESTCASE)

Used to set position of the plot

pltManager = plt.get_current_fig_manager()
pltManager.window.SetPosition((10,10))  # pixels offset from top left corner of display     

#p1, = plt.plot(xdata, ydata,lw=2.5,label='Data for FS fit')   # data to be fitted with a FS
p1, = plt.plot(xdata, ydata,marker='.',ls=':',label='Data for FS fit')   # data to be fitted with a FS
p2, = plt.plot(extx, y2_eval,ls=':',color='red',label='FS fit extended')
p3, = plt.plot(xdata,y1_eval,lw=2.0,color='yellow',label='FS fit')

plt.xlabel('x')
plt.ylabel('y')
plt.legend(handles=[p1,p3,p2])
plt.grid(True)
plt.title(title)
plt.show()

def indexChk(N,Np): print (" N,Np: %4d,%4d"%(N,Np)) if Np > N: print ("ERROR - Np > N not allowed!") sys.exit()
def chkRange(name,value,Range): print (' %s: %d'%(name,value) ) if value not in Range: print("ERROR - %s not in range!"%name) sys.exit()

--------------------------------------------------------------------------------

if name == 'main':

Set maximum value for TESTCASE

testCaseMx = 2

TESTCASE = 1   # sinusoidal functions
#TESTCASE = 2   # step function

chkRange('TESTCASE',TESTCASE,[1,testCaseMx])
if TESTCASE == 1:
    # Set maximum value for nsigMx (maximum number of sinusoidal functions that can be added)
    nsigMx = 3

# Set noise level to be added to signal (or no noise)
noise_sf = 0.0    # no noise
noise_sf = 1.5

if TESTCASE == 1:
    N      = 64   # number of samples in input signal
    # Select NP, the number of samples in input signal to be used for fit
    Np     = 56   # use Np samples of N for FS fit of input signal
    indexChk(N,Np)  

    # Set the number of sinusoidal functions for input signal
    nsigs = 1
    nsigs = 2
    nsigs = 3
    chkRange('nsigs',nsigs,[1,nsigMx])
    # Generate Sinusoidal data        
    n,xdata,ydata = genSinusoidalData(N,Np,nsigs) 
    lenx  = len(xdata)
    # Set Extensions (>= 0)
    extF  = lenx   # extension to frontend
    extB  = lenx   # extension to backend
    title = 'Fourier series (order=%2d) fit to Sinusoids'%n

elif TESTCASE == 2:
    # Generate step function data
    xdata,ydata = genStepData()
    N  = len(xdata)  
    # Select NP, the number of samples in input data to be used for fit
    Np = N   # all samples in input data used for fit
    indexChk(N,Np)

    # Set FS order
    n  = 3
    # Set Extensions (>= 0)
    extF  = 100
    extB  = N+10 
    title = 'Fourier series (order=%2d) fit to Step function'%n 
else:
    # This should never be executed!
    sys.exit()

if noise_sf > 0.0:       
    ydata = addNoise(noise_sf,ydata) 
    title = title + ' + %3.1f*N(0,1)'%noise_sf 
dx = xdata[1] - xdata[0]
print (" dx   (delta between x)      = %10.6f"%dx)

# Generate symbolic Fourier series (FS)
x, y = variables('x, y')
w,   = parameters('w')     # model is nonlinear in this parameter  (base frequency)    
Fourier_model = CallableModel({y: fourier_series(x, f=w, n=n)})
print ('\nFourier Series Model:')
print (Fourier_model)

# Define the Fit object for this model with data ...
fit        = Fit(Fourier_model, x=xdata, y=ydata)
# ... and fit it
fit_result = fit.execute()    
print(fit_result)
# Another way to extract Fourier coefficients
[print("%3s  %20.12f" %(k,v)) for k, v in fit_result.params.items()]    

# Extend xdata    
extx = ext_xdata(Np,dx,extF,extB)   

print (" x for FS fit          in [%10.6f, %10.6f]"%(xdata[0],xdata[-1]))
print (' x for FS fit extended in [%10.6f, %10.6f]'%(extx[0],  extx[-1]))    

y1_eval, = Fourier_model(x=xdata,**fit_result.params)  # evaluate FS model over xdata    
y2_eval, = Fourier_model(x=extx, **fit_result.params)  # evaluate FS model over extx 

# Plot the results
plotResults(TESTCASE,xdata,ydata,extx,y1_eval,y2_eval,title)
pckroon commented 5 years ago

Nice one, thanks. Small question: why do you use a CallableModel instead of a normal Model?

tBuLi commented 5 years ago

Thanks for sharing, it's always nice to have more code examples out there.

Apparently I cannot add markdown to an email reply, so I'll repost the code here with code highlighting:

# coding:utf-8    # this allows utf-8 characters
from   symfit import parameters, variables, sin, cos, Fit, CallableModel
import numpy as np
import matplotlib.pyplot as plt
# For setting size and position of matplotlib figure
import matplotlib
matplotlib.use("WXAgg")
import sys
"""
  Purpose: Example of using callable symfit Fourier Series model to fit
           sinusoidal functions and a step function. Also, shows how to
           extrapolate (frontend and backend) fits.
  Note:
   1. FS model is linear in coefficients of sinusoidal functions but nonlinear in base
      frequency (w)
   2. This code can easily be modified to handlle more test cases and fit types
   3. Code based on example at:
      https://symfit.readthedocs.io/en/stable/examples/ex_fourier_series.html
   4. Python packages used:
      python     - 3.6.5 (v3.6.5:f59c0932b4, Mar 28 2018, 17:00:18)
      matplotlib - vers. 3.0.2
      numpy      - vers. 1.15.2+mk1
      symfit     - vers. 0.4.6

  Coded by: V. Stokes (vs@it.uu.se)
  Date: 2019.02.28
"""
def ext_xdata(orgL,dx,extF,extB):
    """
     Purpose: add x-values to frontend and/or to backend of xdata
      Inputs:
        orgL - original (unextended) length (number) of x-values
        dx   - delta (spacing) between x-values
        extF - number of x-values to be added to frontend  (>= 0)
        extB - number of x-values to be added to backend   (>= 0)
      Otputs:
        extx  - extended xdata
    """
    # Simple numpy array comprehension for extended x-values
    extx = np.fromiter([k*dx for k in np.arange(-extF,orgL+extB)],np.float_)
    return extx

def addNoise(noiseSF,ydata):
    """
     Purpose: add scaled noise, sf*N(0,1) to signal to be modelled
              where, N(0,1) is standardized normal distribution (mean=0,sd=1)
      Inputs:
       noiseSF - scale factor (sf) for noise
       ydata   - TS values to which noise is added
      Otputs:
       ydata   - TS with noise added
    """
    np.random.seed(131)  # for repeatability (same noise series for each execution)
    N         = len(ydata)
    ydata_rms = np.sqrt((1/N)*np.sum(ydata**2))
    noiseTS   = noiseSF*np.random.randn(np.size(ydata))  # white Gaussian noise
    noi_rms   = np.sqrt((1/N)*np.sum(noiseTS**2))
    ydata    +=  noiseTS
    print(" Noise                       = %6.2f*N(0,1)"%noiseSF)
    print(" SNR(dB)                     = %8.3f"%(20*np.log10(ydata_rms/noi_rms)))
    return  ydata

def genStepData():
    """
     Purpose: generate x,y-values for step function centered about x = 0
    """
    xdata = np.linspace(-np.pi, np.pi)
    ydata = np.zeros_like(xdata)
    ydata[xdata > 0] = 20    # step transition from y = 0 to 20 at x = 0
    return xdata,ydata

def genSinusoidalData(N,Np,NSigs):
    """
    Purpose: generate x,y-values for sinusoidal functions
     Inputs:
      N     - number of samples in signal
      Np    - number of samples to be used for fit  (Np <= N)
      NSigs - number of sinusoidal functions added to generate composite signal
     Otputs:
      n     - order of FS model to used for fit
      arg   - 2*np.pi*t  (xdata) (length = Np)
      sig   - composite signal   (length = Np) to be approximated by FS
    """
    # Center frequencies of signal components
    sigs  = NSigs*[0]
    if NSigs == 1:
        fc = [ 3]
        am = [ 5]
        n = 1
    elif NSigs == 2:
        fc = [ 1,  3]
        am = [20,  5]
        n = 3
    elif NSigs == 3:
        fc = [ 1,  2, 3]
        am = [20, 10, 5]
        n = 4

    dt   = 1/N      # time between samples (sec)
    fs   = 1/dt     # sample frequency     (Hz)
    fny  = fs/2     # Nyquist frequency    (Hz)

    t  = dt*np.arange(0.0,Np)    # Np samples to be used for fit
    # print('t = ',t)

    print (" dt   (time between samples) = %10.6f"%dt)
    print (' fs   (Sample freq)          = %6.1f [Hz]'%fs)
    print (' fny  (Nyquist freq)         = %6.1f [Hz]'%fny)

    FMax = max(fc)  # max freq present in in signals
    print (' Fmax (max freq in signal)   = %4d [Hz]'%FMax)
    """
    if FMax >= fny:
        print ('!WARNING -- sample frequency (fs) must be > Fmax for PSD estimation')
        sys.exit()
    """
    arg   = 2*np.pi*t   # (xdata)
    if NSigs   == 1:
        sigs[0]  = am[0]*np.sin(arg*fc[0])
    elif NSigs == 2:
        sigs[0]  = am[0]*np.sin(arg*fc[0])
        sigs[1]  = am[1]*np.cos(arg*fc[1])
    elif NSigs == 3:
        sigs[0]  = am[0]*np.cos(arg*fc[0])
        sigs[1]  = am[1]*np.sin(arg*fc[1])
        sigs[2]  = am[2]*np.cos(arg*fc[2])

    sig     = sum(sigs)      # composite (sum) of sinusoidal functions
    return n,arg,sig

def fourier_series(x, f, n=1):
    """
    Purpose: generate symbolic Fourier series model of order-n
     Inputs:
      x      - independent variable
      f      - base (fundamental) frequency of the Fourier series  (radians)
      n      - order of the Fourier series.
     Otputs:
      series - symbolic equation for FS of order-n
    """
    # 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 FS
    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

def plotResults(TESTCASE,xdata,ydata,extx,y1_eval,y2_eval,title):
    """
     Purpose: plot the results of Fourier series fit and any extensions to fit
      Inputs:
       TESTCASE - number of test case
       xdata    - x-values for original (unextended) data
       ydata    - y-values for original (unextended) data
       extx     - extended x-values
       y1_eval  - y-values for fit to unextended data
       y2_eval  - y-values for extended data from unextended fit
       title    - title of plot figure
    """
    fig,_ = plt.subplots(1,1,constrained_layout=True,figsize=(9.0,5.0))
    fig.canvas.set_window_title('SymFit - Fourier Series Approximation, TESTCASE = %d'%TESTCASE)
    # Used to set position of the plot
    pltManager = plt.get_current_fig_manager()
    pltManager.window.SetPosition((10,10))  # pixels offset from top left corner of display

    #p1, = plt.plot(xdata, ydata,lw=2.5,label='Data for FS fit')   # data to be fitted with a FS
    p1, = plt.plot(xdata, ydata,marker='.',ls=':',label='Data for FS fit')   # data to be fitted with a FS
    p2, = plt.plot(extx, y2_eval,ls=':',color='red',label='FS fit extended')
    p3, = plt.plot(xdata,y1_eval,lw=2.0,color='yellow',label='FS fit')

    plt.xlabel('x')
    plt.ylabel('y')
    plt.legend(handles=[p1,p3,p2])
    plt.grid(True)
    plt.title(title)
    plt.show()

def indexChk(N,Np):
    print (" N,Np: %4d,%4d"%(N,Np))
    if Np > N:
        print ("ERROR - Np > N not allowed!")
        sys.exit()
def chkRange(name,value,Range):
    print (' %s: %d'%(name,value) )
    if value not in Range:
        print("ERROR - %s not in range!"%name)
        sys.exit()
#--------------------------------------------------------------------------------
if __name__ == '__main__':
    # Set maximum value for TESTCASE
    testCaseMx = 2

    TESTCASE = 1   # sinusoidal functions
    #TESTCASE = 2   # step function

    chkRange('TESTCASE',TESTCASE,[1,testCaseMx])
    if TESTCASE == 1:
        # Set maximum value for nsigMx (maximum number of sinusoidal functions that can be added)
        nsigMx = 3

    # Set noise level to be added to signal (or no noise)
    noise_sf = 0.0    # no noise
    noise_sf = 1.5

    if TESTCASE == 1:
        N      = 64   # number of samples in input signal
        # Select NP, the number of samples in input signal to be used for fit
        Np     = 56   # use Np samples of N for FS fit of input signal
        indexChk(N,Np)

        # Set the number of sinusoidal functions for input signal
        nsigs = 1
        nsigs = 2
        nsigs = 3
        chkRange('nsigs',nsigs,[1,nsigMx])
        # Generate Sinusoidal data
        n,xdata,ydata = genSinusoidalData(N,Np,nsigs)
        lenx  = len(xdata)
        # Set Extensions (>= 0)
        extF  = lenx   # extension to frontend
        extB  = lenx   # extension to backend
        title = 'Fourier series (order=%2d) fit to Sinusoids'%n

    elif TESTCASE == 2:
        # Generate step function data
        xdata,ydata = genStepData()
        N  = len(xdata)
        # Select NP, the number of samples in input data to be used for fit
        Np = N   # all samples in input data used for fit
        indexChk(N,Np)

        # Set FS order
        n  = 3
        # Set Extensions (>= 0)
        extF  = 100
        extB  = N+10
        title = 'Fourier series (order=%2d) fit to Step function'%n
    else:
        # This should never be executed!
        sys.exit()

    if noise_sf > 0.0:
        ydata = addNoise(noise_sf,ydata)
        title = title + ' + %3.1f*N(0,1)'%noise_sf
    dx = xdata[1] - xdata[0]
    print (" dx   (delta between x)      = %10.6f"%dx)

    # Generate symbolic Fourier series (FS)
    x, y = variables('x, y')
    w,   = parameters('w')     # model is nonlinear in this parameter  (base frequency)
    Fourier_model = CallableModel({y: fourier_series(x, f=w, n=n)})
    print ('\nFourier Series Model:')
    print (Fourier_model)

    # Define the Fit object for this model with data ...
    fit        = Fit(Fourier_model, x=xdata, y=ydata)
    # ... and fit it
    fit_result = fit.execute()
    print(fit_result)
    # Another way to extract Fourier coefficients
    [print("%3s  %20.12f" %(k,v)) for k, v in fit_result.params.items()]

    # Extend xdata
    extx = ext_xdata(Np,dx,extF,extB)

    print (" x for FS fit          in [%10.6f, %10.6f]"%(xdata[0],xdata[-1]))
    print (' x for FS fit extended in [%10.6f, %10.6f]'%(extx[0],  extx[-1]))

    y1_eval, = Fourier_model(x=xdata,**fit_result.params)  # evaluate FS model over xdata
    y2_eval, = Fourier_model(x=extx, **fit_result.params)  # evaluate FS model over extx

    # Plot the results
    plotResults(TESTCASE,xdata,ydata,extx,y1_eval,y2_eval,title)