def OutMidCorrection(self, correctionType, firOrd, fs):
"""
Method to "correct" the middle outer ear transfer function.
As appears in :
- A. Härmä, and K. Palomäki, ''HUTear – a free Matlab toolbox for modeling of human hearing'',
in Proceedings of the Matlab DSP Conference, pp 96-99, Espoo, Finland 1999.
"""
# Lookup tables for correction
f1 = np.array([20, 25, 30, 35, 40, 45, 50, 55, 60, 70, 80, 90, 100, \
125, 150, 177, 200, 250, 300, 350, 400, 450, 500, 550, \
600, 700, 800, 900, 1000, 1500, 2000, 2500, 2828, 3000, \
3500, 4000, 4500, 5000, 5500, 6000, 7000, 8000, 9000, 10000, \
12748, 15000])
ELC = np.array([ 31.8, 26.0, 21.7, 18.8, 17.2, 15.4, 14.0, 12.6, 11.6, 10.6, \
9.2, 8.2, 7.7, 6.7, 5.3, 4.6, 3.9, 2.9, 2.7, 2.3, \
2.2, 2.3, 2.5, 2.7, 2.9, 3.4, 3.9, 3.9, 3.9, 2.7, \
0.9, -1.3, -2.5, -3.2, -4.4, -4.1, -2.5, -0.5, 2.0, 5.0, \
10.2, 15.0, 17.0, 15.5, 11.0, 22.0])
MAF = np.array([ 73.4, 65.2, 57.9, 52.7, 48.0, 45.0, 41.9, 39.3, 36.8, 33.0, \
29.7, 27.1, 25.0, 22.0, 18.2, 16.0, 14.0, 11.4, 9.2, 8.0, \
6.9, 6.2, 5.7, 5.1, 5.0, 5.0, 4.4, 4.3, 3.9, 2.7, \
0.9, -1.3, -2.5, -3.2, -4.4, -4.1, -2.5, -0.5, 2.0, 5.0, \
10.2, 15.0, 17.0, 15.5, 11.0, 22.0])
f2 = np.array([ 125, 250, 500, 1000, 1500, 2000, 3000, \
4000, 6000, 8000, 10000, 12000, 14000, 16000])
MAP = np.array([ 30.0, 19.0, 12.0, 9.0, 11.0, 16.0, 16.0, \
14.0, 14.0, 9.9, 24.7, 32.7, 44.1, 63.7])
if correctionType == 'ELC':
freqTable = f1
CorrectionTable = ELC
elif correctionType == 'MAF':
freqTable = f1
CorrectionTable = MAF
elif correctionType == 'MAP':
freqTable = f2
CorrectionTable = MAP
else :
print('Unrecongised operation: ELC will be used instead...')
freqTable = f1
CorrectionTable = ELC
freqN = np.arange(0, firOrd) * fs/2. / (firOrd-1)
spline = uspline(freqTable, CorrectionTable)
crc = spline(freqN)
crclin = 10. ** (-crc/ 10.)
return crclin, freqN, crc
python类InterpolatedUnivariateSpline()的实例源码
def estimate_lorentzian_dip(self, x_axis, data, params):
""" Provides an estimator to obtain initial values for lorentzian function.
@param numpy.array x_axis: 1D axis values
@param numpy.array data: 1D data, should have the same dimension as x_axis.
@param lmfit.Parameters params: object includes parameter dictionary which
can be set
@return tuple (error, params):
Explanation of the return parameter:
int error: error code (0:OK, -1:error)
Parameters object params: set parameters of initial values
"""
# check if parameters make sense
error = self._check_1D_input(x_axis=x_axis, data=data, params=params)
# check if input x-axis is ordered and increasing
sorted_indices = np.argsort(x_axis)
if not np.all(sorted_indices == np.arange(len(x_axis))):
x_axis = x_axis[sorted_indices]
data = data[sorted_indices]
data_smooth, offset = self.find_offset_parameter(x_axis, data)
# data_level = data-offset
data_level = data_smooth - offset
# calculate from the leveled data the amplitude:
amplitude = data_level.min()
smoothing_spline = 1 # must be 1<= smoothing_spline <= 5
function = InterpolatedUnivariateSpline(x_axis, data_level,
k=smoothing_spline)
numerical_integral = function.integral(x_axis[0], x_axis[-1])
x_zero = x_axis[np.argmin(data_smooth)]
# according to the derived formula, calculate sigma. The crucial part is
# here that the offset was estimated correctly, then the area under the
# curve is calculated correctly:
sigma = np.abs(numerical_integral / (np.pi * amplitude))
# auxiliary variables
stepsize = x_axis[1] - x_axis[0]
n_steps = len(x_axis)
params['amplitude'].set(value=amplitude, max=-1e-12)
params['sigma'].set(value=sigma, min=stepsize / 2,
max=(x_axis[-1] - x_axis[0]) * 10)
params['center'].set(value=x_zero, min=(x_axis[0]) - n_steps * stepsize,
max=(x_axis[-1]) + n_steps * stepsize)
params['offset'].set(value=offset)
return error, params
def setup_xarray(scan_list,scan_list2,scan_range,scan_range2,log,log2,precis,precis2,scan_type):
if scan_type=='1d':
#for scanning a predefined list of numbers
if scan_list!=[]:
#interpolate to prescribed length
l1_spl=US(range(len(scan_list)),scan_list)
XARRAY=l1_spl(np.linspace(0,len(scan_list)-1,precis))
X2ARRAY=None
#for having a second simultaneously varying variable
if scan_list2!=[]:
#interpolate to prescribed length
l2_spl=US(scan_list,scan_list2)
X2ARRAY=l2_spl(XARRAY)
#log- or linearly spaced 1d scan
if (scan_list==[] and scan_list2==[]):
X2ARRAY=None #np.empty(1,dtype=np.float)
if log:
XARRAY=np.logspace(np.log10(scan_range[0]),np.log10(scan_range[1]),precis)
else:
XARRAY=np.linspace(scan_range[0],scan_range[1],precis)
elif scan_type=='2d':
#for scanning a predefined list of numbers
if scan_list!=[]:
#interpolate to prescribed length
l1_spl=US(range(len(scan_list)),scan_list)
XARRAY=np.tile(l1_spl(np.linspace(0,len(scan_list)-1,precis)),precis).flatten()
#for having a second simultaneously varying variable
if scan_list2!=[]:
#interpolate to prescribed length
l2_spl=US(scan_list,scan_list2)
X2ARRAY=np.repeat(l2_spl(XARRAY),precis)
else:
exit('Error: need to specify two scan_lists for 2d scan!')
#log- or linearly spaced 2d scan
if (scan_list==[] and scan_list2==[]):
if log:
XARRAY=np.tile(np.logspace(np.log10(scan_range[0]),np.log10(scan_range[1]),precis),precis2).flatten()
else:
XARRAY=np.tile(np.linspace(scan_range[0],scan_range[1],precis),precis2).flatten()
if log2:
X2ARRAY=np.repeat(np.logspace(np.log10(scan_range2[0]),np.log10(scan_range2[1]),precis2),precis).flatten()
else:
X2ARRAY=np.repeat(np.linspace(scan_range2[0],scan_range2[1],precis2),precis).flatten()
return XARRAY, X2ARRAY
def ffht(fEM, time, freq, ftarg):
"""Fourier Transform using a Cosine- or a Sine-filter.
It follows the Filter methodology [Anderson_1975]_, see `fht` for more
information.
The function is called from one of the modelling routines in :mod:`model`.
Consult these modelling routines for a description of the input and output
parameters.
This function is based on `get_CSEM1D_TD_FHT.m` from the source code
distributed with [Key_2012]_.
Returns
-------
tEM : array
Returns time-domain EM response of `fEM` for given `time`.
conv : bool
Only relevant for QWE/QUAD.
"""
# Get ftarg values
fftfilt, pts_per_dec, ftkind = ftarg
# Settings depending if cos/sin plus scaling
if ftkind == 'sin':
fEM = -fEM.imag
else:
fEM = fEM.real
if pts_per_dec: # Use pts_per_dec frequencies per decade
# 1. Interpolate in frequency domain
sfEM = iuSpline(np.log10(2*np.pi*freq), fEM)
ifEM = sfEM(np.log10(fftfilt.base/time[:, None]))
# 2. Filter
tEM = np.dot(ifEM, getattr(fftfilt, ftkind))
else: # Standard FHT procedure
# Get new times in frequency domain
_, itime = get_spline_values(fftfilt, time)
# Re-arranged fEM with shape (ntime, nfreq). Each row starts one
# 'freq' higher.
fEM = np.concatenate((np.tile(fEM, itime.size).squeeze(),
np.zeros(itime.size)))
fEM = fEM.reshape(itime.size, -1)[:, :fftfilt.base.size]
# 1. Filter
stEM = np.dot(fEM, getattr(fftfilt, ftkind))
# 2. Interpolate in time domain
itEM = iuSpline(np.log10((itime)[::-1]), stEM[::-1])
tEM = itEM(np.log10(time))
# Return the electromagnetic time domain field
# (Second argument is only for QWE)
return tEM/time, True
def __init__(self, data, categorical_parameters, continuous_parameters, func=None, order=1):
self.key_columns = categorical_parameters
self.parameter_columns = continuous_parameters
self.func = func
if len(self.parameter_columns) not in [1, 2]:
raise ValueError("Only interpolation over 1 or 2 variables is supported")
if len(self.parameter_columns) == 1 and order == 0:
raise ValueError("Order 0 only supported for 2d interpolation")
# These are the columns which the interpolation function will approximate
value_columns = sorted(data.columns.difference(set(self.key_columns)|set(self.parameter_columns)))
if self.key_columns:
# Since there are key_columns we need to group the table by those
# columns to get the sub-tables to fit
sub_tables = data.groupby(self.key_columns)
else:
# There are no key columns so we will fit the whole table
sub_tables = {None: data}.items()
self.interpolations = {}
for key, base_table in sub_tables:
if base_table.empty:
continue
# For each permutation of the key columns build interpolations
self.interpolations[key] = {}
for value_column in value_columns:
# For each value in the table build an interpolation function
if len(self.parameter_columns) == 2:
# 2 variable interpolation
if order == 0:
x = base_table[list(self.parameter_columns)]
y = base_table[value_column]
func = interpolate.NearestNDInterpolator(x=x.values, y=y.values)
else:
index, column = self.parameter_columns
table = base_table.pivot(index=index, columns=column, values=value_column)
x = table.index.values
y = table.columns.values
z = table.values
func = interpolate.RectBivariateSpline(x=x, y=y, z=z, ky=order, kx=order).ev
else:
# 1 variable interpolation
base_table = base_table.sort_values(by=self.parameter_columns[0])
x = base_table[self.parameter_columns[0]]
y = base_table[value_column]
func = interpolate.InterpolatedUnivariateSpline(x, y, k=order)
self.interpolations[key][value_column] = func
