For NL-SIM we need to have a good way to measure, directly, the NL modulation amplitude. This task is very model dependent. The best way would be to use the fourier transform, but unfortunately the signal to noise and number of points is often too low to make the measurement. For instance if we wanted to nyquist sample a signal with a maximum frequency content of $k_{max}$ we would need a sampling inteval of $1/2 k_{max}$. The period of our pattern is such that we sample $2 \pi$ in $n$ points. So basically our frequency, $f$, is 1. To perfectly Nyquist sample our signal (considering a single harmonic) we need three points across our interval of $[0, 2\pi)$ as that gives the required sampling frequency of $\pi$. Each higher harmonic requires 2 more points. So for 15 points we can nyquist sample up to 7 harmonics. The number of points needed for $H$ harmonics is $2H + 1$

But we still have signal to noise concerns. The other option is to non-linearly fit the signal to a sum of sinusoids, but then you have to fit the phase, amplitude of each sinusoid plus an offset and frequency, so the number of points needed is $2H + 2$ to ensure that the system is not underdetermined.

The final option is to use a variant of LPSVD which would independently fit the amplitude, frequency, phase and offset of each sinusoid. In this case we'd need $4H$ points.

```
In [ ]:
```def cosine(x, amp, f, p, o):
'''
Utility function to fit nonlinearly
'''
return amp*np.cos(2*np.pi*f*(x-p))+o
def cosine_sat(x,*params):
o = params[-1]
p = params[-2]
f = params[-3]
amps = params[:len(params)-3]
to_return = zeros_like(x)
for i, amp in enumerate(amps):
#each amp corresponds to the ith harmonic
to_return += cosine(x,amp,f*(i+1),p,0)
return to_return+o
def calc_mod_sat(data, num_harms = 1, nphases = 24,periods = 2):
'''
Need to change this so that it:
- first tries to fit only the amplitude and phase
- if that doesn't work, estimate amp and only fit phase
- then do full fit
'''
#pull internal number of phases
#nphases = self.nphases
#only deal with finite data
#NOTE: could use masked wave here.
finite_args = np.isfinite(data)
data_fixed = data[finite_args]
popt = None
if len(data_fixed) > 4:
#we can't fit data with less than 4 points
#make x-wave
x = np.arange(nphases,dtype=data_fixed.dtype)[finite_args]
#make guesses
#amp of sine wave is sqrt(2) the standard deviation
g_a = np.sqrt(2)*(data_fixed.std())
#offset is mean
g_o = data_fixed.mean()
#frequency is such that `nphases` covers `periods`
g_f = periods/nphases
#guess of phase is from first data point (maybe mean of all?)
g_p = nan
i = 0
while not isfinite(g_p):
g_p = x[i]-np.arccos((data_fixed[i]-g_o)/g_a)/(2*np.pi*g_f)
i+=1
#make guess sequence
if num_harms == 1:
amps = [g_a]
if num_harms == 2:
# https://www.wolframalpha.com/input/?i=expand+1-cos(x)%5E4
amps = [-4/5*g_a, -1/5*g_a]
g_p+=3*pi/2
if num_harms == 3:
# https://www.wolframalpha.com/input/?i=expand+1-cos(x)%5E6
amps = [-15/22*g_a, -6/22*g_a,-1/22*g_a]
g_p+=3*pi/2
if num_harms == 4:
amps = [-56/93*g_a, -28/93*g_a,-8/93*g_a,-1/93*g_a]
g_p+=3*pi/2
pguess = amps + [g_f,g_p,g_o]
try:
popt,pcov = curve_fit(cosine_sat,x,data_fixed,p0=array(pguess))
except RuntimeError as e:
#if fit fails, put nan
print(e)
mod = np.nan
res = nan
except TypeError as e:
print(e)
print(data_fixed)
mod = np.nan
res = nan
else:
if len(popt)==4:
opt_a,opt_f,opt_p,opt_o = popt
opt_a = np.abs(opt_a)
#if any part of the fit is negative, mark as failure
if opt_o - opt_a < 0:
mod = np.nan
else:
#calc mod
mod = 2*opt_a/(opt_o+opt_a)
else:
mod = np.nan
res = (data_fixed-cosine_sat(x,*popt))**2
res = res.sum()
else:
mod = np.nan
res = nan
return popt, res, mod

