In [ ]:
import sys
sys.path.append('../..')

import pyotc

Generate a PSD

pure Lorentzian


In [ ]:
from pyotc.psd import low_pass_filter
from pyotc.psd import gen_filtered_data
from scipy import arange

fs = 1024  # Hz
T_msr = 10

t = arange(0, T_msr, 1/fs)

# this should generate pure lorentzian psds
x = gen_filtered_data(low_pass_filter, fs, T_msr, 1400) * 0.4
y = gen_filtered_data(low_pass_filter, fs, T_msr, 1200) * 0.4
z = gen_filtered_data(low_pass_filter, fs, T_msr, 150) + 4.5

In [ ]:
pyotc.plt.close('all')
fig =pyotc.plt.figure(figsize=(8, 4))
ax = pyotc.add_plot_to_figure(fig, t, x, fmt='-b', linewidth=1, alpha=0.5, label='x')
pyotc.add_plot_to_figure(fig, t, y, fmt='-g', linewidth=1, alpha=0.5, label='y')
pyotc.add_plot_to_figure(fig, t, z, fmt='-', color='orange', linewidth=1, alpha=0.7, label='z',
                         xlabel='Time', ylabel='Signal')

ax.set_xticklabels([])
ax.set_yticklabels([])
fig

Generate an aliased and lp-filtered Lorentzian and hydro PSD


In [ ]:
# create fine data with ~2.5 MHz

fs = 40000  # Hz
fs_full = fs * 25

T_msr = 10

#t = arange(0, T_msr, 1/fs)

# generate a pure lorentzian psd
x_full = gen_filtered_data(low_pass_filter, fs_full, T_msr, 1000)
# aliased sampling
x_aliased = x_full[0::25]

In [ ]:
def lp2_(freq, f3dB_1, f3dB_2, alpha_1, alpha_2):
    f = low_pass_filter(freq, f3dB_1, alpha=alpha_1) * low_pass_filter(freq, f3dB_2, alpha_2)
    return f

y = gen_filtered_data(lp2_, fs, T_msr, 1000, 5000, 0, 0.5)

In [ ]:
from pyotc import k_B, drag

radius = 1.0e-6
temp = 300

D = k_B * temp / drag(radius, temp)
f_c = 1200

height = 10e-6
rho=4500

expset = pyotc.psd.ExpSetting(temp, radius, height=height, material='Titania',
                              temp_unit='K', radius_unit='m', height_unit='m')

psdm = pyotc.PSDMeasurement(exp_setting=expset)

psd = pyotc.gen_PSD_from_time_series(x_aliased, fs, 40, name='x')
psdm.add_psd('x', psd)
psd = pyotc.gen_PSD_from_time_series(y, fs, 40, name='y',)
psdm.add_psd('y', psd)

In [ ]:
from pyotc.psd import hydro_psd, lorentzian_psd

z = gen_filtered_data(hydro_psd, fs, T_msr, D, f_c, radius=radius, height=height, temp=temp, rho=rho) * 1e9  # in nm

psd = pyotc.gen_PSD_from_time_series(z, fs, 40, name='z')
psdm.add_psd('z', psd)

In [ ]:
psdm.plot_psds(figsize=(4.5, 3))

setup fits and fit the data


In [ ]:
pf = pyotc.PSDFit(psdm)
pf.setup_fit(names='x', model='lorentzian', aliasing=True, f_sample=psdm.get_f_sample('x'), N_alias=9)
pf.setup_fit(names='y', model='lorentzian', lp_filter=True, lp_fixed=False, f3dB=10e3, alpha=0.1)
pf.setup_fit(names='z', model='hydro')
#pf.fit_kwargs

In [ ]:
pf.fit_psds()

In [ ]:
pf.plot_fits(plot_data=True, showLegend=False)

fit to user defined function


In [ ]:
from scipy import ones

pf.setup_fit('x', model='other')

kws = pf.fit_kwargs['x']

def fun(freq, D, f_c, parameterx, **other):
    return parameterx * (freq)

# pfun = partial(fun, **other)

kws['model_fun'] = fun
kws['name'] = 'userdefined'
kws['parameterx'] = 10.0
kws['param_names'] = ['D', 'f_c', 'parameterx']
#kws['expon'] = 2
pf.fit_kwargs['x'].update(kws)
pf.fit_kwargs['x']

In [ ]:
pf.fit_psd('x')

In [ ]:
pf.plot_fits(names='x', plot_data=1)

In [ ]: