Python numpy.sinc() Examples

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 test_fourier_trafo_charfun_1d():
    # Characteristic function of [0, 1], its Fourier transform is
    # given by exp(-1j * y / 2) * sinc(y/2)
    def char_interval(x):
        return (x >= 0) & (x <= 1)

    def char_interval_ft(x):
        return np.exp(-1j * x / 2) * sinc(x / 2) / np.sqrt(2 * np.pi)

    # Base version
    discr = odl.uniform_discr(-2, 2, 40, impl='numpy')
    dft_base = FourierTransform(discr)

    # Complex version, should be as good
    discr = odl.uniform_discr(-2, 2, 40, impl='numpy', dtype='complex64')
    dft_complex = FourierTransform(discr)

    # Without shift
    discr = odl.uniform_discr(-2, 2, 40, impl='numpy', dtype='complex64')
    dft_complex_shift = FourierTransform(discr, shift=False)

    for dft in [dft_base, dft_complex, dft_complex_shift]:
        func_true_ft = dft.range.element(char_interval_ft)
        func_dft = dft(char_interval)
        assert (func_dft - func_true_ft).norm() < 5e-6 

def setUp(self):
        X_task_1 = np.random.rand(10, 2)
        y_task_1 = np.sinc(X_task_1 * 10 - 5).sum(axis=1)
        X_task_1 = np.concatenate((X_task_1, np.zeros([10, 1])), axis=1)

        X_task_2 = np.random.rand(10, 2)
        y_task_2 = np.sinc(X_task_2 * 2 - 4).sum(axis=1)
        X_task_2 = np.concatenate((X_task_2, np.ones([10, 1])), axis=1)

        X_task_3 = np.random.rand(10, 2)
        y_task_3 = np.sinc(X_task_3 * 8 - 6).sum(axis=1)
        X_task_3 = np.concatenate((X_task_3, 2 * np.ones([10, 1])), axis=1)

        self.X = np.concatenate((X_task_1, X_task_2, X_task_3), axis=0)
        self.y = np.concatenate((y_task_1, y_task_2, y_task_3), axis=0)
        self.model = WrapperBohamiann()
        self.model.train(self.X, self.y) 

def test_projected_incumbent_estimation(self):
        X = np.random.randn(20, 10)
        y = np.sinc(X).sum(axis=1)

        class DemoModel(object):
            def train(self, X, y):
                self.X = X
                self.y = y

            def predict(self, X):
                return self.y, np.ones(self.y.shape[0])

        model = DemoModel()
        model.train(X, y)
        inc, inc_val = incumbent_estimation.projected_incumbent_estimation(model, X, proj_value=1)
        b = np.argmin(y)

        assert inc[-1] == 1
        assert np.all(inc[:-1] == X[b])
        assert inc_val == y[b] 

def test_fourier_trafo_freq_shifted_charfun_1d():
    # Frequency-shifted characteristic function: mult. with
    # exp(-1j * b * x) corresponds to shifting the FT by b.
    def fshift_char_interval(x):
        return np.exp(-1j * x * np.pi) * ((x >= -0.5) & (x <= 0.5))

    def fshift_char_interval_ft(x):
        return sinc((x + np.pi) / 2) / np.sqrt(2 * np.pi)

    # Number of points is very important here (aliasing)
    discr = odl.uniform_discr(-2, 2, 400, impl='numpy',
                              dtype='complex64')
    dft = FourierTransform(discr)
    func_true_ft = dft.range.element(fshift_char_interval_ft)
    func_dft = dft(fshift_char_interval)
    assert (func_dft - func_true_ft).norm() < 0.001 

def test_dft_with_known_pairs_2d():

    # Frequency-shifted product of characteristic functions
    def fshift_char_rect(x):
        # Characteristic function of the cuboid
        # [-1, 1] x [1, 2]
        return (x[0] >= -1) & (x[0] <= 1) & (x[1] >= 1) & (x[1] <= 2)

    def fshift_char_rect_ft(x):
        # FT is a product of shifted and frequency-shifted sinc functions
        # 1st comp.: 2 * sinc(y)
        # 2nd comp.: exp(-1j * y * 3/2) * sinc(y/2)
        # Overall factor: (2 * pi)^(-1)
        return (
            2 * sinc(x[0])
            * np.exp(-1j * x[1] * 3 / 2) * sinc(x[1] / 2)
            / (2 * np.pi)
        )

    discr = odl.uniform_discr([-2] * 2, [2] * 2, (100, 400), impl='numpy',
                              dtype='complex64')
    dft = FourierTransform(discr)
    func_true_ft = dft.range.element(fshift_char_rect_ft)
    func_dft = dft(fshift_char_rect)
    assert (func_dft - func_true_ft).norm() < 0.001 

def lanczos(dx, a=3):
    """Lanczos kernel

    Parameters
    ----------
    dx: float
        amount to shift image
    a: int
        Lanczos window size parameter

    Returns
    -------
    result: array-like
        1D Lanczos kernel
    """
    if np.abs(dx) > 1:
        raise ValueError("The fractional shift dx must be between -1 and 1")
    window = np.arange(-a + 1, a + 1) + np.floor(dx)
    y = np.sinc(dx - window) * np.sinc((dx - window) / a)
    return y, window.astype(int) 

def test_fourier_trafo_scaling():
    # Test if the FT scales correctly

    # Characteristic function of [0, 1], its Fourier transform is
    # given by exp(-1j * y / 2) * sinc(y/2)
    def char_interval(x):
        return (x >= 0) & (x <= 1)

    def char_interval_ft(x):
        return np.exp(-1j * x / 2) * sinc(x / 2) / np.sqrt(2 * np.pi)

    discr = odl.uniform_discr(-2, 2, 40, impl='numpy', dtype='complex128')
    dft = FourierTransform(discr)

    for factor in (2, 1j, -2.5j, 1 - 4j):
        func_true_ft = factor * dft.range.element(char_interval_ft)
        func_dft = dft(factor * discr.element(char_interval))
        assert (func_dft - func_true_ft).norm() < 1e-6 

def diffuse_covar(num_bins, dist_mat, sr=16000, c=340, diag_eps=0.1):
    """
    Compute covarance matrices of the spherically isotropic noise field
        Gamma(omega)_{ij} = sinc(omega * tau_{ij}) = sinc(2 * pi * f * tau_{ij})
    Arguments:
        num_bins: number of the FFT points
        dist_mat: distance matrix
        sr: sample rate
        c: sound of the speed
    Return:
        covar: covariance matrix, F x N x N
    """
    N, _ = dist_mat.shape
    eps = np.eye(N) * diag_eps
    covar = np.zeros([num_bins, N, N])
    for f in range(num_bins):
        omega = np.pi * f * sr / (num_bins - 1)
        covar[f] = np.sinc(dist_mat * omega / c) + eps
    return covar 

def expmap_to_quaternion(e):
    """
    Convert axis-angle rotations (aka exponential maps) to quaternions.
    Stable formula from "Practical Parameterization of Rotations Using the Exponential Map".
    Expects a tensor of shape (*, 3), where * denotes any number of dimensions.
    Returns a tensor of shape (*, 4).
    """
    assert e.shape[-1] == 3
    
    original_shape = list(e.shape)
    original_shape[-1] = 4
    e = e.reshape(-1, 3)

    theta = np.linalg.norm(e, axis=1).reshape(-1, 1)
    w = np.cos(0.5*theta).reshape(-1, 1)
    xyz = 0.5*np.sinc(0.5*theta/np.pi)*e
    return np.concatenate((w, xyz), axis=1).reshape(original_shape) 

def test_fraunhofer_propagation_rectangular():
	for num_pix in [512, 1024]:
		pupil_grid = make_pupil_grid(num_pix)
		focal_grid = make_focal_grid(16, 8)

		for size in [[1,1], [0.75,1], [0.75,0.75]]:
			aperture = evaluate_supersampled(rectangular_aperture(size), pupil_grid, 8)

			for focal_length in [1, 1.3]:
				prop = FraunhoferPropagator(pupil_grid, focal_grid, focal_length=focal_length)

				for wavelength in [1, 0.8]:
					wf = Wavefront(aperture, wavelength)
					img = prop(wf).electric_field
					img /= img[np.argmax(np.abs(img))]

					k_x, k_y = np.array(size) / wavelength / focal_length
					reference = (np.sinc(k_x * focal_grid.x) * np.sinc(k_y * focal_grid.y))

					if num_pix == 512:
						assert np.abs(img - reference).max() < 5e-5
					elif num_pix == 1024:
						assert np.abs(img - reference).max() < 2e-5
					else:
						assert False # This should never happen. 

