Python numpy 模块，sinc() 实例源码

def generate_mtf(pixel_aperture=1, azimuth=0, num_samples=128):
    ''' generates the 1D diffraction-limited MTF for a given pixel size and azimuth.

    Args:
        pixel_aperture (`float`): aperture of the pixel, in microns.  Pixel is
            assumed to be square.

        azimuth (`float`): azimuth to retrieve the MTF at, in degrees.

        num_samples (`int`): number of samples in the output array.

    Returns:
        `tuple` containing:
            `numpy.ndarray`: array of units, in cy/mm.

            `numpy.ndarray`: array of MTF values (rel. 1.0).

    Notes:
        Azimuth is not actually implemented yet.
    '''
    pitch_unit = pixel_aperture / 1e3
    normalized_frequencies = np.linspace(0, 2, num_samples)
    otf = np.sinc(normalized_frequencies)
    mtf = np.abs(otf)
    return normalized_frequencies / pitch_unit, mtf

def smPhiDP(phiDP, ran):
    # smooth phiDP field and take derivative
    # calculate lanczos filter weights
    numRan = ran.shape[0]
    numK = 31
    fc = 0.015
    kt = np.linspace(-(numK-1)/2, (numK-1)/2, numK)
    w = np.sinc(2.*kt*fc)*(2.*fc)*np.sinc(kt/(numK/2))

    #smoothPhiDP = convolve1d(phiDP, w, axis=1, mode='constant', cval=-999.)
    smoothPhiDP = conv(phiDP, w)
    #smoothPhiDP = np.ma.masked_where(smoothPhiDP==-999., smoothPhiDP)

    return smoothPhiDP

# function for estimating kdp
#----------------------------------

def n_even_fcn(f, o, w, l):
    """Even case."""
    # Variables :
    k = np.array(range(0, int(l) + 1, 1)) + 0.5
    b = np.zeros(k.shape)

    # # Run Loop :
    for s in range(0, len(f), 2):
        m = (o[s + 1] - o[s]) / (f[s + 1] - f[s])
        b1 = o[s] - m * f[s]
        b = b + (m / (4 * np.pi * np.pi) * (np.cos(2 * np.pi * k * f[
            s + 1]) - np.cos(2 * np.pi * k * f[s])) / (
            k * k)) * abs(np.square(w[round((s + 1) / 2)]))
        b = b + (f[s + 1] * (m * f[s + 1] + b1) * np.sinc(2 * k * f[
            s + 1]) - f[s] * (m * f[s] + b1) * np.sinc(2 * k * f[s])) * abs(
            np.square(w[round((s + 1) / 2)]))

    a = (np.square(w[0])) * 4 * b
    h = 0.5 * np.concatenate((np.flipud(a), a))

    return h

def NevenFcn(F, M, W, L):  # N is even
    # Variables :
    k = np.array(range(0, int(L) + 1, 1)) + 0.5
    b = np.zeros(k.shape)

    # # Run Loop :
    for s in range(0, len(F), 2):
        m = (M[s + 1] - M[s]) / (F[s + 1] - F[s])
        b1 = M[s] - m * F[s]
        b = b + (m / (4 * np.pi * np.pi) * (np.cos(2 * np.pi * k * F[
            s + 1]) - np.cos(2 * np.pi * k * F[s])) / (
            k * k)) * abs(np.square(W[round((s + 1) / 2)]))
        b = b + (F[s + 1] * (m * F[s + 1] + b1) * np.sinc(2 * k * F[
          s + 1]) - F[s] * (m * F[s] + b1) * np.sinc(2 * k * F[s])) * abs(
            np.square(W[round((s + 1) / 2)]))

    a = (np.square(W[0])) * 4 * b
    h = 0.5 * np.concatenate((np.flipud(a), a))

    return h


####################################################################
# - Filt the signal :
####################################################################

def lpfls(N,wp,ws,W):
    M = (N-1)/2
    nq = np.arange(0,2*M+1)
    nb = np.arange(0,M+1)
    q = (wp/np.pi)*np.sinc((wp/np.pi)*nq) - W*(ws/np.pi)*np.sinc((ws/np.pi)*nq)
    b = (wp/np.pi)*np.sinc((wp/np.pi)*nb)
    b[0] = wp/np.pi
    q[0] = wp/np.pi + W*(1-ws/np.pi) # since sin(pi*n)/pi*n = 1, not 0
    b = b.transpose()

    Q1 = ln.toeplitz(q[0:M+1])
    Q2 = ln.hankel(q[0:M+1],q[M:])
    Q = Q1+Q2

    a = ln.solve(Q,b)
    h = list(nq)
    for i in nb:
        h[i] = 0.5*a[M-i]
        h[N-1-i] = h[i]
    h[M] = 2*h[M]
    hmax = max(np.absolute(h))
    for i in nq:
        h[i] = (8191/hmax)*h[i]
    return h

