Python fftpack.fft方法代码示例

# 需要导入模块: from scipy import fftpack [as 别名]
# 或者: from scipy.fftpack import fft [as 别名]
def test_high_frequency_completion(self):
        path = dirpath + '/data/test16000.wav'
        fs, x = wavfile.read(path)

        f0rate = 0.5
        shifter = Shifter(fs, f0rate=f0rate)
        mod_x = shifter.f0transform(x, completion=False)
        mod_xc = shifter.f0transform(x, completion=True)
        assert len(mod_x) == len(mod_xc)

        N = 512
        fl = int(fs * 25 / 1000)
        win = np.hanning(fl)
        sts = [1000, 5000, 10000, 20000]
        for st in sts:
            # confirm w/o completion
            f_mod_x = fft(mod_x[st: st + fl] / 2**16 * win)
            amp_mod_x = 20.0 * np.log10(np.abs(f_mod_x))

            # confirm w/ completion
            f_mod_xc = fft(mod_xc[st: st + fl] / 2**16 * win)
            amp_mod_xc = 20.0 * np.log10(np.abs(f_mod_xc))

            assert np.mean(amp_mod_x[N // 4:] < np.mean(amp_mod_xc[N // 4:])) 

# 需要导入模块: from scipy import fftpack [as 别名]
# 或者: from scipy.fftpack import fft [as 别名]
def _centered(arr, newsize):
    # Return the center newsize portion of the array.
    newsize = numpy.asarray(newsize)
    currsize = numpy.array(arr.shape)
    startind = (currsize - newsize) // 2
    endind = startind + newsize
    myslice = [slice(startind[k], endind[k]) for k in range(len(endind))]
    return arr[tuple(myslice)]


#def _rfftn(a, s=None, axes=None):
#    s, axes = _cook_nd_args(a, s, axes);
#    a = fft.rfft(a, s[-1], axes[-1], overwrite_x = True)
#    for ii in range(len(axes)-1):
#        a = fft.fft(a, s[ii], axes[ii], overwrite_x = True)
#    return a

#def _irfftn(a, s=None, axes=None):   
#    #a = asarray(a).astype('complex64')
#    s, axes = _cook_nd_args(a, s, axes, invreal=1)
#    for ii in range(len(axes)-1):
#        a = fft.ifft(a, s[ii], axes[ii], overwrite_x = True);
#    a = fft.ifft(a, s[-1], axes[-1], overwrite_x = True);
#    a = a.real;
#    return a 

# 需要导入模块: from scipy import fftpack [as 别名]
# 或者: from scipy.fftpack import fft [as 别名]
def spectralwhitening(stream):
    """
    Apply spectral whitening to data.
    Data is divided by its smoothed (Default: None) amplitude spectrum.
    """
    stream2 = copy.deepcopy(stream)

    for trace in arange(len(stream2)):
        data = stream2[trace].data

        n = len(data)
        nfft = nextpow2(n)

        spec = fft(data, nfft)
        spec_ampl = sqrt(abs(multiply(spec, conjugate(spec))))

        spec /= spec_ampl  # Do we need to do some smoothing here?
        ret = real(ifft(spec, nfft)[:n])

        stream2[trace].data = ret

    return stream2 

# 需要导入模块: from scipy import fftpack [as 别名]
# 或者: from scipy.fftpack import fft [as 别名]
def spectralwhitening_smooth(stream, N):
    """
    Apply spectral whitening to data.
    Data is divided by its smoothed (Default: None) amplitude spectrum.
    """
    stream2 = copy.deepcopy(stream)

    for trace in arange(len(stream2)):
        data = stream2[trace].data

        n = len(data)
        nfft = nextpow2(n)

        spec = fft(data, nfft)
        spec_ampl = sqrt(abs(multiply(spec, conjugate(spec))))

        spec_ampl = smooth(spec_ampl, N)

        spec /= spec_ampl  # Do we need to do some smoothing here?
        ret = real(ifft(spec, nfft)[:n])

        stream2[trace].data = ret

    return stream2 

# 需要导入模块: from scipy import fftpack [as 别名]
# 或者: from scipy.fftpack import fft [as 别名]
def carr_madan_fraction_fft_call_pricer(N, d_u, d_k, alpha, r, t, S0, q, chf_ln_st):
    rou = (d_u * d_k) / (2 * np.pi)
    beta = np.log(S0) - d_k * N / 2
    u_arr = np.arange(N) * d_u
    k_arr = beta + np.arange(N) * d_k
    delta_arr = np.zeros(N)
    delta_arr[0] = 1
    w_arr = d_u / 3 * (3 + (-1) ** (np.arange(N) + 1) - delta_arr)
    call_chf = (np.exp(-r * t) / ((alpha + 1j * u_arr) * (alpha + 1j * u_arr + 1))) * chf_ln_st(
        u_arr - (alpha + 1) * 1j,
        t, r, q=q, S0=S0)
    x_arr = np.exp(-1j * beta * u_arr) * call_chf * w_arr
    y_arr = np.zeros(2 * N) * 0j
    y_arr[:N] = np.exp(-1j * np.pi * rou * np.arange(N) ** 2) * x_arr
    z_arr = np.zeros(2 * N) * 0j
    z_arr[:N] = np.exp(1j * np.pi * rou * np.arange(N) ** 2)
    z_arr[N:] = np.exp(1j * np.pi * rou * np.arange(N - 1, -1, -1) ** 2)
    ffty = (fft(y_arr))
    fftz = (fft(z_arr))
    fftx = ffty * fftz
    fftpsi = ifft(fftx)
    fft_prices = np.exp(-1j * np.pi * (np.arange(N) ** 2) * rou) * fftpsi[:N]
    call_prices = (np.exp(-alpha * k_arr) / np.pi) * fft_prices.real
    return np.exp(k_arr), call_prices 

# 需要导入模块: from scipy import fftpack [as 别名]
# 或者: from scipy.fftpack import fft [as 别名]
def test_shape_axes_argument(self):
        small_x = [[1,2,3],[4,5,6],[7,8,9]]
        large_x1 = array([[1,2,3,0],
                                  [4,5,6,0],
                                  [7,8,9,0],
                                  [0,0,0,0]])
        # Disable tests with shape and axes of different lengths
        # y = fftn(small_x,shape=(4,4),axes=(-1,))
        # for i in range(4):
        #    assert_array_almost_equal (y[i],fft(large_x1[i]))
        # y = fftn(small_x,shape=(4,4),axes=(-2,))
        # for i in range(4):
        #    assert_array_almost_equal (y[:,i],fft(large_x1[:,i]))
        y = fftn(small_x,shape=(4,4),axes=(-2,-1))
        assert_array_almost_equal(y,fftn(large_x1))
        y = fftn(small_x,shape=(4,4),axes=(-1,-2))
        assert_array_almost_equal(y,swapaxes(
            fftn(swapaxes(large_x1,-1,-2)),-1,-2)) 

# 需要导入模块: from scipy import fftpack [as 别名]
# 或者: from scipy.fftpack import fft [as 别名]
def cut_norm(full_y, dt, area=1.0):
    """Cut out an FFT spectrum from scipy.fftpack.fft() (or numpy.fft.fft())
    result and normalize the integral `int(f) y(f) df = area`.

    full_y : 1d array
        Result of fft(...)
    dt : time step
    area : integral area
    """
    full_faxis = np.fft.fftfreq(full_y.shape[0], dt)
    split_idx = full_faxis.shape[0]/2
    y_out = full_y[:split_idx]
    faxis = full_faxis[:split_idx]
    return faxis, num.norm_int(y_out, faxis, area=area)


###############################################################################
# common settings for 1d and 3d case
############################################################################### 

# 需要导入模块: from scipy import fftpack [as 别名]
# 或者: from scipy.fftpack import fft [as 别名]
def ezfft(y, dt=1.0):
    """Simple FFT function for interactive use.

    Parameters
    ----------
    y : 1d array to fft
    dt : float
        time step

    Returns
    -------
    faxis, fft(y)

    Examples
    --------
    >>> t = linspace(0,1,200)
    >>> x = sin(2*pi*10*t) + sin(2*pi*20*t)
    >>> f,d = signal.ezfft(x, dt=t[1]-t[0])
    >>> plot(f,abs(d))
    """
    assert y.ndim == 1
    faxis = np.fft.fftfreq(len(y), dt)
    split_idx = len(faxis)/2
    return faxis[:split_idx], fft(y)[:split_idx] 

# 需要导入模块: from scipy import fftpack [as 别名]
# 或者: from scipy.fftpack import fft [as 别名]
def fft_1d_loop(arr, axis=-1):
    """Like scipy.fft.pack.fft and numpy.fft.fft, perform fft along an axis.
    Here do this by looping over remaining axes and perform 1D FFTs.

    This was implemented as a low-memory version like
    :func:`~pwtools.crys.smooth` to be used in :func:`~pwtools.pydos.pdos`,
    which fills up the memory for big MD data. But actually it has the same
    memory footprint as the plain scipy fft routine. Keep it here anyway as a
    nice reference for how to loop over remaining axes in the ndarray case.
    """
    if axis < 0:
        axis = arr.ndim - 1
    axes = [ax for ax in range(arr.ndim) if ax != axis]
    # tuple here is 3x faster than generator expression
    #   idxs = (range(arr.shape[ax]) for ax in axes)
    idxs = tuple(range(arr.shape[ax]) for ax in axes)
    out = np.empty(arr.shape, dtype=complex)
    for idx_tup in product(*idxs):
        sl = [slice(None)] * arr.ndim
        for idx,ax in zip(idx_tup, axes):
            sl[ax] = idx
        tsl = tuple(sl)
        out[tsl] = fft(arr[tsl])
    return out 

# 需要导入模块: from scipy import fftpack [as 别名]
# 或者: from scipy.fftpack import fft [as 别名]
def stft(X, fftsize=128, step="half", mean_normalize=True, real=False,
         compute_onesided=True):
    """
    Compute STFT for 1D real valued input X
    """
    if real:
        local_fft = fftpack.rfft
        cut = -1
    else:
        local_fft = fftpack.fft
        cut = None
    if compute_onesided:
        cut = fftsize // 2 + 1
    if mean_normalize:
        X -= X.mean()
    if step == "half":
        X = halfoverlap(X, fftsize)
    else:
        X = overlap(X, fftsize, step)
    size = fftsize
    win = 0.54 - .46 * np.cos(2 * np.pi * np.arange(size) / (size - 1))
    X = X * win[None]
    X = local_fft(X)[:, :cut]
    return X 

# 需要导入模块: from scipy import fftpack [as 别名]
# 或者: from scipy.fftpack import fft [as 别名]
def cheaptrick_get_power_spectrum(waveform, fs, fft_size, f0):
    power_spectrum = np.abs(np.fft.fft(waveform, fft_size)) ** 2
    frequency_axis = np.arange(fft_size) / float(fft_size) * float(fs)
    ind = frequency_axis < (f0 + fs / fft_size)
    low_frequency_axis = frequency_axis[ind]
    low_frequency_replica = interp1d(f0 - low_frequency_axis,
            power_spectrum[ind], kind="linear",
            fill_value="extrapolate")(low_frequency_axis)
    p1 = low_frequency_replica[(frequency_axis < f0)[:len(low_frequency_replica)]]
    p2 = power_spectrum[(frequency_axis < f0)[:len(power_spectrum)]]
    power_spectrum[frequency_axis < f0] = p1 + p2
    lb1 = int(fft_size / 2) + 1
    lb2 = 1
    ub2 = int(fft_size / 2)
    power_spectrum[lb1:] = power_spectrum[lb2:ub2][::-1]
    return power_spectrum 

# 需要导入模块: from scipy import fftpack [as 别名]
# 或者: from scipy.fftpack import fft [as 别名]
def d4c_love_train(x, fs, current_f0, current_position, threshold):
    vuv = 0
    if current_f0 == 0:
        return vuv
    lowest_f0 = 40
    current_f0 = max([current_f0, lowest_f0])
    fft_size = int(2 ** np.ceil(np.log2(3. * fs / lowest_f0 + 1)))
    boundary0 = int(np.ceil(100 / (float(fs) / fft_size)))
    boundary1 = int(np.ceil(4000 / (float(fs) / fft_size)))
    boundary2 = int(np.ceil(7900 / (float(fs) / fft_size)))

    waveform = d4c_get_windowed_waveform(x, fs, current_f0, current_position,
            1.5, 2)
    power_spectrum = np.abs(np.fft.fft(waveform, int(fft_size)) ** 2)
    power_spectrum[0:boundary0 + 1] = 0.
    cumulative_spectrum = np.cumsum(power_spectrum)
    if (cumulative_spectrum[boundary1] / cumulative_spectrum[boundary2]) > threshold:
        vuv = 1
    return vuv 

# 需要导入模块: from scipy import fftpack [as 别名]
# 或者: from scipy.fftpack import fft [as 别名]
def win2mgc(windowed_signal, order=20, alpha=0.35, gamma=-0.41, miniter=2,
        maxiter=30, criteria=0.001, otype=0, verbose=False):
    """
    Accepts 1D or 2D array of windowed signal frames.

    If 2D, assumes time is axis 0.

    Returns mel generalized cepstral coefficients.

    Based on r9y9 Julia code
    https://github.com/r9y9/MelGeneralizedCepstrums.jl
    """
    if len(windowed_signal.shape) == 1:
        sp = np.fft.fft(windowed_signal)
        return _sp2mgc(sp, order=order, alpha=alpha, gamma=gamma,
                miniter=miniter, maxiter=maxiter, criteria=criteria,
                otype=otype, verbose=verbose)
    else:
        raise ValueError("2D input not yet complete for win2mgc") 

# 需要导入模块: from scipy import fftpack [as 别名]
# 或者: from scipy.fftpack import fft [as 别名]
def run_mgc_example():
    import matplotlib.pyplot as plt
    fs, x = wavfile.read("test16k.wav")
    pos = 3000
    fftlen = 1024
    win = np.blackman(fftlen) / np.sqrt(np.sum(np.blackman(fftlen) ** 2))
    xw = x[pos:pos + fftlen] * win
    sp = 20 * np.log10(np.abs(np.fft.rfft(xw)))
    mgc_order = 20
    mgc_alpha = 0.41
    mgc_gamma = -0.35
    mgc_arr = win2mgc(xw, order=mgc_order, alpha=mgc_alpha, gamma=mgc_gamma, verbose=True)
    xwsp = 20 * np.log10(np.abs(np.fft.rfft(xw)))
    sp = mgc2sp(mgc_arr, mgc_alpha, mgc_gamma, fftlen)
    plt.plot(xwsp)
    plt.plot(20. / np.log(10) * np.real(sp), "r")
    plt.xlim(1, len(xwsp))
    plt.show() 

# 需要导入模块: from scipy import fftpack [as 别名]
# 或者: from scipy.fftpack import fft [as 别名]
def calculate_falff(timeseries, min_low_freq, max_low_freq, min_total_freq, max_total_freq, calc_alff):
    ''' this will calculate falff from a timeseries'''

    n = len(timeseries)
    time = (np.arange(n))*2

    # Takes fast Fourier transform of timeseries
    fft_timeseries = fft(timeseries)
    # Calculates frequency scale
    freq_scale = np.fft.fftfreq(n, 1/1)

    # Calculates power of fft
    mag = (abs(fft_timeseries))**0.5

    # Finds low frequency range (0.01-0.08) and total frequency range (0.0-0.25)
    low_ind = np.where((float(min_low_freq) <= freq_scale) & (freq_scale <= float(max_low_freq)))
    total_ind = np.where((float(min_total_freq) <= freq_scale) & (freq_scale <= float(max_total_freq)))

    # Indexes power to low frequency index, total frequency range
    low_power = mag[low_ind]
    total_power = mag[total_ind]
    # Calculates sum of lower power and total power
    low_pow_sum = np.sum(low_power)
    total_pow_sum = np.sum(total_power)

    # Calculates alff as the sum of amplitudes within the low frequency range
    if calc_alff:
        calc = low_pow_sum
    # Calculates falff as the sum of power in low frequnecy range divided by sum of power in the total frequency range
    else:
        calc = np.divide(low_pow_sum, total_pow_sum)

    return calc