def conv_receiver_slit(self):
        """
        compute the rectangular convolution for the receiver slit or SiPSD
        pixel size

        Returns
        -------
        array-like
            the convolver
        """
        me = self.get_function_name()  # the name of the convolver, as a string
        # The receiver slit convolution is a top-hat of angular half-width
        # a=(slit_width/2)/diffractometer_radius
        # which has Fourier transform of sin(a omega)/(a omega)
        # NOTE! numpy's sinc(x) is sin(pi x)/(pi x), not sin(x)/x
        if self.param_dicts[me].get("slit_width", None) is None:
            return None

        epsr = (self.param_dicts["conv_receiver_slit"]["slit_width"] /
                self.param_dicts["conv_global"]["diffractometer_radius"])
        return self.general_tophat(me, epsr) 

def CompensateTSC(w, v):
    """
    Return the Fourier-space kernel that accounts for the convolution of
    the gridded field with the TSC window function in configuration space.

    .. note::
        see equation 18 (with p=3) of
        `Jing et al 2005 <https://arxiv.org/abs/astro-ph/0409240>`_

    Parameters
    ----------
    w : list of arrays
        the list of "circular" coordinate arrays, ranging from
        :math:`[-\pi, \pi)`.
    v : array_like
        the field array
    """
    for i in range(3):
        wi = w[i]
        tmp = (numpy.sinc(0.5 * wi / numpy.pi) ) ** 3
        v = v / tmp
    return v 

def CompensatePCS(w, v):
    """
    Return the Fourier-space kernel that accounts for the convolution of
    the gridded field with the PCS window function in configuration space.

    .. note::
        see equation 18 (with p=4) of
        `Jing et al 2005 <https://arxiv.org/abs/astro-ph/0409240>`_

    Parameters
    ----------
    w : list of arrays
        the list of "circular" coordinate arrays, ranging from
        :math:`[-\pi, \pi)`.
    v : array_like
        the field array
    """
    for i in range(3):
        wi = w[i]
        tmp = (numpy.sinc(0.5 * wi / numpy.pi) ) ** 4
        v = v / tmp
    return v 

def CompensateCIC(w, v):
    """
    Return the Fourier-space kernel that accounts for the convolution of
    the gridded field with the CIC window function in configuration space

    .. note::
        see equation 18 (with p=2) of
        `Jing et al 2005 <https://arxiv.org/abs/astro-ph/0409240>`_

    Parameters
    ----------
    w : list of arrays
        the list of "circular" coordinate arrays, ranging from
        :math:`[-\pi, \pi)`.
    v : array_like
        the field array
    """
    for i in range(3):
        wi = w[i]
        tmp = (numpy.sinc(0.5 * wi / numpy.pi) ) ** 2
        tmp[wi == 0.] = 1.
        v = v / tmp
    return v 

def val_exp(B_val):
    """
    Fast implementation of the translation and rotation specific exp function
    """
    t_val = imt_func(B_val, no.value)

    phiP_val = B_val - mult_with_ninf(t_val)
    phi = np.sqrt(-float(gmt_func(phiP_val, phiP_val)[0]))
    P_val = phiP_val / phi

    P_n_val = gmt_func(P_val, I3.value)
    t_nor_val = gmt_func(imt_func(t_val, P_n_val), P_n_val)
    t_par_val = t_val - t_nor_val

    coef_val = np.sin(phi) * P_val
    coef_val[0] += np.cos(phi)

    R_val = coef_val + gmt_func(coef_val, mult_with_ninf(t_nor_val)) + \
        np.sinc(phi/np.pi) * mult_with_ninf(t_par_val)
    return R_val 