def bpfls(N,ws1,wp1,wp2,ws2,W):
    M = (N-1)/2
    nq = np.arange(0,2*M+1)
    nb = np.arange(0,M+1)
    q = W*np.sinc(nq) - (W*ws2/np.pi) * np.sinc(nq* (ws2/np.pi)) + (wp2/np.pi) * np.sinc(nq*(wp2/np.pi)) - (wp1/np.pi) * np.sinc(nq*(wp1/np.pi)) + (W*ws1/np.pi) * np.sinc(nq*(ws1/np.pi))
    b = (wp2/np.pi)*np.sinc((wp2/np.pi)*nb) - (wp1/np.pi)*np.sinc((wp1/np.pi)*nb)
    b[0] = wp2/np.pi - wp1/np.pi
    q[0] = W - W*ws2/np.pi + wp2/np.pi - wp1/np.pi + W*ws1/np.pi # since sin(pi*n)/pi*n = 1, not 0
    b = b.transpose()

    Q1 = ln.toeplitz(q[0:M+1])
    Q2 = ln.hankel(q[0:M+1],q[M:])
    Q = Q1+Q2

    a = ln.solve(Q,b)
    h = list(nq)
    for i in nb:
        h[i] = 0.5*a[M-i]
        h[N-1-i] = h[i]
    h[M] = 2*h[M]
    hmax = max(np.absolute(h))
    for i in nq:
        h[i] = (8191/hmax)*h[i]
    return h

def hpfls(N,ws,wp,W):
    M = (N-1)/2
    nq = np.arange(0,2*M+1)
    nb = np.arange(0,M+1)
    b = 1 - (wp/np.pi)* np.sinc(nb * wp/np.pi)
    b[0] = 1- wp/np.pi
    q = 1 - (wp/np.pi)* np.sinc(nq * wp/np.pi) + W * (ws/np.pi) * np.sinc(nq * ws/np.pi) # since sin(pi*n)/pi*n = 1, not 0
    q[0] = b[0] + W* ws/np.pi
    b = b.transpose()

    Q1 = ln.toeplitz(q[0:M+1])
    Q2 = ln.hankel(q[0:M+1],q[M:])
    Q = Q1+Q2

    a = ln.solve(Q,b)
    h = list(nq)
    for i in nb:
        h[i] = 0.5*a[M-i]
        h[N-1-i] = h[i]
    h[M] = 2*h[M]
    hmax = max(np.absolute(h))
    for i in nq:
        h[i] = (8191/hmax)*h[i]
    return h

def lanczosSubPixKernel( subPixShift, kernelShape=3, lobes=None  ):
    """
    Generate a kernel suitable for ni.convolve to subpixally shift an image.
    """
    kernelShape = np.array( [kernelShape], dtype='int' )
    if kernelShape.ndim == 1: # make it 2-D
        kernelShape = np.array( [kernelShape[0], kernelShape[0]], dtype='int' )

    if lobes is None:
        lobes = (kernelShape[0]+1)/2

    x_range = np.arange(-kernelShape[1]/2,kernelShape[1]/2)+1.0-subPixShift[1]
    x_range = ( 2.0 / kernelShape[1] ) * x_range 
    y_range = np.arange(-kernelShape[1]/2,kernelShape[0]/2)+1.0-subPixShift[0]
    y_range = ( 2.0 /kernelShape[0] ) * y_range
    [xmesh,ymesh] = np.meshgrid( x_range, y_range )


    lanczos_filt = np.sinc(xmesh * lobes) * np.sinc(xmesh) * np.sinc(ymesh * lobes) * np.sinc(ymesh)

    lanczos_filt = lanczos_filt / np.sum(lanczos_filt) # Normalize filter output
    return lanczos_filt

def sincinterp(x):
        """
        Sinc interpolation for computation of fractional transformations.
        As appears in :
        -https://github.com/audiolabs/frft/
        ----------
        Args:
            f       : (array) Complex valued input array
            a       : (float) Alpha factor
        Returns:
            ret     : (array) Real valued synthesised data
        """
        N = len(x)
        y = np.zeros(2 * N - 1, dtype=x.dtype)
        y[:2 * N:2] = x
        xint = fftconvolve( y[:2 * N], np.sinc(np.arange(-(2 * N - 3), (2 * N - 2)).T / 2),)
        return xint[2 * N - 3: -2 * N + 3]

