Python scipy.interpolate.splev() Examples

def sample_613():
    """
    6.1.3 插值
    :return:
    """
    from scipy.interpolate import interp1d, splrep, splev

    # 示例两种插值计算方式
    _, axs = plt.subplots(nrows=1, ncols=2, figsize=(14, 5))

    # 线性插值
    linear_interp = interp1d(x, y)
    # axs[0]左边的
    axs[0].set_title('interp1d')
    # 在相同坐标系下，同样的x，插值的y值使r.绘制（红色点）
    axs[0].plot(x, y, '', x, linear_interp(x), 'r.')

    # B-spline插值
    splrep_interp = splrep(x, y)
    # axs[1]右边的
    axs[1].set_title('splrep')
    # #在相同坐标系下，同样的x，插值的y值使g.绘制（绿色点）
    axs[1].plot(x, y, '', x, splev(x, splrep_interp), 'g.')
    plt.show() 

def curve_length(self, start=None, end=None, precision=0.01):
        '''
        Calculates the length of the curve by dividing the curve up
        into pieces of parameterized-length <precision>.
        '''
        if start is None: start = self.t[0]
        if end is None: end = self.t[-1]
        from scipy import interpolate
        if self.order == 1:
            # we just want to add up along the steps...
            ii = [ii for (ii,t) in enumerate(self.t) if start < t and t < end]
            ts = np.concatenate([[start], self.t[ii], [end]])
            xy = np.vstack([[self(start)], self.coordinates[:,ii].T, [self(end)]])
            return np.sum(np.sqrt(np.sum((xy[1:] - xy[:-1])**2, axis=1)))
        else:
            t = np.linspace(start, end, int(np.ceil((end-start)/precision)))
            dt = t[1] - t[0]
            dx = interpolate.splev(t, self.splrep[0], der=1)
            dy = interpolate.splev(t, self.splrep[1], der=1)
            return np.sum(np.sqrt(dx**2 + dy**2)) * dt 

def residual(pars, spforcorr, sptarget=None):
        # unpack parameters:
        #  extract .value attribute for each parameter
        xs = pars['xshift'].value
        yx = pars['yshift'].value
        yx2 = pars['yshiftcarr'].value
        spcorr = np.zeros((shape(spforcorr)))
    
        # we need to resample the spectra to compare them
        tck = interpolate.splrep(spforcorr[:,0]-xs,spforcorr[:,1]*yx+spforcorr[:,1]**2*yx2,s=0)
        spcorr[:,0] = spforcorr[:,0]    
        spcorr[:,1] = interpolate.splev(spforcorr[:,0],tck,der=0)
        
        if sptarget is None: #in such case we return the corrected spectrum
            return spcorr
        return (spcorr[:,1] - sptarget[:,1])
    ###########################################################################        
        
    # Now we choose the portion of spectra to fit 

def residual(pars, spforcorr, sptarget=None):
            # unpack parameters:
            #  extract .value attribute for each parameter
            xs = pars['xshift'].value
            yx = pars['yshift'].value
            yx2 = pars['yshiftcarr'].value
            spcorr = np.zeros((shape(spforcorr)))
        
            # we need to resample the spectra to compare them
            tck = interpolate.splrep(spforcorr[:,0]-xs,spforcorr[:,1]*yx+spforcorr[:,1]**2*yx2,s=0)
            spcorr[:,0] = spforcorr[:,0]    
            spcorr[:,1] = interpolate.splev(spcorr[:,0],tck,der=0)
            
            if sptarget is None: #in such case we return the corrected spectrum
                return spcorr
            return (spcorr[:,1] - sptarget[:,1])
        #######################################################################        
        
        # Now we choose the portion of spectra to fit 

def get_velocities(positions, times, tol):
    positions = np.atleast_2d(positions)
    n = len(positions)
    deg = min(3, n - 1)

    good_inds = np.r_[True, (abs(times[1:] - times[:-1]) >= 1e-6)]
    good_positions = positions[good_inds]
    good_times = times[good_inds]

    if len(good_inds) == 1:
        return np.zeros(positions[0:1].shape)

    (tck, _) = si.splprep(good_positions.T, s=tol ** 2 * (n + 1), u=good_times, k=deg)
    # smooth_positions = np.r_[si.splev(times,tck,der=0)].T
    velocities = np.r_[si.splev(times, tck, der=1)].T
    return velocities 

