Python numpy.angle() Examples

The following are 30 code examples for showing how to use numpy.angle(). These examples are extracted from open source projects. You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example.

You may check out the related API usage on the sidebar.

You may also want to check out all available functions/classes of the module numpy , or try the search function .

Example 1
Project: clock-recovery   Author: mossmann   File: dsss-bpsk-reverse.py    License: MIT License 7 votes vote down vote up
def extract_chip_samples(samples):
    a = array(samples)
    f = scipy.fft(a*a)
    p = find_clock_frequency(abs(f))
    if 0 == p:
        return []
    cycles_per_sample = (p*1.0)/len(f)
    clock_phase = 0.25 + numpy.angle(f[p])/(tau)
    if clock_phase <= 0.5:
        clock_phase += 1
    chip_samples = []
    for i in range(len(a)):
        if clock_phase >= 1:
            clock_phase -= 1
            chip_samples.append(a[i])
        clock_phase += cycles_per_sample
    return chip_samples

# input: complex valued samples, FFT bin number of chip rate
#        input signal must be centered at 0 frequency
# output: number of chips found in repetitive chip sequence 
Example 2
Project: neuropythy   Author: noahbenson   File: retinotopy.py    License: GNU Affero General Public License v3.0 6 votes vote down vote up
def _retinotopic_field_sign_triangles(m, retinotopy):
    t = m.tess if isinstance(m, geo.Mesh) or isinstance(m, geo.Topology) else m
    # get the polar angle and eccen data as a complex number in degrees
    if pimms.is_str(retinotopy):
        (x,y) = as_retinotopy(retinotopy_data(m, retinotopy), 'geographical')
    elif retinotopy is Ellipsis:
        (x,y) = as_retinotopy(retinotopy_data(m, 'any'),      'geographical')
    else:
        (x,y) = as_retinotopy(retinotopy,                     'geographical')
    # Okay, now we want to make some coordinates...
    coords = np.asarray([x, y])
    us = coords[:, t.indexed_faces[1]] - coords[:, t.indexed_faces[0]]
    vs = coords[:, t.indexed_faces[2]] - coords[:, t.indexed_faces[0]]
    (us,vs) = [np.concatenate((xs, np.full((1, t.face_count), 0.0))) for xs in [us,vs]]
    xs = np.cross(us, vs, axis=0)[2]
    xs[np.isclose(xs, 0)] = 0
    return np.sign(xs) 
Example 3
Project: OpenFermion-Cirq   Author: quantumlib   File: fermionic_simulation.py    License: Apache License 2.0 6 votes vote down vote up
def _eigen_components(self):
        # projector onto subspace spanned by basis states with
        # Hamming weight != 2
        zero_component = np.diag(
            [int(bin(i).count('1') != 2) for i in range(16)])

        state_pairs = (('0110', '1001'), ('0101', '1010'), ('0011', '1100'))

        plus_minus_components = tuple(
            (-abs(weight) * sign / np.pi,
             state_swap_eigen_component(state_pair[0], state_pair[1], sign,
                                        np.angle(weight)))
            for weight, state_pair in zip(self.weights, state_pairs)
            for sign in (-1, 1))

        return ((0, zero_component),) + plus_minus_components 
Example 4
Project: galario   Author: mtazzari   File: test_galario.py    License: GNU Lesser General Public License v3.0 6 votes vote down vote up
def test_get_coords_meshgrid(nxy, inc, dxy, Dx, Dy, real_type, tol, acc_lib):

    ncol, nrow = nxy, nxy

    # create the referencemesh grid
    inc_cos = np.cos(inc)
    x = (np.linspace(0.5, -0.5 + 1./float(ncol), ncol, dtype=real_type)) * dxy * ncol
    y = (np.linspace(0.5, -0.5 + 1./float(nrow), nrow, dtype=real_type)) * dxy * nrow

    # we shrink the x axis, since PA is the angle East of North of the
    # the plane of the disk (orthogonal to the angular momentum axis)
    # PA=0 is a disk with vertical orbital node (aligned along North-South)
    x_m, y_m = np.meshgrid((x - Dx)/ inc_cos, y - Dy)
    R_m = np.sqrt(x_m ** 2. + y_m ** 2.)

    x_test, y_test, x_m_test, y_m_test, R_m_test = acc_lib.get_coords_meshgrid(nrow, ncol, dxy, inc, Dx=Dx, Dy=Dy, origin='upper')

    assert_allclose(x, x_test, atol=0, rtol=tol)
    assert_allclose(y, y_test, atol=0, rtol=tol)
    assert_allclose(x_m, x_m_test, atol=0, rtol=tol)
    assert_allclose(y_m, y_m_test, atol=0, rtol=tol)
    assert_allclose(R_m, R_m_test, atol=0, rtol=tol) 
Example 5
Project: ocelot   Author: ocelot-collab   File: wave.py    License: GNU General Public License v3.0 6 votes vote down vote up
def psi(self):
        # psi angle 0 - horizontal, pi/2 - vertical
        with np.errstate(divide='ignore'):
            psi = np.arctan(self.s2 / self.s1) / 2

        idx1 = np.where((self.s1 < 0) & (self.s2 > 0))
        idx2 = np.where((self.s1 < 0) & (self.s2 < 0))
        if np.size(psi) == 1:
            # continue
            # psi = psi
            if np.size(idx1): psi += np.pi / 2
            if np.size(idx2): psi -= np.pi / 2
        else:
            psi[idx1] += np.pi / 2
            psi[idx2] -= np.pi / 2
        return psi 
Example 6
Project: pyleecan   Author: Eomys   File: comp_angle_opening.py    License: Apache License 2.0 6 votes vote down vote up
def comp_angle_opening(self):
    """Compute the average opening angle of the Slot

    Parameters
    ----------
    self : Slot
        A Slot object

    Returns
    -------
    alpha: float
        Average opening angle of the slot [rad]

    """

    line_list = self.build_geometry()
    Z1 = line_list[0].get_begin()
    Z2 = line_list[-1].get_end()

    return angle(Z2) - angle(Z1) 
Example 7
Project: Griffin_lim   Author: candlewill   File: audio.py    License: MIT License 5 votes vote down vote up
def _griffin_lim(S):
    angles = np.exp(2j * np.pi * np.random.rand(*S.shape))
    S_complex = np.abs(S).astype(np.complex)
    for i in range(hparams.griffin_lim_iters):
        if i > 0:
            angles = np.exp(1j * np.angle(_stft(y)))
        y = _istft(S_complex * angles)
    return y 
Example 8
Project: neuropythy   Author: noahbenson   File: retinotopy.py    License: GNU Affero General Public License v3.0 5 votes vote down vote up
def calc_registration(preregistration_map, anchors,
                      max_steps=2000, max_step_size=0.05, method='random'):
    '''
    calc_registration is a calculator that creates the registration coordinates.
    '''
    # if max steps is a tuple (max, stride) then a trajectory is saved into
    # the registered_map meta-data
    pmap = preregistration_map
    if is_tuple(max_steps) or is_list(max_steps):
        (max_steps, stride) = max_steps
        traj = [preregistration_map.coordinates]
        x = preregistration_map.coordinates
        for s in np.arange(0, max_steps, stride):
            x = mesh_register(
                preregistration_map,
                [['edge',      'harmonic',      'scale', 1.0],
                 ['angle',     'infinite-well', 'scale', 1.0],
                 ['perimeter', 'harmonic'],
                 anchors],
                initial_coordinates=x,
                method=method,
                max_steps=stride,
                max_step_size=max_step_size)
            traj.append(x)
        pmap = pmap.with_meta(trajectory=np.asarray(traj))
    else:
        x = mesh_register(
            preregistration_map,
            [['edge',      'harmonic',      'scale', 1.0],
             ['angle',     'infinite-well', 'scale', 1.0],
             ['perimeter', 'harmonic'],
             anchors],
            method=method,
            max_steps=max_steps,
            max_step_size=max_step_size)
    return pmap.copy(coordinates=x) 
Example 9
Project: clock-recovery   Author: mossmann   File: wpcr.py    License: MIT License 5 votes vote down vote up
def wpcr(a):
    if len(a) < 4:
        return []
    b = (a > midpoint(a)) * 1.0
    d = numpy.diff(b)**2
    if len(numpy.argwhere(d > 0)) < 2:
        return []
    f = scipy.fft(d, len(a))
    p = find_clock_frequency(abs(f))
    if p == 0:
        return []
    cycles_per_sample = (p*1.0)/len(f)
    clock_phase = 0.5 + numpy.angle(f[p])/(tau)
    if clock_phase <= 0.5:
        clock_phase += 1
    symbols = []
    for i in range(len(a)):
        if clock_phase >= 1:
            clock_phase -= 1
            symbols.append(a[i])
        clock_phase += cycles_per_sample
    if debug:
        print("peak frequency index: %d / %d" % (p, len(f)))
        print("samples per symbol: %f" % (1.0/cycles_per_sample))
        print("clock cycles per sample: %f" % (cycles_per_sample))
        print("clock phase in cycles between 1st and 2nd samples: %f" % (clock_phase))
        print("clock phase in cycles at 1st sample: %f" % (clock_phase - cycles_per_sample/2))
        print("symbol count: %d" % (len(symbols)))
    return symbols

# convert soft symbols into bits (assuming binary symbols) 
Example 10
Project: OpenFermion-Cirq   Author: quantumlib   File: bogoliubov_transform.py    License: Apache License 2.0 5 votes vote down vote up
def _slater_basis_change(qubits: Sequence[cirq.Qid],
                         transformation_matrix: numpy.ndarray,
                         initially_occupied_orbitals: Optional[Sequence[int]]
                         ) -> cirq.OP_TREE:
    n_qubits = len(qubits)

    if initially_occupied_orbitals is None:
        decomposition, diagonal = givens_decomposition_square(
                transformation_matrix)
        circuit_description = list(reversed(decomposition))
        # The initial state is not a computational basis state so the
        # phases left on the diagonal in the decomposition matter
        yield (cirq.rz(rads=numpy.angle(diagonal[j])).on(qubits[j])
               for j in range(n_qubits))
    else:
        initially_occupied_orbitals = cast(
                Sequence[int], initially_occupied_orbitals)
        transformation_matrix = transformation_matrix[
                list(initially_occupied_orbitals)]
        n_occupied = len(initially_occupied_orbitals)
        # Flip bits so that the first n_occupied are 1 and the rest 0
        initially_occupied_orbitals_set = set(initially_occupied_orbitals)
        yield (cirq.X(qubits[j]) for j in range(n_qubits)
               if (j < n_occupied) != (j in initially_occupied_orbitals_set))
        circuit_description = slater_determinant_preparation_circuit(
                transformation_matrix)

    yield _ops_from_givens_rotations_circuit_description(
            qubits, circuit_description) 
Example 11
Project: OpenFermion-Cirq   Author: quantumlib   File: bogoliubov_transform.py    License: Apache License 2.0 5 votes vote down vote up
def _gaussian_basis_change(qubits: Sequence[cirq.Qid],
                           transformation_matrix: numpy.ndarray,
                           initially_occupied_orbitals: Optional[Sequence[int]]
                           ) -> cirq.OP_TREE:
    n_qubits = len(qubits)

    # Rearrange the transformation matrix because the OpenFermion routine
    # expects it to describe annihilation operators rather than creation
    # operators
    left_block = transformation_matrix[:, :n_qubits]
    right_block = transformation_matrix[:, n_qubits:]
    transformation_matrix = numpy.block(
            [numpy.conjugate(right_block), numpy.conjugate(left_block)])

    decomposition, left_decomposition, _, left_diagonal = (
        fermionic_gaussian_decomposition(transformation_matrix))

    if (initially_occupied_orbitals is not None and
            len(initially_occupied_orbitals) == 0):
        # Starting with the vacuum state yields additional symmetry
        circuit_description = list(reversed(decomposition))
    else:
        if initially_occupied_orbitals is None:
            # The initial state is not a computational basis state so the
            # phases left on the diagonal in the Givens decomposition matter
            yield (cirq.rz(rads=
                       numpy.angle(left_diagonal[j])).on(qubits[j])
                   for j in range(n_qubits))
        circuit_description = list(reversed(decomposition + left_decomposition))

    yield _ops_from_givens_rotations_circuit_description(
            qubits, circuit_description) 
Example 12
Project: OpenFermion-Cirq   Author: quantumlib   File: fermionic_simulation.py    License: Apache License 2.0 5 votes vote down vote up
def _arg(x):
    if x == 0:
        return 0
    if cirq.is_parameterized(x):
        return sympy.arg(x)
    return np.angle(x) 