def collect_pixel(self, pixel_i, frame_i,  k,j,i):
        #print pixel_i, k,j,i
        t0 = time.time()
        #px_data = np.random.rand()
        #px_data = t0 - self.prev_px
        x_hw = self.stage.settings.x_position.read_from_hardware(send_signal=False)
        y_hw = self.stage.settings.y_position.read_from_hardware(send_signal=False)

        theta = np.pi/10. * frame_i
        x = x_hw*np.cos(theta) - y_hw*np.sin(theta)
        y = x_hw*np.sin(theta) + y_hw*np.cos(theta)

        px_data = np.sinc((x-50)*0.05)**2 * np.sinc(0.05*(y-50))**2 #+ 0.05*np.random.random()
        #px_data = (x-xhw)**2 + ( y-yhw)**2
        #if px_data > 1:
        #    print('hw', x, xhw, y, yhw)
        self.display_image_map[k,j,i] = px_data
        if self.settings['save_h5']:
            self.test_data[frame_i, k,j,i] = px_data 
        time.sleep(self.settings['pixel_time'])
        #self.prev_px = t0

def design_windowed_sinc_lpf(fc, bw):
        N = Filter.get_filter_length_from_bandwidth(bw)

        # Compute sinc filter impulse response
        h = np.sinc(2 * fc * (np.arange(N) - (N - 1) / 2.))

        # We use blackman window function
        w = np.blackman(N)

        # Multiply sinc filter with window function
        h = h * w

        # Normalize to get unity gain
        h_unity = h / np.sum(h)

        return h_unity

def intensitiesFFIntraAtom(nq, dq, partNrs, nameList, ffDict):
    """ uses atomistic form factors """
    dIntensity = np.zeros((nq, 2), dtype=float)
    partInt = np.zeros((nq, len(partNrs) + 1))
    qList = np.zeros(nq)
    qList[:] = [float(i * dq) for i in range(nq)]
    dIntensity[:, 0] = qList[:]
    partInt[:, 0] = qList[:]
    # for j in range(0,len(dIntegrand)):
    #    r=dIntegrand[j,0]
    #    sinc=j0(q*r)
    k = 0
    formFacProd = np.zeros((nq, len(partNrs) + 1))
    # partNrsProd=getPartNrsProd(partNrs)
    # print "partNrsProd ", partNrsProd
    for i in range(nq):
        for k in range(1, len(partNrs) + 1):
            # print k
            formFacProd[i, k] = ff.fiveGaussian(
                ffDict[nameList[k - 1]], qList[i]) ** 2
            partInt[i, k] += partNrs[k - 1] * formFacProd[i, k]
    for i in range(nq):
        dIntensity[i, 1] = partInt[i, 1:].sum()
    return partInt, dIntensity

def iterate_l1(L, alpha, arg, beta, K, N, rr):
    oversample_ratio = (1.0 * K / N)

    for l1 in range(-L, L + 1):
        alf = alpha[abs(l1)] * 1.0
        if l1 < 0:
            alf = numpy.conj(alf)
    #             r1 = numpy.sinc(1.0*(arg+1.0*l1*beta)/(1.0*K/N))
        input_array = (arg + 1.0 * l1 * beta) / oversample_ratio
        r1 = dirichlet(input_array.astype(numpy.float32))
        rr = iterate_sum(rr, alf, r1)
    return rr

def nufft_r(om, N, J, K, alpha, beta):
    '''
    equation (30) of Fessler's paper

    '''
    def iterate_sum(rr, alf, r1):
        rr = rr + alf * r1
        return rr
    def iterate_l1(L, alpha, arg, beta, K, N, rr):
        oversample_ratio = (1.0 * K / N)
        import time
        t0=time.time()
        for l1 in range(-L, L + 1):
            alf = alpha[abs(l1)] * 1.0
    #         if l1 < 0:
    #             alf = numpy.conj(alf)
        #             r1 = numpy.sinc(1.0*(arg+1.0*l1*beta)/(1.0*K/N))
            input_array = (arg + 1.0 * l1 * beta) / oversample_ratio
            r1 = dirichlet(input_array)
            rr = iterate_sum(rr, alf, r1)
        return rr

    M = numpy.size(om)  # 1D size
    gam = 2.0 * numpy.pi / (K * 1.0)
    nufft_offset0 = nufft_offset(om, J, K)  # om/gam -  nufft_offset , [M,1]
    dk = 1.0 * om / gam - nufft_offset0  # om/gam -  nufft_offset , [M,1]
    arg = outer_sum(-numpy.arange(1, J + 1) * 1.0, dk)
    L = numpy.size(alpha) - 1
#     print('alpha',alpha)
    rr = numpy.zeros((J, M), dtype=numpy.float32)
    rr = iterate_l1(L, alpha, arg, beta, K, N, rr)
    return (rr, arg)

def analytic_ft(self, unit_x, unit_y):
        ''' Analytic fourier transform of a pixel aperture.

        Args:
            unit_x (`numpy.ndarray`): sample points in x axis.

            unit_y (`numpy.ndarray`): sample points in y axis.

        Returns:
            `numpy.ndarray`: 2D numpy array containing the analytic fourier transform.

        '''
        xq, yq = np.meshgrid(unit_x, unit_y)
        return (sinc(xq * self.size_x / 1e3) *
                sinc(yq * self.size_y / 1e3)).astype(config.precision)