def _get_taper_eigenvalues(tapers, half_bandwidth, time_index):
    '''Finds the eigenvalues of the original spectral concentration
    problem using the autocorr sequence technique from Percival and Walden,
    1993 pg 390

    Parameters
    ----------
    tapers : array, shape (n_tapers, n_time_samples_per_window)
    half_bandwidth : float
    time_index : array, (n_time_samples_per_window,)

    Returns
    -------
    eigenvalues : array, shape (n_tapers,)

    '''

    ideal_filter = 4 * half_bandwidth * np.sinc(
        2 * half_bandwidth * time_index)
    ideal_filter[0] = 2 * half_bandwidth
    n_time_samples_per_window = len(time_index)
    return np.dot(
        _auto_correlation(tapers)[:, :n_time_samples_per_window],
        ideal_filter) 

def fractional_delay(t0):
    """
    Creates a fractional delay filter using a windowed sinc function.
    The length of the filter is fixed by the module wide constant
    `frac_delay_length` (default 81).

    Parameters
    ----------
    t0: float
        The delay in fraction of sample. Typically between -1 and 1.

    Returns
    -------
    numpy array
        A fractional delay filter with specified delay.
    """

    N = constants.get("frac_delay_length")

    return np.hanning(N) * np.sinc(np.arange(N) - (N - 1) / 2 - t0) 

def main():
    """Test driver"""
    # From pp. 149-150
    x = np.ones(21)
    p = q = 1
    print('x: {}\np: {}\nq: {}'.format(x,p,q))
    b,a,err = prony(x, p, q)
    print('a: {}\nb: {}\nerr: {}'.format(a,b,err))
 
    # From pp. 152-153
    # Note that these results don't match the book, but they do match the
    # MATLAB version. So I'm either setting things up wrong or this is an
    # errata in the book.
    p = q = 5
    nd = 5
    n = np.arange(11)
    i = np.sinc((n-nd)/2)/2
    b,a,err = prony(i, p, q)
    print('a: {}\nb: {}\nerr: {}'.format(a,b,err)) 

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 sinc(x):
    # numpy.sinc scales by pi, we don't want that
    return np.sinc(x / np.pi)


# ---- DiscreteFourierTransform ---- # 

def test_sinx_div_x(self):

        def fun(x):
            return np.sin(x) / x

        for path in ['radial', 'spiral']:
            lim_f = Limit(fun, path=path, full_output=True)

            x = np.arange(-10, 10) / np.pi
            lim_f0, err = lim_f(x * np.pi)
            assert_array_almost_equal(lim_f0, np.sinc(x))
            assert np.all(err.error_estimate < 1.0e-14) 

def _interp_kernel_ft(norm_freqs, interp):
    """Scaled FT of a one-dimensional interpolation kernel.

    For normalized frequencies ``-1/2 <= xi <= 1/2``, this
    function returns::

        sinc(pi * xi)**k / sqrt(2 * pi)

    where ``k=1`` for 'nearest' and ``k=2`` for 'linear' interpolation.

    Parameters
    ----------
    norm_freqs : `numpy.ndarray`
        Normalized frequencies between -1/2 and 1/2
    interp : {'nearest', 'linear'}
        Type of interpolation kernel

    Returns
    -------
    ker_ft : `numpy.ndarray`
        Values of the kernel FT at the given frequencies
    """
    # Numpy's sinc(x) is equal to the 'math' sinc(pi * x)
    ker_ft = np.sinc(norm_freqs)
    interp_ = str(interp).lower()
    if interp_ == 'nearest':
        pass
    elif interp_ == 'linear':
        ker_ft *= ker_ft
    else:
        raise ValueError("`interp` '{}' not understood".format(interp))

    ker_ft /= np.sqrt(2 * np.pi)
    return ker_ft 

def test_sinc(self):
        q = [0., 3690., -270., 690.] * u.deg
        out = np.sinc(q)
        expected = np.sinc(q.to_value(u.radian)) * u.one
        assert isinstance(out, u.Quantity)
        assert np.all(out == expected)
        with pytest.raises(u.UnitsError):
            np.sinc(1.*u.one) 