def test_Taitel_Dukler_splines():
    from scipy.interpolate import splrep, splev
    import numpy as np
    from fluids.two_phase import Dukler_XA_tck, Dukler_XC_tck, Dukler_XD_tck
    Dukler_XA_tck2 = splrep(np.log10(Dukler_XA_Xs), np.log10(Dukler_XA_As), s=5e-3, k=3)
    [assert_close(i, j) for i, j in zip(Dukler_XA_tck, Dukler_XA_tck2)]
#     XA_interp = UnivariateSpline(np.log10(Dukler_XA_Xs), np.log10(Dukler_XA_As), s=5e-3, k=3) # , ext='const'
#    XA_interp_obj = lambda x: 10**float(splev(log10(x), Dukler_XA_tck))
    
    Dukler_XD_tck2 = splrep(np.log10(Dukler_XD_Xs), np.log10(Dukler_XD_Ds), s=1e-2, k=3)
    [assert_close(i, j) for i, j in zip(Dukler_XD_tck, Dukler_XD_tck)]
#     XD_interp = UnivariateSpline(np.log10(Dukler_XD_Xs), np.log10(Dukler_XD_Ds), s=1e-2, k=3) # , ext='const'
#    XD_interp_obj = lambda x: 10**float(splev(log10(x), Dukler_XD_tck))
    
    Dukler_XC_tck2 = splrep(np.log10(Dukler_XC_Xs), np.log10(Dukler_XC_Cs), s=1e-3, k=3)
    [assert_close(i, j) for i, j in zip(Dukler_XC_tck, Dukler_XC_tck2)]
#     XC_interp = UnivariateSpline(np.log10(Dukler_XC_Xs), np.log10(Dukler_XC_Cs), s=1e-3, k=3) # ext='const'
#    XC_interp_obj = lambda x: 10**float(splev(log10(x), Dukler_XC_tck))
    
    # Curves look great to 1E-4! Also to 1E4. 

def pla(data, period=15):
    N = int(len(data)/period)
    orig_x = range(0,len(data))
    tck = splrep(orig_x, data,s=0)
    test_xs = np.linspace(0,len(data),N)
    spline_ys = splev(test_xs, tck)
    spline_yps = splev(test_xs, tck, der=1)
    xi = np.unique(tck[0])
    yi = [[splev(x, tck, der=j) for j in xrange(3)] for x in xi]
    P = interpolate.PiecewisePolynomial(xi,yi,orders=1)
    test_ys = P(test_xs)
    #inter_y = interp0(test_xs, test_ys, orig_x)
    inter_y = interp1(test_xs, test_ys, orig_x)
    
    mae = sqrt(mean_absolute_error(inter_y, data))
    #       mae = np.var(inter_y-data)
    return mae

#def paa(data, period=15): 

def pla(data, period=15):
    N = int(len(data)/period)
    orig_x = range(0,len(data))
    tck = splrep(orig_x, data,s=0)
    test_xs = np.linspace(0,len(data),N)
    spline_ys = splev(test_xs, tck)
    spline_yps = splev(test_xs, tck, der=1)
    xi = np.unique(tck[0])
    yi = [[splev(x, tck, der=j) for j in xrange(3)] for x in xi]
    P = interpolate.PiecewisePolynomial(xi,yi,orders=1)
    test_ys = P(test_xs)
    #inter_y = interp0(test_xs, test_ys, orig_x)
    inter_y = interp1(test_xs, test_ys, orig_x)
    
    mae = sqrt(mean_absolute_error(inter_y, data))
    #       mae = np.var(inter_y-data)
    return mae

#def paa(data, period=15): 

def interpolate(x, x_arr, y_arr, type='interp', order=3, left=None,
                right=None):
    """Perform interpolation on value x using arrays x_arr
       and y_arr.  The type of interpolate is defined by interp

       type:
       interp--use numpy.interp
       spline--use scipy.splrep and splev

       return
    """
    if type == 'interp':
        y = np.interp(x, x_arr, y_arr, left=left, right=right)
    if type == 'spline':
        if left is None:
            y_arr[0] = left
        if right is None:
            y_arr[-1] = right

        tk = scint.splrep(x_arr, y_arr, k=order)
        y = scint.splev(x, tk, der=0)

    return y 