def make_bin_weights(self):
        erb_max = hz2erb(self.sr/2.0)
        ngrid = 1000
        erb_grid = np.arange(ngrid) * erb_max / (ngrid - 1)
        hz_grid = (np.exp(erb_grid / 9.26) - 1) / 0.00437
        resp = np.zeros((ngrid, self.n_bins))
        for b in range(self.n_bins):
            w = self.widths[b]
            r =  (2.0 * w + 1.0) / self.sr * (hz_grid - self.hz_freqs[b])
            resp[:,b] = np.square(np.sinc(r)+ 0.5 * np.sinc(r + 1.0) + 0.5 * np.sinc(r  - 1.0));
        self.weights = np.dot(linalg.pinv(resp), np.ones((ngrid,1)))

def autocorrelation_function(self, r):
        """compute the real space autocorrelation function for the Gaussian random field model

        """
        # compute the cut-level parameter beta
        beta = np.sqrt(2) * erfinv(2 * (1-self.frac_volume) - 1)

        # the covariance of the GRF
        acf_psi = (np.exp(-r/self.corr_length) * (1 + r / self.corr_length)
                   * np.sinc(2*r/self.repeat_distance))

        # integral discretization. henning says: the resolution 1e-2 is ad hoc, test required,
        # the integrand has a (integrable) singularity for t=1 and acf_psi = 1, so an adaptive
        # discretization seems preferable -> TODO
        dt = 1e-2
        t = np.arange(0, 1, dt)

        # the gridded integrand, via change of integration variable
        # compared to the wp-2 docu, to enable array-based computation
        t_gridded, acf_psi_gridded = np.meshgrid(t, acf_psi)
        integrand_gridded = (acf_psi_gridded / np.sqrt(1 - (t_gridded * acf_psi_gridded)**2)
                             * np.exp( - beta**2 / (1 + t_gridded * acf_psi_gridded)))

        acf = 1.0 / (2 * np.pi) * np.trapz(integrand_gridded, x=t_gridded)

        return acf

    # ft not known analytically: deligate
    # def ft_autocorrelation_function(self, k):

def autocorrelation_function(self, r):
        """compute the real space autocorrelation function for the Teubner Strey model

        """

        acf = self.corr_func_at_origin * np.exp(-r/self.corr_length) * np.sinc(2*r/self.repeat_distance)
        return acf

def n_odd_fcn(f, o, w, l):
    """Odd case."""
    # Variables :
    b0 = 0
    m = np.array(range(int(l + 1)))
    k = m[1:len(m)]
    b = np.zeros(k.shape)

    # Run Loop :
    for s in range(0, len(f), 2):
        m = (o[s + 1] - o[s]) / (f[s + 1] - f[s])
        b1 = o[s] - m * f[s]
        b0 = b0 + (b1 * (f[s + 1] - f[s]) + m / 2 * (
            f[s + 1] * f[s + 1] - f[s] * f[s])) * abs(
            np.square(w[round((s + 1) / 2)]))
        b = b + (m / (4 * np.pi * np.pi) * (
            np.cos(2 * np.pi * k * f[s + 1]) - np.cos(2 * np.pi * k * f[s])
        ) / (k * k)) * abs(np.square(w[round((s + 1) / 2)]))
        b = b + (f[s + 1] * (m * f[s + 1] + b1) * np.sinc(2 * k * f[
            s + 1]) - f[s] * (m * f[s] + b1) * np.sinc(2 * k * f[s])) * abs(
            np.square(w[round((s + 1) / 2)]))

    b = np.insert(b, 0, b0)
    a = (np.square(w[0])) * 4 * b
    a[0] = a[0] / 2
    aud = np.flipud(a[1:len(a)]) / 2
    a2 = np.insert(aud, len(aud), a[0])
    h = np.concatenate((a2, a[1:] / 2))

    return h