def _resample_window_oct(p, q):
    """Port of Octave code to Python"""

    gcd = np.gcd(p, q)
    if gcd > 1:
        p /= gcd
        q /= gcd

    # Properties of the antialiasing filter
    log10_rejection = -3.0
    stopband_cutoff_f = 1. / (2 * max(p, q))
    roll_off_width = stopband_cutoff_f / 10

    # Determine filter length
    rejection_dB = -20 * log10_rejection
    L = np.ceil((rejection_dB - 8) / (28.714 * roll_off_width))

    # Ideal sinc filter
    t = np.arange(-L, L + 1)
    ideal_filter = 2 * p * stopband_cutoff_f \
        * np.sinc(2 * stopband_cutoff_f * t)

    # Determine parameter of Kaiser window
    if (rejection_dB >= 21) and (rejection_dB <= 50):
        beta = 0.5842 * (rejection_dB - 21)**0.4 \
            + 0.07886 * (rejection_dB - 21)
    elif rejection_dB > 50:
        beta = 0.1102 * (rejection_dB - 8.7)
    else:
        beta = 0.0

    # Apodize ideal filter response
    h = np.kaiser(2 * L + 1, beta) * ideal_filter

    return h 

def sinc2D(y, x):
    """
    2-D sinc function based on the product of 2 1-D sincs

    Parameters
    ----------
        x, y: arrays
            Coordinates where to evaluate the 2-D sinc
    Returns
    -------
    result: array
        2-D sinc evaluated in x and y
    """
    return np.dot(np.sinc(y), np.sinc(x)) 

def sinc(bw):
    width = 100 / bw
    return Waveform(bounds=(-0.5 * width, 0.5 * width, +np.inf),
                    seq=(_zero, _basic_wave(SINC, bw), _zero)) 

def low_pass_dirac(t0, alpha, Fs, N):
    """
    Creates a vector containing a lowpass Dirac of duration T sampled at Fs
    with delay t0 and attenuation alpha.

    If t0 and alpha are 2D column vectors of the same size, then the function
    returns a matrix with each line corresponding to pair of t0/alpha values.
    """

    return alpha * np.sinc(np.arange(N) - Fs * t0) 

def spec_correction(array, alpha=1.00, fix=True, cutoff=0.05, func=(None, None)):
        if array.ndim == 1:
            array = array[np.newaxis, ...]
        corr = np.empty_like(array)
        cent = [None] * len(array)
        for k, arr in enumerate(array):
            spec = arr / np.mean(arr)
            spec = filters.uniform_filter(spec, len(array) // 16, mode='wrap')
            offset = np.argmax(spec)
            spec = np.roll(spec, -offset)
            up_bound = 0
            while up_bound < len(spec) and spec[up_bound] > cutoff * spec[0]:
                up_bound += 1
            lo_bound = -1
            while lo_bound > -len(spec) and spec[lo_bound] > cutoff * spec[0]:
                lo_bound -= 1

            if func[0] == "Hamming" or func[0] is None:
                w = func[1] - (1.0 - func[1]) * np.cos(
                    2 * np.pi * np.arange(up_bound - lo_bound) / (up_bound - lo_bound - 1))
            elif func[0] == "Lanczos":
                w = np.sinc(2 * np.arange(up_bound - lo_bound) / (up_bound - lo_bound - 1) - 1)

            if fix is True:
                corr[k, ...] = 1.0 / spec
            else:
                corr[k, ...] = 1.0

            corr[k, up_bound - 1:lo_bound + 1] = 0.0

            corr[k, lo_bound:] *= w[0:-lo_bound]
            corr[k, :up_bound] *= w[-up_bound:]
            corr[k, ...] = np.roll(corr[k, ...], offset)
            corr[k, ...] /= np.mean(corr[k, ...])

            if offset > len(spec) // 2:
                offset -= len(spec)
            cent[k] = (up_bound + lo_bound) // 2 + offset
        return np.squeeze(corr), cent 

def extractRotorComponents(R):
    """
    Extracts the translation and rotation information from a TR rotor

    From :cite:`wareham-interpolation`.
    """
    phi = np.arccos(R[()])             # scalar
    phi2 = phi * phi                  # scalar
    # Notice: np.sinc(pi * x)/(pi x)
    phi_sinc = np.sinc(phi/np.pi)             # scalar
    phiP = ((R(2)*ninf)|ep)/(phi_sinc)
    t_normal_n = -((phiP * R(4))/(phi2 * phi_sinc))
    t_perpendicular_n = -(phiP * (phiP * R(2))(2))/(phi2 * phi_sinc)
    return phiP, t_normal_n, t_perpendicular_n