Python fft.rfft方法代码示例

# 需要导入模块: from numpy import fft [as 别名]
# 或者: from numpy.fft import rfft [as 别名]
def pink(N):
    """
    Pink noise.

    :param N: Amount of samples.

    Pink noise has equal power in bands that are proportionally wide.
    Power density decreases with 3 dB per octave.

    """
    # This method uses the filter with the following coefficients.
    # b = np.array([0.049922035, -0.095993537, 0.050612699, -0.004408786])
    # a = np.array([1, -2.494956002, 2.017265875, -0.522189400])
    # return lfilter(B, A, np.random.randn(N))
    # Another way would be using the FFT
    # x = np.random.randn(N)
    # X = rfft(x) / N
    uneven = N % 2
    X = np.random.randn(N // 2 + 1 + uneven) + 1j * np.random.randn(N // 2 + 1 + uneven)
    S = np.sqrt(np.arange(len(X)) + 1.)  # +1 to avoid divide by zero
    y = (irfft(X / S)).real
    if uneven:
        y = y[:-1]
    return normalize(y) 

# 需要导入模块: from numpy import fft [as 别名]
# 或者: from numpy.fft import rfft [as 别名]
def mass(query, ts):
    """
    >>> np.round(mass(np.array([0.0, 1.0, -1.0, 0.0]), np.array([-1, 1, 0, 0, -1, 1])), 3)
    array([ 2.   ,  2.828,  2.   ])
    """
    query = zNormalize(query)
    m = len(query)
    n = len(ts)

    stdv = movstd(ts, m)
    query = query[::-1]
    query = np.pad(query, (0, n - m), 'constant')
    dots = fft.irfft(fft.rfft(ts) * fft.rfft(query))
    return np.sqrt(2 * (m - (dots[m - 1 :] / stdv))) 

# 需要导入模块: from numpy import fft [as 别名]
# 或者: from numpy.fft import rfft [as 别名]
def bench_random(self):
        from numpy.fft import rfft as numpy_rfft
        print()
        print('Fast Fourier Transform (real data)')
        print('==================================')
        print(' size |  scipy  |  numpy  ')
        print('----------------------------------')
        for size,repeat in [(100,7000),(1000,2000),
                            (256,10000),
                            (512,10000),
                            (1024,1000),
                            (2048,1000),
                            (2048*2,500),
                            (2048*4,500),
                            ]:
            print('%5s' % size, end=' ')
            sys.stdout.flush()

            x = random([size]).astype(double)
            print('|%8.2f' % measure('rfft(x)',repeat), end=' ')
            sys.stdout.flush()

            print('|%8.2f' % measure('numpy_rfft(x)',repeat), end=' ')
            sys.stdout.flush()

            print(' (secs for %s calls)' % (repeat))
        sys.stdout.flush() 

# 需要导入模块: from numpy import fft [as 别名]
# 或者: from numpy.fft import rfft [as 别名]
def test_fourier_gaussian_real01(self):
        for shape in [(32, 16), (31, 15)]:
            for type in [numpy.float32, numpy.float64]:
                a = numpy.zeros(shape, type)
                a[0, 0] = 1.0
                a = fft.rfft(a, shape[0], 0)
                a = fft.fft(a, shape[1], 1)
                a = ndimage.fourier_gaussian(a, [5.0, 2.5],
                                                       shape[0], 0)
                a = fft.ifft(a, shape[1], 1)
                a = fft.irfft(a, shape[0], 0)
                assert_almost_equal(ndimage.sum(a), 1) 

# 需要导入模块: from numpy import fft [as 别名]
# 或者: from numpy.fft import rfft [as 别名]
def test_fourier_uniform_real01(self):
        for shape in [(32, 16), (31, 15)]:
            for type in [numpy.float32, numpy.float64]:
                a = numpy.zeros(shape, type)
                a[0, 0] = 1.0
                a = fft.rfft(a, shape[0], 0)
                a = fft.fft(a, shape[1], 1)
                a = ndimage.fourier_uniform(a, [5.0, 2.5],
                                                      shape[0], 0)
                a = fft.ifft(a, shape[1], 1)
                a = fft.irfft(a, shape[0], 0)
                assert_almost_equal(ndimage.sum(a), 1.0) 

# 需要导入模块: from numpy import fft [as 别名]
# 或者: from numpy.fft import rfft [as 别名]
def test_fourier_shift_real01(self):
        for shape in [(32, 16), (31, 15)]:
            for dtype in [numpy.float32, numpy.float64]:
                expected = numpy.arange(shape[0] * shape[1], dtype=dtype)
                expected.shape = shape
                a = fft.rfft(expected, shape[0], 0)
                a = fft.fft(a, shape[1], 1)
                a = ndimage.fourier_shift(a, [1, 1], shape[0], 0)
                a = fft.ifft(a, shape[1], 1)
                a = fft.irfft(a, shape[0], 0)
                assert_array_almost_equal(a[1:, 1:], expected[:-1, :-1])
                assert_array_almost_equal(a.imag, numpy.zeros(shape)) 

# 需要导入模块: from numpy import fft [as 别名]
# 或者: from numpy.fft import rfft [as 别名]
def test_fourier_ellipsoid_real01(self):
        for shape in [(32, 16), (31, 15)]:
            for type in [numpy.float32, numpy.float64]:
                a = numpy.zeros(shape, type)
                a[0, 0] = 1.0
                a = fft.rfft(a, shape[0], 0)
                a = fft.fft(a, shape[1], 1)
                a = ndimage.fourier_ellipsoid(a, [5.0, 2.5],
                                              shape[0], 0)
                a = fft.ifft(a, shape[1], 1)
                a = fft.irfft(a, shape[0], 0)
                assert_almost_equal(ndimage.sum(a), 1.0) 

# 需要导入模块: from numpy import fft [as 别名]
# 或者: from numpy.fft import rfft [as 别名]
def test_fourier_gaussian_real01(self):
        for shape in [(32, 16), (31, 15)]:
            for type_, dec in zip([numpy.float32, numpy.float64], [6, 14]):
                a = numpy.zeros(shape, type_)
                a[0, 0] = 1.0
                a = fft.rfft(a, shape[0], 0)
                a = fft.fft(a, shape[1], 1)
                a = ndimage.fourier_gaussian(a, [5.0, 2.5], shape[0], 0)
                a = fft.ifft(a, shape[1], 1)
                a = fft.irfft(a, shape[0], 0)
                assert_almost_equal(ndimage.sum(a), 1, decimal=dec) 

# 需要导入模块: from numpy import fft [as 别名]
# 或者: from numpy.fft import rfft [as 别名]
def test_fourier_uniform_real01(self):
        for shape in [(32, 16), (31, 15)]:
            for type_, dec in zip([numpy.float32, numpy.float64], [6, 14]):
                a = numpy.zeros(shape, type_)
                a[0, 0] = 1.0
                a = fft.rfft(a, shape[0], 0)
                a = fft.fft(a, shape[1], 1)
                a = ndimage.fourier_uniform(a, [5.0, 2.5], shape[0], 0)
                a = fft.ifft(a, shape[1], 1)
                a = fft.irfft(a, shape[0], 0)
                assert_almost_equal(ndimage.sum(a), 1.0, decimal=dec) 

# 需要导入模块: from numpy import fft [as 别名]
# 或者: from numpy.fft import rfft [as 别名]
def test_fourier_shift_real01(self):
        for shape in [(32, 16), (31, 15)]:
            for type_, dec in zip([numpy.float32, numpy.float64], [4, 11]):
                expected = numpy.arange(shape[0] * shape[1], dtype=type_)
                expected.shape = shape
                a = fft.rfft(expected, shape[0], 0)
                a = fft.fft(a, shape[1], 1)
                a = ndimage.fourier_shift(a, [1, 1], shape[0], 0)
                a = fft.ifft(a, shape[1], 1)
                a = fft.irfft(a, shape[0], 0)
                assert_array_almost_equal(a[1:, 1:], expected[:-1, :-1],
                                          decimal=dec)
                assert_array_almost_equal(a.imag, numpy.zeros(shape),
                                          decimal=dec) 

# 需要导入模块: from numpy import fft [as 别名]
# 或者: from numpy.fft import rfft [as 别名]
def test_fourier_ellipsoid_real01(self):
        for shape in [(32, 16), (31, 15)]:
            for type_, dec in zip([numpy.float32, numpy.float64], [5, 14]):
                a = numpy.zeros(shape, type_)
                a[0, 0] = 1.0
                a = fft.rfft(a, shape[0], 0)
                a = fft.fft(a, shape[1], 1)
                a = ndimage.fourier_ellipsoid(a, [5.0, 2.5],
                                              shape[0], 0)
                a = fft.ifft(a, shape[1], 1)
                a = fft.irfft(a, shape[0], 0)
                assert_almost_equal(ndimage.sum(a), 1.0, decimal=dec) 

# 需要导入模块: from numpy import fft [as 别名]
# 或者: from numpy.fft import rfft [as 别名]
def PowerSpectrum(X, Fs):
  """
  Calculate the power spectrum of real-valued signal X measured with a
  sampling frequency Fs. Result is in the same units as Fs.
  """

  X = X - X.mean()
  F = rfft(X)[:-1]
  freq = np.arange(0, 0.5*Fs, Fs*1.0/len(X))
  pw = np.abs(F)**2
  return freq, pw/max(pw) 

# 需要导入模块: from numpy import fft [as 别名]
# 或者: from numpy.fft import rfft [as 别名]
def lpc_python(x, order):
    def levinson(R,order):
        """ input: autocorrelation and order, output: LPC coefficients """
        a   = np.zeros(order+2)
        a[0] = -1
        k = np.zeros(order+1)
        # step 1: initialize prediction error "e" to R[0]
        e = R[0]
        # step 2: iterate over [1:order]
        for i in range(1,order+1):
            # step 2-1: calculate PARCOR coefficients
            k[i]= (R[i] - np.sum(a[1:i] * R[i-1:0:-1])) / e
            # step 2-2: update LPCs
            a[i] = np.copy(k[i])
            a_old = np.copy(a) 
            for j in range(1,i):
                a[j] -=  k[i]* a_old[i-j] 
            # step 2-3: update prediction error "e" 
            e = e * (1.0 - k[i]**2)
        return -1*a[0:order+1], e, -1*k[1:]
    # N: compute next power of 2
    n = x.shape[0]
    N = int(np.power(2, np.ceil(np.log2(2 * n - 1))))
    # calculate autocorrelation using FFT
    X = rfft(x, N)
    r = irfft(abs(X) ** 2)[:order+1]
    return levinson(r, order) 

# 需要导入模块: from numpy import fft [as 别名]
# 或者: from numpy.fft import rfft [as 别名]
def multS(s, v, L, TP=False):

    N = s.shape[0]
    vp = prepare_v(v, N, L, TP=TP)
    p = irfft(rfft(vp) * rfft(s))
    if not TP:
        return p[:L]
    return p[L - 1:] 

# 需要导入模块: from numpy import fft [as 别名]
# 或者: from numpy.fft import rfft [as 别名]
def autocorr(x):
    """ Fast autocorrelation computation using the FFT """
    
    X = fft.rfft(x, n=2*x.shape[0])
    r = fft.irfft(np.abs(X)**2)
    
    return r[:x.shape[0]] / (x.shape[0] - 1) 

# 需要导入模块: from numpy import fft [as 别名]
# 或者: from numpy.fft import rfft [as 别名]
def toeplitz_multiplication(c, r, A, **kwargs):
    """ 
    Compute numpy.dot(scipy.linalg.toeplitz(c,r), A) using the FFT.

    Parameters
    ----------
    c: ndarray
        the first column of the Toeplitz matrix
    r: ndarray
        the first row of the Toeplitz matrix
    A: ndarray
        the matrix to multiply on the right
    """
    
    m = c.shape[0]
    n = r.shape[0]

    fft_len = int(2**np.ceil(np.log2(m+n-1)))
    zp = fft_len - m - n + 1
    
    if A.shape[0] != n:
        raise ValueError('A dimensions not compatible with toeplitz(c,r)')
    
    x = np.concatenate((c, np.zeros(zp, dtype=c.dtype), r[-1:0:-1]))
    xf = np.fft.rfft(x, n=fft_len)
    
    Af = np.fft.rfft(A, n=fft_len, axis=0)
    
    return np.fft.irfft((Af.T*xf).T, n=fft_len, axis=0)[:m,] 

# 需要导入模块: from numpy import fft [as 别名]
# 或者: from numpy.fft import rfft [as 别名]
def analysis(self, x):
        """
        Perform frequency analysis of a real input using DFT.

        Parameters
        ----------
        x : numpy array
            Real signal in time domain.

        Returns
        -------
        numpy array
            DFT of input.
        """

        # check for valid input
        if x.shape != self.x.shape:
            raise ValueError('Invalid input dimensions! Got (%d, %d), expecting (%d, %d).' 
                % (x.shape[0], x.shape[1],
                self.x.shape[0], self.x.shape[1]))

        # apply window if needed
        if self.analysis_window is not None:
            np.multiply(self.analysis_window, x, x)

        # apply DFT
        if self.transform == 'fftw':
            self.a[:,] = x
            self.X[:,] = self._forward()
        elif self.transform == 'mkl':
            self.X[:,] = mkl_fft.rfft_numpy(x, self.nfft, axis=self.axis)
        else:
            self.X[:,] = rfft(x, self.nfft, axis=self.axis)

        return self.X 

# 需要导入模块: from numpy import fft [as 别名]
# 或者: from numpy.fft import rfft [as 别名]
def deconvolve(y, s, length=None, thresh=0.0):
    '''
    Deconvolve an excitation signal from an impulse response

    Parameters
    ----------

    y : ndarray
        The recording
    s : ndarray
        The excitation signal
    length: int, optional
        the length of the impulse response to deconvolve
    thresh : float, optional
        ignore frequency bins with power lower than this
    '''

    # FFT length including zero padding
    n = y.shape[0] + s.shape[0] - 1

    # let FFT size be even for convenience
    if n % 2 != 0:
        n += 1

    # when unknown, pick the filter size as size of test signal
    if length is None:
        length = n

    # Forward transforms
    Y  = fft.rfft(np.array(y, dtype=np.float32), n=n) / np.sqrt(n)
    S = fft.rfft(np.array(s, dtype=np.float32), n=n) / np.sqrt(n)

    # Only do the division where S is large enough
    H = np.zeros(*Y.shape, dtype=Y.dtype)
    I = np.where(np.abs(S) > thresh)
    H[I] = Y[I] / S[I]

    # Inverse transform
    h = fft.irfft(H, n=n)

    return h[:length] 

# 需要导入模块: from numpy import fft [as 别名]
# 或者: from numpy.fft import rfft [as 别名]
def correlate(x1, x2, interp=1, phat=False):
    '''
    Compute the cross-correlation between x1 and x2

    Parameters
    ----------
    x1,x2: array_like
        The data arrays
    interp: int, optional
        The interpolation factor for the output array, default 1.
    phat: bool, optional
        Apply the PHAT weighting (default False)

    Returns
    -------
    The cross-correlation between the two arrays
    '''

    N1 = x1.shape[0]
    N2 = x2.shape[0]

    N = N1 + N2 - 1

    X1 = fft.rfft(x1, n=N)
    X2 = fft.rfft(x2, n=N)

    if phat:
        eps1 = np.mean(np.abs(X1)) * 1e-10
        X1 /= (np.abs(X1) + eps1)
        eps2 = np.mean(np.abs(X2)) * 1e-10
        X2 /= (np.abs(X2) + eps2)

    m = np.minimum(N1, N2)

    out = fft.irfft(X1*np.conj(X2), n=int(N*interp))

    return np.concatenate([out[-interp*(N2-1):], out[:(interp*N1)]]) 

# 需要导入模块: from numpy import fft [as 别名]
# 或者: from numpy.fft import rfft [as 别名]
def slidingDotProduct(query,ts):
    """
    Calculate the dot product between a query and all subsequences of length(query) in the timeseries ts. Note that we use Numpy's rfft method instead of fft.

    Parameters
    ----------
    query: Specific time series query to evaluate.
    ts: Time series to calculate the query's sliding dot product against.
    """

    m = len(query)
    n = len(ts)


    #If length is odd, zero-pad time time series
    ts_add = 0
    if n%2 ==1:
        ts = np.insert(ts,0,0)
        ts_add = 1

    q_add = 0
    #If length is odd, zero-pad query
    if m%2 == 1:
        query = np.insert(query,0,0)
        q_add = 1

    #This reverses the array
    query = query[::-1]


    query = np.pad(query,(0,n-m+ts_add-q_add),'constant')

    #Determine trim length for dot product. Note that zero-padding of the query has no effect on array length, which is solely determined by the longest vector
    trim = m-1+ts_add

    dot_product = fft.irfft(fft.rfft(ts)*fft.rfft(query))


    #Note that we only care about the dot product results from index m-1 onwards, as the first few values aren't true dot products (due to the way the FFT works for dot products)
    return dot_product[trim :] 

# 需要导入模块: from numpy import fft [as 别名]
# 或者: from numpy.fft import rfft [as 别名]
def processFrames(self, windowedSamples):
        self.complexMixtureSpectrogram[:] = rfft(windowedSamples * self.windowFunction, axis=1).astype(np.complex64)
        self.spectrogram.set_value(self.complexMixtureSpectrogram)
        #self.spectrogram.set_value( rfft(windowedSamples * self.windowFunction, axis=1).astype(np.complex64) )
        
        realGCC = self.getComplexGCC()[0].real
        if self.separationEnabled:
            [inputMask, coefficientMask] = self.getTFMask(realGCC)
            outputSpectrogram = inputMask * self.complexMixtureSpectrogram
            
            if self.coefficientMaskHistories:
                self.coefficientMaskHistories[self.dictionarySize].set(1-coefficientMask)
        else:
            outputSpectrogram = self.complexMixtureSpectrogram.copy()
        
        if self.inputSpectrogramHistory:
            self.inputSpectrogramHistory.set( -np.mean(np.abs(self.complexMixtureSpectrogram), axis=0) ** (1/3.0) )
        if self.gccPHATHistory:
            self.gccPHATHistory.set( np.nanmean(realGCC, axis=0).T )
        if self.tdoaHistory:
            if self.localizationEnabled:
                gccPHATHistory = self.gccPHATHistory.getUnraveledArray()
                tdoaIndex = np.argmax( np.nanmean(gccPHATHistory[:, -self.localizationWindowSize:], axis=-1) )
                #tdoaIndex = (self.targetTDOAIndex.get_value() + 1) % self.numTDOAs
                #tdoaIndex = np.random.randint(0, self.numTDOAs+1)
                self.targetTDOAIndex.set_value(tdoaIndex)
            self.tdoaHistory.set( np.array( [[self.targetTDOAIndex.get_value()]] ) )
        if self.outputSpectrogramHistory:
            self.outputSpectrogramHistory.set( -np.nanmean(np.abs(outputSpectrogram), axis=0) ** (1/3.0) )
        
        return np.fft.irfft(outputSpectrogram, axis=1) * self.synthesisWindowFunction 

# 需要导入模块: from numpy import fft [as 别名]
# 或者: from numpy.fft import rfft [as 别名]
def mkl_toeplitz_multiplication(c, r, A, A_padded=False, out=None, fft_len=None):
    """ 
    Compute numpy.dot(scipy.linalg.toeplitz(c,r), A) using the FFT from the mkl_fft package.

    Parameters
    ----------
    c: ndarray
        the first column of the Toeplitz matrix
    r: ndarray
        the first row of the Toeplitz matrix
    A: ndarray
        the matrix to multiply on the right
    A_padded: bool, optional
        the A matrix can be pre-padded with zeros by the user, if this is the case
        set to True
    out: ndarray, optional
        an ndarray to store the output of the multiplication
    fft_len: int, optional
        specify the length of the FFT to use
    """

    if not has_mklfft:
        raise ValueError('Import mkl_fft package unavailable. Install from https://github.com/LCAV/mkl_fft')

    m = c.shape[0]
    n = r.shape[0]

    if fft_len is None:
        fft_len = int(2**np.ceil(np.log2(m+n-1)))
    zp = fft_len - m - n + 1
    
    if (not A_padded and A.shape[0] != n) or (A_padded and A.shape[0] != fft_len):
        raise ValueError('A dimensions not compatible with toeplitz(c,r)')
    
    x = np.concatenate((c, np.zeros(zp, dtype=c.dtype), r[-1:0:-1]))
    xf = fft.rfft(x, n=fft_len)
    
    if out is not None:
        fft.rfft(A, n=fft_len, axis=0, out=out)
    else:
        out = fft.rfft(A, n=fft_len, axis=0)

    out *= xf[:,None]

    if A_padded:
        fft.irfft(out, n=fft_len, axis=0, out=A)
    else:
        A = fft.irfft(out, n=fft_len, axis=0)
    
    return A[:m,] 

# 需要导入模块: from numpy import fft [as 别名]
# 或者: from numpy.fft import rfft [as 别名]
def roll_invariant_euclidean_distances(X, Y=None, squared=False):
    """
    Considering the rows of X (and Y=X) as vectors, compute the
    distance matrix between each pair of vectors.
    The distance is the minimum of the euclidean distance over all rolls:

        dist(x, y) = min_\tau(||x(t) - y(t - \tau)||^2)

    Parameters
    ----------
    X : array, shape (n_samples_1, n_features)

    Y : array, shape (n_samples_2, n_features)

    squared : boolean
        Not used. Only for API compatibility.

    Returns
    -------
    distances : array, shape (n_samples_1, n_samples_2)

    """
    X = np.atleast_2d(X)
    if Y is not None:
        Y = np.atleast_2d(Y)
    X, Y = check_pairwise_arrays(X, Y)
    n_samples_1, n_features = X.shape
    n_samples_2, n_features = Y.shape

    X_norm = np.power(np.linalg.norm(X, axis=1), 2)
    Y_norm = np.power(np.linalg.norm(Y, axis=1), 2)

    # n_pads = 0
    # n_fft = next_fast_len(n_features + n_pads)
    n_fft = n_features  # not fast but otherwise the distance is wrong
    X_hat = rfft(X, n_fft, axis=1)
    Y_hat = rfft(Y, n_fft, axis=1).conj()

    # # broadcasting can have a huge memory cost
    # XY_hat = X_hat[:, None, :] * Y_hat[None, :, :]
    # XY = irfft(XY_hat, n_fft, axis=2).max(axis=2)
    # distances = X_norm[:, None] + Y_norm[None, :] - 2 * XY

    distances = np.zeros((n_samples_1, n_samples_2))
    if n_samples_2 > 1:
        print('RIED on %s samples, this might be slow' % (distances.shape, ))
    for ii in range(n_samples_1):
        for jj in range(n_samples_2):
            XY = irfft(X_hat[ii] * Y_hat[jj], n_fft).max()
            distances[ii, jj] = X_norm[ii] + Y_norm[jj] - 2 * XY

    distances += 1e-12

    return distances