def NoddFcn(F, M, W, L):  # N is odd
    # Variables :
    b0 = 0
    m = np.array(range(int(L + 1)))
    k = m[1:len(m)]
    b = np.zeros(k.shape)

    # Run Loop :
    for s in range(0, len(F), 2):
        m = (M[s + 1] - M[s]) / (F[s + 1] - F[s])
        b1 = M[s] - m * F[s]
        b0 = b0 + (b1 * (F[s + 1] - F[s]) + m / 2 * (
            F[s + 1] * F[s + 1] - F[s] * F[s])) * abs(
            np.square(W[round((s + 1) / 2)]))
        b = b + (m / (4 * np.pi * np.pi) * (
            np.cos(2 * np.pi * k * F[s + 1]) - np.cos(2 * np.pi * k * F[s])
        ) / (k * k)) * abs(np.square(W[round((s + 1) / 2)]))
        b = b + (F[s + 1] * (m * F[s + 1] + b1) * np.sinc(2 * k * F[
          s + 1]) - F[s] * (m * F[s] + b1) * np.sinc(2 * k * F[s])) * abs(
            np.square(W[round((s + 1) / 2)]))

    b = np.insert(b, 0, b0)
    a = (np.square(W[0])) * 4 * b
    a[0] = a[0] / 2
    aud = np.flipud(a[1:len(a)]) / 2
    a2 = np.insert(aud, len(aud), a[0])
    h = np.concatenate((a2, a[1:] / 2))

    return h


# Even case

def lpfls2notch(N,wp,ws,wn1,wn2,W):
    M = (N-1)/2
    nq = np.arange(0,2*M+1)
    nb = np.arange(0,M+1)
    q = (wp/np.pi)*np.sinc((wp/np.pi)*nq) - W*(ws/np.pi)*np.sinc((ws/np.pi)*nq)
    b = (wp/np.pi)*np.sinc((wp/np.pi)*nb)
    q[0] = wp/np.pi + W*(1-ws/np.pi) # since sin(pi*n)/pi*n = 1, not 0
    b = np.asmatrix(b)
    b = b.transpose()

    Q1 = ln.toeplitz(q[0:M+1])
    Q2 = ln.hankel(q[0:M+1],q[M:])
    Q = Q1+Q2

    G1 = np.cos(wn1*nb)
    G2 = np.cos(wn2*nb)
    G = np.matrix([G1,G2])

    d = np.array([0,0])
    d = np.asmatrix(d)
    d = d.transpose()

    c = np.asmatrix(ln.solve(Q,b))

    mu = ln.solve(G*ln.inv(Q)*G.transpose(),G*c - d)

    a = c - ln.solve(Q,G.transpose()*mu)
    h = np.zeros(N)
    for i in nb:
        h[i] = 0.5*a[M-i]
        h[N-1-i] = h[i]
    h[M] = 2*h[M]
    hmax = max(np.absolute(h))
    for i in nq:
        h[i] = (8191/hmax)*h[i]
    return h

def lpfls1notch(N,wp,ws,wn1,W):
    M = (N-1)/2
    nq = np.arange(0,2*M+1)
    nb = np.arange(0,M+1)
    q = (wp/np.pi)*np.sinc((wp/np.pi)*nq) - W*(ws/np.pi)*np.sinc((ws/np.pi)*nq)
    b = (wp/np.pi)*np.sinc((wp/np.pi)*nb)
    q[0] = wp/np.pi + W*(1-ws/np.pi) # since sin(pi*n)/pi*n = 1, not 0
    b = np.asmatrix(b)
    b = b.transpose()

    Q1 = ln.toeplitz(q[0:M+1])
    Q2 = ln.hankel(q[0:M+1],q[M:])
    Q = Q1+Q2

    G1 = np.cos(wn1*nb)
    G = np.matrix([G1])

    d = np.array([0])
    d = np.asmatrix(d)

    c = np.asmatrix(ln.solve(Q,b))

    mu = ln.solve(G*ln.inv(Q)*G.transpose(),G*c - d)

    a = c - ln.solve(Q,G.transpose()*mu)
    h = np.zeros(N)
    for i in nb:
        h[i] = 0.5*a[M-i]
        h[N-1-i] = h[i]
    h[M] = 2*h[M]
    hmax = max(np.absolute(h))
    for i in nq:
        h[i] = (8191/hmax)*h[i]
    return h

def rcosfir(beta, sps, span=None):
    """Generates a raised cosine FIR filter.

    :param beta: shape of the raised cosine filter (0-1)
    :param sps: number of samples per symbol
    :param span: length of the filter in symbols (None => automatic selection)

    >>> import arlpy
    >>> rc = arlpy.comms.rcosfir(0.25, 6)
    >>> bb = arlpy.comms.modulate(arlpy.comms.random_data(100), arlpy.comms.psk())
    >>> pb = arlpy.comms.upconvert(bb, 6, 27000, 18000, rc)
    """
    if beta < 0 or beta > 1:
        raise ValueError('Beta must be between 0 and 1')
    if span is None:
        # from http://www.commsys.isy.liu.se/TSKS04/lectures/3/MichaelZoltowski_SquareRootRaisedCosine.pdf
        # since this recommendation is for root raised cosine filter, it is conservative for a raised cosine filter
        span = 33-int(44*beta) if beta < 0.68 else 4
    delay = int(span*sps/2)
    t = _np.arange(-delay, delay+1, dtype=_np.float)/sps
    denom = 1 - (2*beta*t)**2
    eps = _np.finfo(float).eps
    idx1 = _np.nonzero(_np.abs(denom) > _sqrt(eps))
    b = _np.full_like(t, beta*_sin(_pi/(2*beta))/(2*sps))
    b[idx1] = _np.sinc(t[idx1]) * _cos(_pi*beta*t[idx1])/denom[idx1] / sps
    b /= _sqrt(_np.sum(b**2))
    return b