def interpolate(self, s_freq=4, interp_method='cubic'):
        rri = self.rri
        irregular_time = self.time
        step = 1 / float(s_freq)
        regular_time = arange(0, irregular_time[-1] + step, step)
        
        if interp_method == 'cubic':
            tck = splrep(irregular_time, rri, s=0)
            rri_interp = splev(regular_time, tck, der=0)
            
        elif interp_method == 'linear':
            rri_interp = interp(regular_time, irregular_time, rri)
            
        self.time_interp = regular_time
        self.rri_interp = rri_interp
        self.s_freq = s_freq 

def interp_cubic_spline_qrs(qrs_index, qrs_amp, fs):
	time_qrs = qrs_index / float(ECG_RAW_FREQUENCY)
	time_qrs = time_qrs - time_qrs[0]
	time_qrs_interp = np.arange(0, 60, 1 / float(fs))
	tck = interpolate.splrep(time_qrs, qrs_amp, s=0)
	qrs_interp = interpolate.splev(time_qrs_interp, tck, der=0)
	return qrs_interp 

def test_radial_func1(self):
        r = np.linspace(0, 4, 300)
        f = np.exp(-r)
        o = SphericalOrbital(1, (r, f))
        str(o)
        def i_univariate(r, f):
            return interp.UnivariateSpline(r, f, k=5, s=0, ext=1, check_finite=False)
        def i_interp1d(r, f):
            return interp.interp1d(r, f, kind='cubic', fill_value=(f[0], 0.), bounds_error=False)
        def i_spline(r, f):
            from functools import partial
            tck = interp.splrep(r, f, k=5, s=0)
            return partial(interp.splev, tck=tck, der=0, ext=1)

        # Interpolation radius
        R = np.linspace(0, 5, 400)

        assert np.allclose(o.f(r), f)
        f_default = o.f(R)

        o.set_radial(r, f, interp=i_univariate)
        assert np.allclose(o.f(r), f)
        f_univariate = o.f(R)
        o.set_radial(r, f, interp=i_interp1d)
        assert np.allclose(o.f(r), f)
        f_interp1d = o.f(R)

        o.set_radial(r, f, interp=i_spline)
        assert np.allclose(o.f(r), f)
        f_spline = o.f(R)

        # Checks that they are equal
        assert np.allclose(f_univariate, f_interp1d)
        assert np.allclose(f_univariate, f_spline)
        assert np.allclose(f_univariate, f_default) 

def data_interp(self,x,y,xnew,k=3,der=0,s=0,**kwargs):
    '''Interpolates a dataset y(x) to y(xnew) using B-Splines of order k.'''
    from scipy import interpolate
    tck = interpolate.splrep(x,y,s=s,k=k)
    ynew = interpolate.splev(xnew,tck,der=der)
    
    return ynew 

def spline_fill(x, y, newindex, k=3, s=0):
        """
        s gives the smoothness, k is the degree of the spline
        k=1 "linear interpolation"; k=2 "quadratic spline"; k=3 "cubic spline"
        documentation says that use cubic spline by default

        This function will work for both interpolation and extrapolation
        However, using cubic or quadratic spline for extrapolation can give values far from
        the origional range.
        """
        # First line creates the relationship
        tck = interpolate.splrep(x, y, k=k, s=s)
        # Final line passes the new index to fill
        return interpolate.splev(newindex, tck, der=0) 

def compute_blackbody_response(self, Teffs=None):
        """
        Computes blackbody intensities across the entire range of
        effective temperatures. It does this for two regimes, energy-weighted
        and photon-weighted. It then fits a cubic spline to the log(I)-Teff
        values and exports the interpolation functions _log10_Inorm_bb_energy
        and _log10_Inorm_bb_photon.

        Arguments
        ----------
        * `Teffs` (array, optional, default=None): an array of effective
            temperatures. If None, a default array from ~300K to ~500000K with
            97 steps is used. The default array is uniform in log10 scale.
        """

        if Teffs is None:
            log10Teffs = np.linspace(2.5, 5.7, 97) # this corresponds to the 316K-501187K range.
            Teffs = 10**log10Teffs

        # Energy-weighted intensities:
        log10ints_energy = np.array([np.log10(self._bb_intensity(Teff, photon_weighted=False)) for Teff in Teffs])
        self._bb_func_energy = interpolate.splrep(Teffs, log10ints_energy, s=0)
        self._log10_Inorm_bb_energy = lambda Teff: interpolate.splev(Teff, self._bb_func_energy)

        # Photon-weighted intensities:
        log10ints_photon = np.array([np.log10(self._bb_intensity(Teff, photon_weighted=True)) for Teff in Teffs])
        self._bb_func_photon = interpolate.splrep(Teffs, log10ints_photon, s=0)
        self._log10_Inorm_bb_photon = lambda Teff: interpolate.splev(Teff, self._bb_func_photon)

        if 'blackbody:Inorm' not in self.content:
            self.content.append('blackbody:Inorm') 