Example 13
Project: NeuroKit   Author: neuropsychology   File: signal_synchrony.py    License: MIT License 5 votes vote down vote up
def _signal_synchrony_hilbert(signal1, signal2):

    hill1 = scipy.signal.hilbert(signal1)
    hill2 = scipy.signal.hilbert(signal2)

    phase1 = np.angle(hill1, deg=False)
    phase2 = np.angle(hill2, deg=False)
    synchrony = 1 - np.sin(np.abs(phase1 - phase2) / 2)

    return synchrony 
Example 14
Project: python-control   Author: python-control   File: discrete_test.py    License: BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def test_discrete_bode(self):
        # Create a simple discrete time system and check the calculation
        sys = TransferFunction([1], [1, 0.5], 1)
        omega = [1, 2, 3]
        mag_out, phase_out, omega_out = bode(sys, omega)
        H_z = list(map(lambda w: 1./(np.exp(1.j * w) + 0.5), omega))
        np.testing.assert_array_almost_equal(omega, omega_out)
        np.testing.assert_array_almost_equal(mag_out, np.absolute(H_z))
        np.testing.assert_array_almost_equal(phase_out, np.angle(H_z)) 
Example 15
Project: Tacotron-pytorch   Author: soobinseo   File: data.py    License: Apache License 2.0 5 votes vote down vote up
def _griffin_lim(S):
    '''librosa implementation of Griffin-Lim
    Based on https://github.com/librosa/librosa/issues/434
    '''
    angles = np.exp(2j * np.pi * np.random.rand(*S.shape))
    S_complex = np.abs(S).astype(np.complex)
    y = _istft(S_complex * angles)
    for i in range(hp.griffin_lim_iters):
        angles = np.exp(1j * np.angle(_stft(y)))
        y = _istft(S_complex * angles)
    return y 
Example 16
Project: recruit   Author: Frank-qlu   File: test_core.py    License: Apache License 2.0 5 votes vote down vote up
def test_basic_ufuncs(self):
        # Test various functions such as sin, cos.
        (x, y, a10, m1, m2, xm, ym, z, zm, xf) = self.d
        assert_equal(np.cos(x), cos(xm))
        assert_equal(np.cosh(x), cosh(xm))
        assert_equal(np.sin(x), sin(xm))
        assert_equal(np.sinh(x), sinh(xm))
        assert_equal(np.tan(x), tan(xm))
        assert_equal(np.tanh(x), tanh(xm))
        assert_equal(np.sqrt(abs(x)), sqrt(xm))
        assert_equal(np.log(abs(x)), log(xm))
        assert_equal(np.log10(abs(x)), log10(xm))
        assert_equal(np.exp(x), exp(xm))
        assert_equal(np.arcsin(z), arcsin(zm))
        assert_equal(np.arccos(z), arccos(zm))
        assert_equal(np.arctan(z), arctan(zm))
        assert_equal(np.arctan2(x, y), arctan2(xm, ym))
        assert_equal(np.absolute(x), absolute(xm))
        assert_equal(np.angle(x + 1j*y), angle(xm + 1j*ym))
        assert_equal(np.angle(x + 1j*y, deg=True), angle(xm + 1j*ym, deg=True))
        assert_equal(np.equal(x, y), equal(xm, ym))
        assert_equal(np.not_equal(x, y), not_equal(xm, ym))
        assert_equal(np.less(x, y), less(xm, ym))
        assert_equal(np.greater(x, y), greater(xm, ym))
        assert_equal(np.less_equal(x, y), less_equal(xm, ym))
        assert_equal(np.greater_equal(x, y), greater_equal(xm, ym))
        assert_equal(np.conjugate(x), conjugate(xm)) 
Example 17
Project: recruit   Author: Frank-qlu   File: test_umath_complex.py    License: Apache License 2.0 5 votes vote down vote up
def test_simple(self):
        x = np.array([1+0j, 1+2j])
        y_r = np.log(np.abs(x)) + 1j * np.angle(x)
        y = np.log(x)
        for i in range(len(x)):
            assert_almost_equal(y[i], y_r[i]) 