```
In [ ]:
```def cosine(x, amp, f, p, o):
'''
Utility function to fit nonlinearly
'''
return amp*np.cos(2*np.pi*f*(x-p))+o
def cosine_sat(x,*params):
o = params[-1]
p = params[-2]
f = params[-3]
amps = params[:len(params)-3]
to_return = zeros_like(x)
for i, amp in enumerate(amps):
#each amp corresponds to the ith harmonic
to_return += cosine(x,amp,f*(i+1),p,0)
return to_return+o
#Testing the function
x = linspace(0,23,1024)
#popt = array([ -9.23688107e+02, -2.15499865e+01, 1/12, 12*pi, 2.31774918e+03])
popt = array([ -56/128, -28/128, -8/128, -1/128, 1/12, 0, 93/128])
fig, ax = subplots(1,1,figsize=(12,12))
fit = cosine_sat(x,*popt)
ax.plot(x,fit,'r-', label = 'Cosine Sat Func')
amps = popt[:len(popt)-3]
# plot each harmonic seperately
for i, amp in enumerate(amps):
#each amp corresponds to the ith harmonic
o = popt[-1]
p = popt[-2]
f = popt[-3]
mod = 2*abs(amp)/(o+abs(amp))
ax.plot(x,cosine(x,amp,f*(i+1),p,-o*(amp/abs(amps).sum())),'--',label='Mod = {:.2f}\nAmp = {:.2f}'.format(mod, amp))
fit_direct = 1 - (cosine(x, 1/2, f, p , 1/2))**4
ax.plot(x, fit_direct,'k.', label= '$\cos(x)^4$')
ax.legend()
def calc_mod_sat(data, num_harms = 1, nphases = 30,periods = 2):
'''
Need to change this so that it:
- first tries to fit only the amplitude and phase
- if that doesn't work, estimate amp and only fit phase
- then do full fit
'''
#pull internal number of phases
#nphases = self.nphases
#only deal with finite data
#NOTE: could use masked wave here.
finite_args = np.isfinite(data)
data_fixed = data[finite_args]
popt = ones(4)*nan
pguess = ones(4)*nan
res = np.nan
if len(data_fixed) > 3+num_harms:
#we can't fit data if number of parameters exceeds number of points
#make x-wave
x = np.arange(nphases,dtype=data_fixed.dtype)[finite_args]
#make guess sequence
if num_harms == 1:
#make guesses
#amp of sine wave is sqrt(2) the standard deviation
g_a = np.sqrt(2)*(data_fixed.std())
#offset is mean
g_o = data_fixed.mean()
#frequency is such that `nphases` covers `periods`
g_f = periods/nphases
#guess of phase is from median, note that amp is made negative due to the saturation model.
g_p = median((x-np.arccos(-(data_fixed-g_o)/g_a)/(2*np.pi*g_f))[:nphases//periods])
amps = [1]
if num_harms == 2:
amps = [4/5, 1/5]
if num_harms == 3:
amps = [15/22, 6/22,1/22]
if num_harms == 4:
amps = [56/93, 28/93,8/93,1/93]
if num_harms > 1:
#first do the fit for num_harms 1
popt, res, pguess = calc_mod_sat(data, num_harms = 1, nphases = nphases,periods = periods)
g_a, g_f, g_p, g_o = popt
#if the fitted amp is positive that just means we need to adjust the phase
# and make the amp negative
if g_a > 0:
g_p += nphases/periods/2
g_a = -g_a
# for saturated data we expect the peaks to be negative, so shift phase and make amplitudes negative.
pguess = concatenate((array(amps)*g_a, [g_f, g_p, g_o]))
if isfinite(pguess.all()):
try:
popt,pcov = curve_fit(cosine_sat,x,data_fixed,p0=pguess)
except RuntimeError as e:
#if fit fails, put nan
pass
except TypeError as e:
print(e)
print(data_fixed)
else:
res = ((data_fixed-cosine_sat(x, *popt))**2).sum()
return popt, res, pguess
def find_best_num_harms(y, max_harm = 4, verbose =False):
'''
A function to find the best number of harmonics to fit the data
'''
# set up the stopping criterion
amps_valid = False
mods_valid = False
res_old = np.inf
popt_old = np.nan
# loop through having 2 to 4 harmonics, fitting each one
for i in list(range(2, max_harm+1))+[1]:
popt, res, pguess = calc_mod_sat(y, i)
if verbose:
print('pguess', pguess)
print('popt', popt)
# pull the offset
o = popt[-1]
# and the amps
amps = popt[:len(popt)-3]
# check to see if amps are positive
amps_pos = (amps < 0).all()
# check to see that the amplitudes are ordered such that the strength decreasing with increasing harmonics
amps_valid = amps_pos and (amps.argsort() == arange(len(amps))).all()
# check the modulation depths as well.
mods = 2*abs(amps)/(o+abs(amps))
#make sure modulation depths are between 1 and 0, inclusive
mods_valid = np.logical_and(mods <= 1.0, mods >= 0.0).all()
if res > res_old:
# if the residuals have increased, don't update
if verbose:
print('Residuals have increased')
# break
elif mods_valid and (amps_valid or len(popt) == 4):
# if the residuals have decreased and the parameter are valid, update
if verbose:
print('Residuals have decreased and popt valid, updating...')
res_old = res
popt_old = popt
else:
if verbose:
print('Residuals have decreased but popt invalid')
if verbose:
print('{:.3e}'.format(res_old))
if verbose:
print('Final residuals = {:.3e}'.format(res_old))
return popt_old, res_old, pguess

Probably the best thing to do is to have a function which calculates the saturation directly from a relevant model then you only have to fit:

- amplitude
- phase
- frequency
- saturation factor
- offset

For any saturation factor or number of harmonics

The following model functions should be used: power law ($a x^{-b}$), exponential function ($a e^{-b x}$) and saturation ($a/(1 + x/b)$)

```
In [ ]:
```