def __call__(self, t, derivative=0):
        from scipy import interpolate
        xint = interpolate.splev(t, self.splrep[0], der=derivative, ext=0)
        yint = interpolate.splev(t, self.splrep[1], der=derivative, ext=0)
        return np.asarray([xint,yint]) 

def getCoonsPatchPointBez(bounds,x,y,width,height, densities = None):
    
    p00 = np.array(interpolate.splev( 0.0,bounds['s_top'])).flatten()
    p10 = np.array(interpolate.splev( 1.0,bounds['s_top'])).flatten()
    p11 = np.array(interpolate.splev( 0.0,bounds['s_bottom'])).flatten()
    p01 = np.array(interpolate.splev( 1.0,bounds['s_bottom'])).flatten()
    
    u = 1.0 * x / (width-1)
    v = 1.0 * y / (height-1)
    iu = 1.0 - u
    iv = 1.0 - v
    if densities is None:
        pu0 = np.array(interpolate.splev( u,bounds['s_top'])).flatten()
        pu1 = np.array(interpolate.splev(iu,bounds['s_bottom'])).flatten()
        pv0 = np.array(interpolate.splev(iv,bounds['s_left'])).flatten()
        pv1 = np.array(interpolate.splev( v,bounds['s_right'])).flatten()
    else:
        ut = 0.0
        ub = 0.0
        for i in range(x):
            ut+=densities['top'][i]
            ub+=densities['bottom'][i]
        vl = 0.0
        vr = 0.0
        for i in range(y):
            vl+=densities['left'][i]
            vr+=densities['right'][i]
        
        pu0 = np.array(interpolate.splev( ut,bounds['s_top'])).flatten()
        pu1 = np.array(interpolate.splev(1.0-ub,bounds['s_bottom'])).flatten()
        pv0 = np.array(interpolate.splev(1-0-vl,bounds['s_left'])).flatten()
        pv1 = np.array(interpolate.splev( vr,bounds['s_right'])).flatten()   
        
    return iv * pu0 + v * pu1 + iu * pv0 + u * pv1 - iu * iv * p00 - u * iv * p10 - iu * v * p01 - u * v * p11 

def get_delay(self,f_data,z_data,delay=None,ignoreslope=True,guess=True):
		'''
		ignoreslope option not used here
		retrieves the cable delay assuming the ideal resonance has a circular shape
		modifies the cable delay until the shape Im(S21) vs Re(S21) is circular
		see "do_calibration"
		'''
		maxval = np.max(np.absolute(z_data))
		z_data = z_data/maxval
		A1, A2, A3, A4, fr, Ql = self._fit_skewed_lorentzian(f_data,z_data)
		if self.df_error/fr > 0.0001 or self.dQl_error/Ql>0.1:
			#print("WARNING: Calibration using Lorentz fit failed, trying phase fit...")
			A1 = np.mean(np.absolute(z_data))
			A2 = 0.
			A3 = 0.
			A4 = 0.
			#fr = np.mean(f_data)
			f = splrep(f_data,np.unwrap(np.angle(z_data)),k=5,s=self.phasefitsmooth)
			fr = f_data[np.argmax(np.absolute(splev(f_data,f,der=1)))]
			Ql = 1e4
		if ignoreslope==True:
			A2 = 0.
		else:
			A2 = 0.
			print("WARNING: The ignoreslope option is ignored! Corrections to the baseline should be done manually prior to fitting.")
			print("see also: resonator_tools.calibration.fit_baseline_amp() etc. for help on fitting the baseline.")
			print("There is also an example ipython notebook for using this function.")
			print("However, make sure to understand the impact of the baseline (parasitic coupled resonances etc.) on your system.")
			#z_data = (np.absolute(z_data)-A2*(f_data-fr)) * np.exp(np.angle(z_data)*1j)  #usually not necessary
		if delay is None:
			if guess==True:
				delay = self._guess_delay(f_data,z_data)
			else:
				delay=0.
			delay = self._fit_delay(f_data,z_data,delay,maxiter=200)
		params = [A1, A2, A3, A4, fr, Ql]
		return delay, params 