Example 18
Project: galario   Author: mtazzari   File: utils.py    License: GNU Lesser General Public License v3.0 5 votes vote down vote up
def py_sampleProfile(intensity, Rmin, dR, nxy, dxy, udat, vdat, dRA=0., dDec=0., PA=0, inc=0.):
    """
    Python implementation of sampleProfile.

    """
    inc_cos = np.cos(inc)

    nrad = len(intensity)
    gridrad = np.linspace(Rmin, Rmin + dR * (nrad - 1), nrad)

    ncol, nrow = nxy, nxy
    # create the mesh grid
    x = (np.linspace(0.5, -0.5 + 1./float(ncol), ncol)) * dxy * ncol
    y = (np.linspace(0.5, -0.5 + 1./float(nrow), nrow)) * dxy * nrow

    # we shrink the x axis, since PA is the angle East of North of the
    # the plane of the disk (orthogonal to the angular momentum axis)
    # PA=0 is a disk with vertical orbital node (aligned along North-South)
    x_axis, y_axis = np.meshgrid(x / inc_cos, y)
    x_meshgrid = np.sqrt(x_axis ** 2. + y_axis ** 2.)

    # convert to Jansky
    sr_to_px = dxy**2.
    intensity *= sr_to_px
    f = interp1d(gridrad, intensity, kind='linear', fill_value=0.,
                 bounds_error=False, assume_sorted=True)
    intensmap = f(x_meshgrid)

    intensmap[nrow//2, ncol//2] = central_pixel(intensity, Rmin, dR, dxy)

    vis = py_sampleImage(intensmap, dxy, udat, vdat, PA=PA, dRA=dRA, dDec=dDec)

    return vis 
Example 19
Project: galario   Author: mtazzari   File: test_galario.py    License: GNU Lesser General Public License v3.0 5 votes vote down vote up
def test_interpolate(size, real_type, complex_type, rtol, atol, acc_lib):
    """
    Test the interpolation of the output FT.

    """
    nsamples = 10000
    maxuv = 1000.

    reference_image = create_reference_image(size=size, dtype=real_type)
    udat, vdat = create_sampling_points(nsamples, maxuv/2.2)
    # this factor has to be > than 2 because the matrix cover between -maxuv/2 to +maxuv/2,
    # therefore the sampling points have to be contained inside.

    udat = udat.astype(real_type)
    vdat = vdat.astype(real_type)

    # no rotation
    du = maxuv/size
    uroti, vroti = uv_idx_r2c(udat, vdat, du, size/2.)

    uroti = uroti.astype(real_type)
    vroti = vroti.astype(real_type)

    ft = np.fft.fftshift(np.fft.fft2(np.fft.fftshift(reference_image))).astype(complex_type, order='C')

    ReInt = int_bilin_MT(ft.real, uroti, vroti)
    ImInt = int_bilin_MT(ft.imag, uroti, vroti)
    AmpInt = int_bilin_MT(np.abs(ft), uroti, vroti)
    uneg = udat < 0.
    ImInt[uneg] *= -1.
    PhaseInt = np.angle(ReInt + 1j*ImInt)
    
    ReInt = AmpInt * np.cos(PhaseInt)
    ImInt = AmpInt * np.sin(PhaseInt)

    complexInt = acc_lib.interpolate(ft, du,
                                     udat.astype(real_type),
                                     vdat.astype(real_type))

    assert_allclose(ReInt, complexInt.real, rtol, atol)
    assert_allclose(ImInt, complexInt.imag, rtol, atol) 
Example 20
Project: magpy   Author: geomagpy   File: emd.py    License: BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
def calc_inst_info(modes,samplerate):
    """
Calculate the instantaneous frequency, amplitude, and phase of
each mode.
    """

    amp=np.zeros(modes.shape,np.float32)
    phase=np.zeros(modes.shape,np.float32)
    f=np.zeros(modes.shape,np.float32)

    print("Mode 1:", len(modes), samplerate)

    for m in range(len(modes)):
        h=scipy.signal.hilbert(modes[m])
        print(len(modes[m]))
        print("Mean Amplitude of mode ", m, np.mean(np.abs(h)))
        print("Mean Phase of mode ", m, np.mean(np.angle(h)))
        phase[m,:]=np.angle(h)
        print("Frequ", np.diff(np.unwrap(phase[:,np.r_[0,0:len(modes[m])]]))/(2*np.pi)*samplerate)
        amp[m,:]=np.abs(h)
        phase[m,:]=np.angle(h)
        f[m,:] = np.r_[np.nan,
                      0.5*(np.angle(-h[2:]*np.conj(h[0:-2]))+np.pi)/(2*np.pi) * samplerate,
                      np.nan]
        print("Mean Frequ of mode ", m, np.mean(np.diff(np.unwrap(phase[:,np.r_[0,0:len(modes[0])]]))/(2*np.pi)*samplerate))

        #f(m,:) = [nan 0.5*(angle(-h(t+1).*conj(h(t-1)))+pi)/(2*pi) * sr nan];

    # calc the freqs (old way)
    #f=np.diff(np.unwrap(phase[:,np.r_[0,0:len(modes[0])]]))/(2*np.pi)*samplerate

    # clip the freqs so they don't go below zero
    #f = f.clip(0,f.max())

    return f,amp,phase 
Example 21
Project: pygriffinlim   Author: rbarghou   File: griffinlim.py    License: GNU General Public License v3.0 5 votes vote down vote up
def griffin_lim_generator(
        spectrogram,
        iterations=10,
        approximated_signal=None,
        stft_kwargs={},
        istft_kwargs={}):
    """
    Implements the basic Griffin Lim algorithm.

    Returns a generator that outputs the approximated signal at the various iterations.
    :param spectrogram: The Spectrogram from which reconstruction should begin
    :param iterations: The number of iterations you want to perform reconstruction with.
    :param approximated_signal: if you want to begin with an existing approximated signal
    :param stft_kwargs: The arguments to pass to STFT as defined by the librosa
    implementation of these functions.
    :param istft_kwargs: The arguments to pass to ISTFT as defined by the librosa
    implementation of these functions.
    :return generator:  This is a generator of approximated signals at each iteration of
    the griffin lim algorithm
    """
    _M = spectrogram
    for k in range(iterations):
        if approximated_signal is None:
            _P = np.random.randn(*_M.shape)
        else:
            _D = librosa.stft(approximated_signal, **stft_kwargs)
            _P = np.angle(_D)

        _D = _M * np.exp(1j * _P)
        approximated_signal = librosa.istft(_D, **istft_kwargs)
        yield approximated_signal 
Example 22
Project: pygriffinlim   Author: rbarghou   File: fastgriffinlim.py    License: GNU General Public License v3.0 5 votes vote down vote up
def fast_griffin_lim_generator(
        spectrogram,
        iterations=10,
        approximated_signal=None,
        alpha=.1,
        stft_kwargs={},
        istft_kwargs={}):
    """

    :param spectrogram:
    :param iterations:
    :param approximated_signal:
    :param alpha:
    :param stft_kwargs:
    :param istft_kwargs:
    :return:
    """
    _M = spectrogram
    for k in range(iterations):
        if approximated_signal is None:
            _P = np.random.randn(*_M.shape)
        else:
            _D = librosa.stft(approximated_signal, **stft_kwargs)
            _P = np.angle(_D)

        _D = _M * np.exp(1j * _P)
        _M = spectrogram + (alpha * np.abs(_D))
        approximated_signal = librosa.istft(_D, **istft_kwargs)
        yield approximated_signal 
Example 23
Project: pygriffinlim   Author: rbarghou   File: modifiedfastgriffinlim.py    License: GNU General Public License v3.0 5 votes vote down vote up
def modified_fast_griffin_lim_generator(
        spectrogram,
        iterations,
        approximated_signal=None,
        alpha_loc=.1,
        alpha_scale=.4,
        stft_kwargs={},
        istft_kwargs={}):
    """

    :param spectrogram:
    :param iterations:
    :param approximated_signal:
    :param alpha_loc:
    :param alpha_scale:
    :param stft_kwargs:
    :param istft_kwargs:
    :return:
    """
    _M = spectrogram
    for k in range(iterations):
        if approximated_signal is None:
            _P = np.random.randn(*_M.shape)
        else:
            _D = librosa.stft(approximated_signal, **stft_kwargs)
            _P = np.angle(_D)

        _D = _M * np.exp(1j * _P)
        alpha = np.random.normal(alpha_loc, alpha_scale)
        _M = spectrogram + (alpha * np.abs(_D))
        approximated_signal = librosa.istft(_D, **istft_kwargs)
        yield approximated_signal 
Example 24
Project: parallel-wavenet-vocoder   Author: andabi   File: audio.py    License: MIT License 5 votes vote down vote up
def wav2spec(wav, n_fft, win_length, hop_length, time_first=True):
    """
    Get magnitude and phase spectrogram from waveforms.

    Parameters
    ----------
    wav : np.ndarray [shape=(n,)]
        The real-valued waveform.

    n_fft : int > 0 [scalar]
        FFT window size.

    win_length  : int <= n_fft [scalar]
        The window will be of length `win_length` and then padded
        with zeros to match `n_fft`.

    hop_length : int > 0 [scalar]
        Number audio of frames between STFT columns.

    time_first : boolean. optional.
        if True, time axis is followed by bin axis. In this case, shape of returns is (t, 1 + n_fft/2)

    Returns
    -------
    mag : np.ndarray [shape=(t, 1 + n_fft/2) or (1 + n_fft/2, t)]
        Magnitude spectrogram.

    phase : np.ndarray [shape=(t, 1 + n_fft/2) or (1 + n_fft/2, t)]
        Phase spectrogram.

    """
    stft = librosa.stft(y=wav, n_fft=n_fft, hop_length=hop_length, win_length=win_length)
    mag = np.abs(stft)
    phase = np.angle(stft)

    if time_first:
        mag = mag.T
        phase = phase.T

    return mag, phase 
Example 25
Project: sprocket   Author: k2kobayashi   File: ms.py    License: MIT License 5 votes vote down vote up
def logpowerspec(self, data, phase=False):
        """Calculate log power spectrum in each dimension

        Parameters
        ----------
        data : array, shape (`T`, `dim`)
            Array of data sequence

        Returns
        -------
        logpowerspec : array, shape (`dim`, `fftsize // 2 + 1`)
            Log power spectrum of sequence
        phasespec : array, shape (`dim`, `fftsize // 2 + 1`), optional
            Phase spectrum of sequences
        """

        # create zero padded data
        T, dim = data.shape
        padded_data = np.zeros((self.fftsize, dim))
        padded_data[:T] += data

        # calculate log power spectum of data
        complex_spec = np.fft.fftn(padded_data, axes=(0, 1))
        logpowerspec = 2 * np.log(np.absolute(complex_spec))

        if phase:
            phasespec = np.angle(complex_spec)
            return logpowerspec, phasespec
        else:
            return logpowerspec 
Example 26
Project: ludwig   Author: uber   File: audio_utils.py    License: Apache License 2.0 5 votes vote down vote up
def get_phase_stft_magnitude(raw_data, sampling_rate_in_hz, window_length_in_s,
                             window_shift_in_s, num_fft_points, window_type):
    stft = _get_stft(raw_data, sampling_rate_in_hz, window_length_in_s,
                     window_shift_in_s, num_fft_points, window_type=window_type)
    abs_stft = np.abs(stft)
    phase = np.angle(stft)
    stft_phase = np.concatenate((phase, abs_stft), axis=1)
    return np.transpose(stft_phase) 
Example 27
Project: tenpy   Author: tenpy   File: site.py    License: GNU General Public License v3.0 5 votes vote down vote up
def rename_op(self, old_name, new_name):
        """Rename an added operator.

        Parameters
        ----------
        old_name : str
            The old name of the operator.
        new_name : str
            The new name of the operator.
        """
        if old_name == new_name:
            return
        if new_name in self.opnames:
            raise ValueError("new_name already exists")
        old_hc_name = self.hc_ops.get(old_name, None)
        op = getattr(self, old_name)
        need_JW = old_name in self.need_JW_string
        hc_op_name = self.get_hc_op_name(old_name)
        self.remove_op(old_name)
        setattr(self, new_name, op)
        self.opnames.add(new_name)
        if need_JW:
            self.need_JW_string.add(new_name)
        if new_name == 'JW':
            self.JW_exponent = np.real_if_close(np.angle(np.diag(op.to_ndarray())) / np.pi)
        if old_hc_name is not None:
            if old_hc_name == old_name:
                self.hc_ops[new_name] = new_name
            else:
                self.hc_ops[new_name] = old_hc_name
                self.hc_ops[old_hc_name] = new_name 
Example 28
Project: pylops   Author: equinor   File: sparsity.py    License: GNU Lesser General Public License v3.0 5 votes vote down vote up
def _softthreshold(x, thresh):
    r"""Soft thresholding.

    Applies soft thresholding to vector ``x`` (equal to the proximity
    operator for :math:`||\mathbf{x}||_1`) as shown in [1]_.

    .. [1] Chen, Y., Chen, K., Shi, P., Wang, Y., “Irregular seismic
       data reconstruction using a percentile-half-thresholding algorithm”,
       Journal of Geophysics and Engineering, vol. 11. 2014.

    Parameters
    ----------
    x : :obj:`numpy.ndarray`
        Vector
    thresh : :obj:`float`
        Threshold

    Returns
    -------
    x1 : :obj:`numpy.ndarray`
        Tresholded vector

    """
    #if np.iscomplexobj(x):
    #    # https://stats.stackexchange.com/questions/357339/soft-thresholding-
    #    # for-the-lasso-with-complex-valued-data
    #    x1 = np.maximum(np.abs(x) - thresh, 0.) * np.exp(1j * np.angle(x))
    #else:
    #    x1 = np.maximum(np.abs(x)-thresh, 0.) * np.sign(x)
    x1 = x - thresh * x / np.abs(x)
    x1[np.abs(x) <= thresh] = 0
    return x1 
Example 29
Project: tf-pose   Author: SrikanthVelpuri   File: functions.py    License: Apache License 2.0 5 votes vote down vote up
def removePeriodic(data, f0=60.0, dt=None, harmonics=10, samples=4):
    if (hasattr(data, 'implements') and data.implements('MetaArray')):
        data1 = data.asarray()
        if dt is None:
            times = data.xvals('Time')
            dt = times[1]-times[0]
    else:
        data1 = data
        if dt is None:
            raise Exception('Must specify dt for this data')
    
    ft = np.fft.fft(data1)
    
    ## determine frequencies in fft data
    df = 1.0 / (len(data1) * dt)
    freqs = np.linspace(0.0, (len(ft)-1) * df, len(ft))
    
    ## flatten spikes at f0 and harmonics
    for i in xrange(1, harmonics + 2):
        f = f0 * i # target frequency
        
        ## determine index range to check for this frequency
        ind1 = int(np.floor(f / df))
        ind2 = int(np.ceil(f / df)) + (samples-1)
        if ind1 > len(ft)/2.:
            break
        mag = (abs(ft[ind1-1]) + abs(ft[ind2+1])) * 0.5
        for j in range(ind1, ind2+1):
            phase = np.angle(ft[j])   ## Must preserve the phase of each point, otherwise any transients in the trace might lead to large artifacts.
            re = mag * np.cos(phase)
            im = mag * np.sin(phase)
            ft[j] = re + im*1j
            ft[len(ft)-j] = re - im*1j
            
    data2 = np.fft.ifft(ft).real
    
    if (hasattr(data, 'implements') and data.implements('MetaArray')):
        return metaarray.MetaArray(data2, info=data.infoCopy())
    else:
        return data2 
Example 30
Project: tf-pose   Author: SrikanthVelpuri   File: Filters.py    License: Apache License 2.0 5 votes vote down vote up
def processData(self, data):
        times = data.xvals('Time')
        dt = times[1]-times[0]
        
        data1 = data.asarray()
        ft = np.fft.fft(data1)
        
        ## determine frequencies in fft data
        df = 1.0 / (len(data1) * dt)
        freqs = np.linspace(0.0, (len(ft)-1) * df, len(ft))
        
        ## flatten spikes at f0 and harmonics
        f0 = self.ctrls['f0'].value()
        for i in xrange(1, self.ctrls['harmonics'].value()+2):
            f = f0 * i # target frequency
            
            ## determine index range to check for this frequency
            ind1 = int(np.floor(f / df))
            ind2 = int(np.ceil(f / df)) + (self.ctrls['samples'].value()-1)
            if ind1 > len(ft)/2.:
                break
            mag = (abs(ft[ind1-1]) + abs(ft[ind2+1])) * 0.5
            for j in range(ind1, ind2+1):
                phase = np.angle(ft[j])   ## Must preserve the phase of each point, otherwise any transients in the trace might lead to large artifacts.
                re = mag * np.cos(phase)
                im = mag * np.sin(phase)
                ft[j] = re + im*1j
                ft[len(ft)-j] = re - im*1j
                
        data2 = np.fft.ifft(ft).real
        
        ma = metaarray.MetaArray(data2, info=data.infoCopy())
        return ma