def lanczos(x, freq, radius=3):
    l = 0.5*freq * np.sinc(0.5*freq*x) * np.sinc(x/radius)
    l[(x < -radius) | (x > radius)] = 0.0
    return l

def pollData(self):
        """Poll for new data.  This method sleeps in order to ensure
        that self.pollSize observations are generated at a realistic rate.
        """
        # figure time between polls using exponentially weighted moving average
        curTime = time.time()
        if self.pollDelay >= 0.0:
            self.pollDelay = 0.8*self.pollDelay + 0.2*(curTime - self.lastPollTime)
        else:
            self.pollDelay = 0.0
            self.lastPollTime = curTime - self.shift
        sleepTime = np.max((0.0, self.shift - self.pollDelay))
        self.lastPollTime = curTime + sleepTime
        time.sleep(sleepTime)

        with self.lock:
            # generate some random data
            data = np.random.uniform(-self.scale.value, self.scale.value,
                                     size=(self.pollSize,self.nChan))

            erpEnd = self.erpStart.value +\
                self.erpSpeed.value *\
                data.shape[0]/float(self.sampRate)

            erp = np.linspace(self.erpStart.value, erpEnd, data.shape[0])
            erp = np.repeat(erp, data.shape[1]).reshape((-1,data.shape[1]))
            erp = erp * 0.5*(np.arange(data.shape[1])+1.0)
            erp = np.sinc(erp)

            data += erp

            self.erpStart.value = erpEnd

        return data

def lanczos(n, radius=3):
    taps = np.linspace(-radius,radius, n)
    return np.sinc(taps/radius)

def sinc(n, radius=3, freq=1.0):
    taps = np.linspace(-radius,radius, n)
    return freq * np.sinc(freq*taps)

def initImpulseResponse(self, window):
        if self.bandType == 'allpass':
            self.impulseResponse = windows.kroneckerDelta(self.order+1)

        elif self.bandType == 'allstop':
            self.impulseResponse = np.zeros_like(window)

        elif self.bandType == 'lowpass':
            hightaps = self.high*self.taps
            self.impulseResponse = self.high*np.sinc(hightaps) * window

        elif self.bandType == 'highpass':
            lowtaps = self.low*self.taps
            self.impulseResponse = (-self.low*np.sinc(lowtaps) * window +
                windows.kroneckerDelta(self.order+1))

        elif self.bandType == 'bandpass':
            lowtaps = self.low*self.taps
            hightaps = self.high*self.taps
            self.impulseResponse = (self.high*np.sinc(hightaps) -
                self.low*np.sinc(lowtaps)) * window

        elif self.bandType == 'bandstop':
            lowtaps = self.low*self.taps
            hightaps = self.high*self.taps
            self.impulseResponse = ((self.high*np.sinc(hightaps) -
                self.low*np.sinc(lowtaps)) * window +
                windows.kroneckerDelta(self.order+1))

        else:
            raise Exception('Invalid bandType: ' + str(self.bandType))

        self.impulseResponse = self.impulseResponse.astype(self.dtype, copy=False)

def BandpassFilter(lowFreq, highFreq, sampRate=1.0, order=None, filtType='butter', **kwargs):
    filtType = filtType.lower()
    if filtType in ('butter', 'cheby1', 'cheby2', 'ellip', 'bessel'):
        if order is None: order = 3
        return BandpassFilterIIR(lowFreq=lowFreq, highFreq=highFreq,
                    sampRate=sampRate, order=order, filtType=filtType, **kwargs)
    elif filtType in ('lanczos', 'sinc-blackman', 'sinc-hamming', 'sinc-hann'):
        if order is None: order = 20
        return BandpassFilterFIR(lowFreq=lowFreq, highFreq=highFreq,
                    sampRate=sampRate, order=order, filtType=filtType, **kwargs)
    else:
        raise Exception('Invalid filter type: ' + str(filtType) + '.')

def rootRaisedCosine(t):
    bit_period = 1/BIT_FREQUENCY

    if (t== bit_period/(4*BETA)):
        return (BETA/(np.pi*np.sqrt(2*bit_period)) * \
                ((np.pi + 2)*np.sin(np.pi/(4*BETA)) + (np.pi - 2)*np.cos(np.pi/(4*BETA))))
    else:
        return (4 * BETA / np.pi / np.sqrt(bit_period) * \
            (np.cos((1 + BETA) * np.pi * t / bit_period) + \
             (1 - BETA) * np.pi / (4 * BETA) * np.sinc((1-BETA)*t/bit_period)) / \
                (1 - (4*BETA*t/bit_period)**2))