def getBoundaryPoints(x , y):
    tck,u = interpolate.splprep([x, y], s=0, per=1)
    unew = np.linspace(u.min(), u.max(), 10000)
    xnew,ynew = interpolate.splev(unew, tck, der=0)
    tup = c_[xnew.astype(int),ynew.astype(int)].tolist()
    coord = list(set(tuple(map(tuple, tup))))
    coord = np.array([list(elem) for elem in coord])
    return coord[:,0],coord[:,1] 

def getBoundaryPoints(x = [], y = []):
	tck,u = interpolate.splprep([x, y], s=0, per=1)
	unew = np.linspace(u.min(), u.max(), 1000)
	xnew,ynew = interpolate.splev(unew, tck, der=0)
	tup = c_[xnew.astype(int),ynew.astype(int)].tolist()
	coord = list(set(tuple(map(tuple, tup))))
	coord = np.array([list(elem) for elem in coord])
	return coord[:,0],coord[:,1] 

def spline_up(x_lr, r):
  x_lr = x_lr.flatten()
  x_hr_len = len(x_lr) * r
  x_sp = np.zeros(x_hr_len)
  
  i_lr = np.arange(x_hr_len, step=r)
  i_hr = np.arange(x_hr_len)
  
  f = interpolate.splrep(i_lr, x_lr)

  x_sp = interpolate.splev(i_hr, f)

  return x_sp 

def spline_up(x_lr, r):
  x_lr = x_lr.flatten()
  x_hr_len = len(x_lr) * r
  x_sp = np.zeros(x_hr_len)
  
  i_lr = np.arange(x_hr_len, step=r)
  i_hr = np.arange(x_hr_len)
  
  f = interpolate.splrep(i_lr, x_lr)

  x_sp = interpolate.splev(i_hr, f)

  return x_sp 

def spline_up(x_lr, r):
  x_lr = x_lr.flatten()
  x_hr_len = len(x_lr) * r
  x_sp = np.zeros(x_hr_len)
  
  i_lr = np.arange(x_hr_len, step=r)
  i_hr = np.arange(x_hr_len)
  
  f = interpolate.splrep(i_lr, x_lr)

  x_sp = interpolate.splev(i_hr, f)

  return x_sp 

def spline_up(x_lr, r):
  x_lr = x_lr.flatten()
  x_hr_len = len(x_lr) * r
  x_sp = np.zeros(x_hr_len)
  
  i_lr = np.arange(x_hr_len, step=r)
  i_hr = np.arange(x_hr_len)
  
  f = interpolate.splrep(i_lr, x_lr)

  x_sp = interpolate.splev(i_hr, f)

  return x_sp 

def upsample(x_lr, r):
  x_lr = x_lr.flatten()
  x_hr_len = len(x_lr) * r
  x_sp = np.zeros(x_hr_len)
  
  i_lr = np.arange(x_hr_len, step=r)
  i_hr = np.arange(x_hr_len)
  
  f = interpolate.splrep(i_lr, x_lr)

  x_sp = interpolate.splev(i_hr, f)

  return x_sp 

def get_boundary_points(x, y):
    tck, u = interpolate.splprep([x, y], s=0, per=1)
    unew = np.linspace(u.min(), u.max(), 1000)
    xnew, ynew = interpolate.splev(unew, tck, der=0)
    tup = c_[xnew.astype(int), ynew.astype(int)].tolist()
    coord = list(set(tuple(map(tuple, tup))))
    coord = np.array([list(elem) for elem in coord])
    return np.array(coord[:, 0], dtype=np.int32), np.array(coord[:, 1], dtype=np.int32) 

def get_boundary_points(x, y):
    tck, u = interpolate.splprep([x, y], s=0, per=1)
    unew = np.linspace(u.min(), u.max(), 1000)
    xnew, ynew = interpolate.splev(unew, tck, der=0)
    tup = c_[xnew.astype(int), ynew.astype(int)].tolist()
    coord = list(set(tuple(map(tuple, tup))))
    coord = np.array([list(elem) for elem in coord])
    return np.array(coord[:, 0], dtype=np.int32), np.array(coord[:, 1], dtype=np.int32) 

