python类sinc()的实例源码

detector.py 文件源码 项目:prysm 作者: brandondube 项目源码 文件源码 阅读 36 收藏 0 点赞 0 评论 0
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
nws_ppi.py 文件源码 项目:nexrad_basemap 作者: rskschrom 项目源码 文件源码 阅读 32 收藏 0 点赞 0 评论 0
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
#----------------------------------
filtering.py 文件源码 项目:tensorpac 作者: EtienneCmb 项目源码 文件源码 阅读 51 收藏 0 点赞 0 评论 0
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
_filtering.py 文件源码 项目:brainpipe 作者: EtienneCmb 项目源码 文件源码 阅读 28 收藏 0 点赞 0 评论 0
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 :
####################################################################
LSFIR.py 文件源码 项目:Least-Squared-Error-Based-FIR-Filters 作者: fourier-being 项目源码 文件源码 阅读 25 收藏 0 点赞 0 评论 0
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
LSFIR.py 文件源码 项目:Least-Squared-Error-Based-FIR-Filters 作者: fourier-being 项目源码 文件源码 阅读 46 收藏 0 点赞 0 评论 0
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
LSFIR.py 文件源码 项目:Least-Squared-Error-Based-FIR-Filters 作者: fourier-being 项目源码 文件源码 阅读 28 收藏 0 点赞 0 评论 0
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
zorro_util.py 文件源码 项目:zorro 作者: C-CINA 项目源码 文件源码 阅读 27 收藏 0 点赞 0 评论 0
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
TFMethods.py 文件源码 项目:ASP 作者: TUIlmenauAMS 项目源码 文件源码 阅读 22 收藏 0 点赞 0 评论 0
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]
base_raster_scan_test.py 文件源码 项目:ScopeFoundry 作者: ScopeFoundry 项目源码 文件源码 阅读 24 收藏 0 点赞 0 评论 0
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
Filter.py 文件源码 项目:urh 作者: jopohl 项目源码 文件源码 阅读 47 收藏 0 点赞 0 评论 0
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
pyntensities.py 文件源码 项目:capriqorn 作者: bio-phys 项目源码 文件源码 阅读 32 收藏 0 点赞 0 评论 0
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
pynufft_gpu.py 文件源码 项目:pynufft 作者: jyhmiinlin 项目源码 文件源码 阅读 31 收藏 0 点赞 0 评论 0
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
helper.py 文件源码 项目:pynufft 作者: jyhmiinlin 项目源码 文件源码 阅读 22 收藏 0 点赞 0 评论 0
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)
detector.py 文件源码 项目:prysm 作者: brandondube 项目源码 文件源码 阅读 24 收藏 0 点赞 0 评论 0
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)
qerbt.py 文件源码 项目:untwist 作者: IoSR-Surrey 项目源码 文件源码 阅读 32 收藏 0 点赞 0 评论 0
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)))
gaussian_random_field.py 文件源码 项目:smrt 作者: smrt-model 项目源码 文件源码 阅读 26 收藏 0 点赞 0 评论 0
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):
teubner_strey.py 文件源码 项目:smrt 作者: smrt-model 项目源码 文件源码 阅读 31 收藏 0 点赞 0 评论 0
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
filtering.py 文件源码 项目:tensorpac 作者: EtienneCmb 项目源码 文件源码 阅读 25 收藏 0 点赞 0 评论 0
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
_filtering.py 文件源码 项目:brainpipe 作者: EtienneCmb 项目源码 文件源码 阅读 26 收藏 0 点赞 0 评论 0
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
LSFIR.py 文件源码 项目:Least-Squared-Error-Based-FIR-Filters 作者: fourier-being 项目源码 文件源码 阅读 25 收藏 0 点赞 0 评论 0
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
LSFIR.py 文件源码 项目:Least-Squared-Error-Based-FIR-Filters 作者: fourier-being 项目源码 文件源码 阅读 22 收藏 0 点赞 0 评论 0
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
comms.py 文件源码 项目:arlpy 作者: org-arl 项目源码 文件源码 阅读 32 收藏 0 点赞 0 评论 0
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
func.py 文件源码 项目:cebl 作者: idfah 项目源码 文件源码 阅读 33 收藏 0 点赞 0 评论 0
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
erptest.py 文件源码 项目:cebl 作者: idfah 项目源码 文件源码 阅读 26 收藏 0 点赞 0 评论 0
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
windows.py 文件源码 项目:cebl 作者: idfah 项目源码 文件源码 阅读 31 收藏 0 点赞 0 评论 0
def lanczos(n, radius=3):
    taps = np.linspace(-radius,radius, n)
    return np.sinc(taps/radius)
windows.py 文件源码 项目:cebl 作者: idfah 项目源码 文件源码 阅读 25 收藏 0 点赞 0 评论 0
def sinc(n, radius=3, freq=1.0):
    taps = np.linspace(-radius,radius, n)
    return freq * np.sinc(freq*taps)
bandpass.py 文件源码 项目:cebl 作者: idfah 项目源码 文件源码 阅读 25 收藏 0 点赞 0 评论 0
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)
bandpass.py 文件源码 项目:cebl 作者: idfah 项目源码 文件源码 阅读 24 收藏 0 点赞 0 评论 0
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) + '.')
r_tuple_former.py 文件源码 项目:pdc-project 作者: ealain 项目源码 文件源码 阅读 36 收藏 0 点赞 0 评论 0
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))


问题


面经


文章

微信
公众号

扫码关注公众号