def rootRaisedCosine(t):
    bit_period = 1/BIT_FREQUENCY

    if (t== bit_period/(4*BETA)):
        return (BETA/(np.pi*np.sqrt(2*bit_period)) * \
                ((np.pi + 2)*np.sin(np.pi/(4*BETA)) + (np.pi - 2)*np.cos(np.pi/(4*BETA))))
    else:
        return (4 * BETA / np.pi / np.sqrt(bit_period) * \
            (np.cos((1 + BETA) * np.pi * t / bit_period) + \
             (1 - BETA) * np.pi / (4 * BETA) * np.sinc((1-BETA)*t/bit_period)) / \
                (1 - (4*BETA*t/bit_period)**2))

def get_CML(self, q, t):
        """
        Calculate C, M, L forming the elements of T matrix
        :param q: a shifted coordinate grid
        :param t: time
        :return: tuple C, M, L
        """
        assert q is self.x_plus or q is self.x_minus, \
            "the shifted coordinate (q) must be either x_plus or x_minus"

        # get the difference of adiabatic potential curves
        Vg_minus_Ve = (self._Vg_plus_Ve_x_plus if q is self.x_plus else self._Vg_minus_Ve_x_minus)

        Veg = self.Veg(q, t)

        D = Veg**2 + 0.25*Vg_minus_Ve**2
        np.sqrt(D, out=D)

        S = np.sinc(D * self.dt / np.pi)
        S *= self.dt

        C = D * self.dt
        np.cos(C, out=C)

        M = S * Vg_minus_Ve
        M *= 0.5

        L = S * Veg

        return C, M, L

def ef_cascade(screenpos, i, nletters):
    v = np.array([0, -1])
    d = lambda t : 1 if t < 0 else abs(np.sinc(t) / (1 + t**4))
    return lambda t: screenpos + v * 400 * d(t - 0.15 * i)

def collect_pixel(self, pixel_i, k,j,i):
        #print pixel_i, k,j,i
        t0 = time.time()
        #px_data = np.random.rand()
        #px_data = t0 - self.prev_px
        x0,y0 = self.pos
        x_set = self.stage.settings['x_position']
        y_set = self.stage.settings['y_position']
        x_hw = self.stage.settings.x_position.read_from_hardware(send_signal=False)
        y_hw = self.stage.settings.y_position.read_from_hardware(send_signal=False)
        if np.abs(x_hw - x0) > 1:
            self.log.debug('='*60)
            self.log.debug('pos      {} {}'.format(x0, y0))
            self.log.debug('settings {} {}'.format(x_set, y_set))
            self.log.debug('hw       {} {}'.format(x_hw, y_hw))            
            self.log.debug('settings value delta {} {}'.format(x_set-x0, y_set-y0))
            self.log.debug('read_hw  value delta {} {}'.format(x_hw-x0, y_hw-y0))
            self.log.debug('='*60)

        x = x_hw
        y = y_hw

        px_data = np.sinc((x-50)*0.05)**2 * np.sinc(0.05*(y-50))**2 #+ 0.05*np.random.random()
        #px_data = (x-xhw)**2 + ( y-yhw)**2
        #if px_data > 1:
        #    print('hw', x, xhw, y, yhw)
        self.display_image_map[k,j,i] = px_data
        if self.settings['save_h5']:
            self.test_data[k,j,i] = px_data 
        time.sleep(self.settings['pixel_time'])
        #self.prev_px = t0

def Sinc(freq=440, amp=1.0, offset=0):
    """Makes a Sinc function.

    freq: float frequency in Hz
    amp: float amplitude, 1.0 is nominal max
    offset: float phase offset in radians

    returns: Sinusoid object
    """
    return Sinusoid(freq, amp, offset, func=np.sinc)

def intensitiesFFFaster(nq, dq, dIntegrand, keys, ffDict, dr):
    """ uses atomistic form factors """
    dIntensity = np.zeros((nq, 2), dtype=float)
    partInt = np.zeros((nq, len(dIntegrand[0])))
    nameList = [k.split(",") for k in keys]
    # print "nameList=",nameList
    qList = np.zeros(nq)
    qList[:] = [float(i * dq) for i in range(nq)]
    dIntensity[:, 0] = qList[:]
    partInt[:, 0] = qList[:]
    rList = dIntegrand[:, 0]
    # for j in range(0,len(dIntegrand)):
    #    r=dIntegrand[j,0]
    #    sinc=j0(q*r)
    formFacProd = np.zeros((nq, len(dIntegrand[0])))
    for i in range(nq):
        sincList = np.sinc(rList * qList[i] / math.pi) * dr
        for k in range(1, len(dIntegrand[0])):
            # print k
            formFacProd[i, k] = ff.fiveGaussian(ffDict[nameList[k - 1][0]], qList[i])\
                * ff.fiveGaussian(ffDict[nameList[k - 1][1]], qList[i])
            partInt[i, k] += (sincList[:] * dIntegrand[:, k]
                              ).sum() * formFacProd[i, k]
    for i in range(nq):
        dIntensity[i, 1] = partInt[i, 1:].sum()
    return partInt, dIntensity