def calcspec(x2, v, *a, **kwargs):
   """
   Doppler shifts the template and applies the polynom a using the Horner schema.
   p(x,a) * f(x, v)

   v : Radial velocity.
   a : Polynomial coefficients.
   fmod : If present, polynomial is applied on this model, otherwise it is loaded again.
   retpoly : If present the computed polynom is returned (instead of fmod * poly).
       This is needed e.g. to normalise the errors.
   """
   #fmod=cubicSpline.evalSpline(globvar.xt,globvar.yt,globvar.kt,x2*(1.0-v/c))
   #fmod=interpolate.splev(x2*(1.0-v/c),globvar.tck,der=0)
   #fmod=interpolate.splev(dopshift(x2,v),globvar.tck,der=0)
   #fmod = spline_ev(dopshift(x2,v), globvar.tck)
   # static variables calcspec.tck calcspec.wcen must be preset
   if not kwargs.get('retpoly'):
      fmod = kwargs.get('fmod', calcspec.tpl(dopshift(x2,v))) # parsed fmod or new
   x2 = x2 - calcspec.wcen
   # Use a Horner schema to evaluate the polynom
   # poly = a_0 + a_1 x + a_2 x^2 + ... = a_0 + x (a_1 + x (a_2+ ...))
   poly = a[-1]
   #for b in a[-2::-1]: poly = b+x2*poly  # add up the contributions
   for b in a[-2::-1]:  # add up the contributions p = ((a_3*x+a_2)*x+a_1)*x+a_0
      poly *= x2
      poly += b
   return poly if kwargs.get('retpoly') else fmod * poly 

def bspline(x, y, xref, smooth, deg_bspline=3, pz=1):
    """
    Bspline interpolation.

    The smoothing factor (s) is calculated based on an empirical formula (made by JCA, based on
    preliminary results) and is a function of pz, density of points and an input smoothing parameter (smooth). The
    formula is adjusted such that the parameter (smooth) produces similar smoothing results than a Hanning window with
    length smooth, as implemented in linear().

    :param x:
    :param y:
    :param xref:
    :param smooth: float: Smoothing factor. 0: no smoothing, 5: moderate smoothing, 50: large smoothing
    :param deg_bspline: int: Degree of spline
    :param pz: float: dimension of pixel along superior-inferior direction (z, assuming RPI orientation)
    :return:
    """
    # TODO: add flag to enforce boundaries, using weight flag in bspline function
    from scipy import interpolate
    if len(x) <= deg_bspline:
        deg_bspline -= 2
    density = (float(len(x)) / len(xref)) ** 2
    s = density * smooth * pz / float(3)
    logger.debug('Smoothing factor: smooth={}'.format(s))
    # Then, run bspline interpolation
    tck = interpolate.splrep(x, y, s=s, k=deg_bspline)
    y_fit = interpolate.splev(xref, tck, der=0)
    y_fit_der = interpolate.splev(xref, tck, der=1)
    return y_fit, y_fit_der 

def _evalfr(self, omega):
        """Evaluate a transfer function at a single angular frequency."""
        # Preallocate the output.
        if getattr(omega, '__iter__', False):
            out = empty((self.outputs, self.inputs, len(omega)), dtype=complex)
        else:
            out = empty((self.outputs, self.inputs), dtype=complex)

        if self.ifunc is None:
            try:
                out = self.fresp[:, :, where(self.omega == omega)[0][0]]
            except Exception:
                raise ValueError(
                    "Frequency %f not in frequency list, try an interpolating"
                    " FRD if you want additional points" % omega)
        else:
            if getattr(omega, '__iter__', False):
                for i in range(self.outputs):
                    for j in range(self.inputs):
                        for k, w in enumerate(omega):
                            frraw = splev(w, self.ifunc[i, j], der=0)
                            out[i, j, k] = frraw[0] + 1.0j * frraw[1]
            else:
                for i in range(self.outputs):
                    for j in range(self.inputs):
                        frraw = splev(omega, self.ifunc[i, j], der=0)
                        out[i, j] = frraw[0] + 1.0j * frraw[1]

        return out

    # Method for generating the frequency response of the system