def dirichlet(x):
    return numpy.sinc(x)

def dirichlet(x):
    return numpy.sinc(x)

def initLanczos(self, filtOrder):
        self.filtOrder = filtOrder

        if self.filtOrder % 2 != 0:
             raise Exception('Invalid filtOrder: ' + str(self.filtOrder) +
                 ' Must be an even integer.')

        radius = self.filtOrder // 2
        win = np.sinc(np.linspace(-radius, radius, self.filtOrder+1) / float(radius)) # lanczos
        #win = spsig.hamming(self.filtOrder+1) # sinc-hamming

        # this should be automated somehow XXX - idfah
        if self.filtOrder <= 6:
            cutoff = 2*0.570

        elif self.filtOrder <= 8:
            cutoff = 2*0.676

        elif self.filtOrder <= 12:
            cutoff = 2*0.781

        elif self.filtOrder <= 16:
            cutoff = 2*0.836

        elif self.filtOrder <= 32:
            cutoff = 2*0.918

        elif self.filtOrder <= 64:
            cutoff = 2*0.959

        # need to fix for multiple pool sizes XXX - idfah
        cutoff /= float(self.poolSize)

        taps = cutoff * np.linspace(-radius, radius, self.filtOrder+1, dtype=self.dtype)
        impulseResponse = cutoff * np.sinc(taps) * win

        self.filters = []
        nReadoutLayers = 1 if self.nHidden is None else 2
        for ni, no in self.layerDims[:-nReadoutLayers]:
            noEmb = no*(self.filtOrder+1) # no outs after filter embedding

            filtMat = np.zeros(noEmb*2, dtype=self.dtype)
            filtMat[noEmb-1::-no] = impulseResponse

            # filters strided for embedding
            sz = filtMat.itemsize
            filtMat = npst.as_strided(filtMat, (no,noEmb), strides=(sz,sz))[::-1].T
            self.filters.append(filtMat.copy())

def rrc(t):
    '''
    Input: T, evaluation point (seconds)
    Output: value of root-raised-cosine at time T
    '''
    # Delay between two bits
    bit_period = 1/BIT_FREQUENCY
    # Total amount of bits to transmit
    nb_bits = len(LIST_OF_BITS)
    # To be returned (sum of contributions)
    s = 0.0
    # Max value of rrc
    m = 4*BETA/np.pi/np.sqrt(bit_period) + (1-BETA)/np.sqrt(bit_period) + sum(abs(2*rootRaisedCosine(i*bit_period)) for i in range(1, TRUNCATION))

    if(t < - TRUNCATION * bit_period  or t >= (nb_bits + TRUNCATION) * bit_period):
        # T out of support
        r = 0.0
    else:
        # Bits that will affect function at time T
        relevant_bits = np.zeros(2*TRUNCATION+1)
        for i in range(2*TRUNCATION+1):
            j = t/bit_period + i - TRUNCATION
            j = int(j) if int(j) <= j else int(j) - 1
            if(j >= 0 and j < nb_bits):
                relevant_bits[i] = -1 if LIST_OF_BITS[j] == '0' else 1

        for i in range(2*TRUNCATION+1):
            tt = t/bit_period
            tt = t - int(tt)*bit_period if int(tt) <= tt else t - (int(tt)-1)*bit_period
            if(t == bit_period * (1 / 4 / BETA + (i - TRUNCATION))):
                # L'Hospital's rule because of potential discontinuity
                s += relevant_bits[i] * BETA / np.pi / np.sqrt(2*bit_period) * 1 / m * \
                     ((np.pi + 2) * np.sin(np.pi/4/BETA) + \
                      (np.pi - 2) * np.cos(np.pi/4/BETA))
            else:
                # General case formula
                s += relevant_bits[i] * 4*BETA/np.pi/np.sqrt(bit_period) * 1 / m * \
                     (np.cos((1 + BETA) * np.pi * ((tt / bit_period - (i-TRUNCATION)))) + \
                      (1 - BETA) * np.pi / 4 / BETA * \
                      np.sinc((1 - BETA) * (tt / bit_period - (i-TRUNCATION))))/ \
                      (1 - (4*BETA*(tt / bit_period - (i-TRUNCATION)))**2)

    return s

###   ###   ###   ###   ###   ###   ###