Python numpy.fft.rfft() Examples

The following are 23 code examples for showing how to use numpy.fft.rfft(). 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.fft , or try the search function .

Example 1
Project: liblsl-Python   Author: labstreaminglayer   File: PerformanceTest.py    License: MIT License 6 votes vote down vote up
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) 
Example 2
Project: pyMatrixProfile   Author: ZiyaoWei   File: util.py    License: MIT License 5 votes vote down vote up
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))) 
Example 3
Project: Computable   Author: ktraunmueller   File: bench_basic.py    License: MIT License 5 votes vote down vote up
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() 
Example 4
Project: Computable   Author: ktraunmueller   File: test_ndimage.py    License: MIT License 5 votes vote down vote up
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) 
Example 5
Project: Computable   Author: ktraunmueller   File: test_ndimage.py    License: MIT License 5 votes vote down vote up
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) 
Example 6
Project: Computable   Author: ktraunmueller   File: test_ndimage.py    License: MIT License 5 votes vote down vote up
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)) 
Example 7
Project: Computable   Author: ktraunmueller   File: test_ndimage.py    License: MIT License 5 votes vote down vote up
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) 
Example 8
Project: GraphicDesignPatternByPython   Author: Relph1119   File: test_ndimage.py    License: MIT License 5 votes vote down vote up
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) 
Example 9
Project: GraphicDesignPatternByPython   Author: Relph1119   File: test_ndimage.py    License: MIT License 5 votes vote down vote up
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) 
Example 10
Project: GraphicDesignPatternByPython   Author: Relph1119   File: test_ndimage.py    License: MIT License 5 votes vote down vote up
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) 
Example 11
Project: GraphicDesignPatternByPython   Author: Relph1119   File: test_ndimage.py    License: MIT License 5 votes vote down vote up
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) 
Example 12
Project: ComputationalNeurodynamics   Author: pmediano   File: PowerSpectrum.py    License: GNU General Public License v3.0 5 votes vote down vote up
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) 
Example 13
Project: pyvocoder   Author: shamidreza   File: lpc.py    License: GNU General Public License v2.0 5 votes vote down vote up
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) 
Example 14
Project: mnnpy   Author: chriscainx   File: irlb.py    License: BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
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:] 
Example 15
Project: pyroomacoustics   Author: LCAV   File: util.py    License: MIT License 5 votes vote down vote up
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) 
Example 16
Project: pyroomacoustics   Author: LCAV   File: util.py    License: MIT License 5 votes vote down vote up
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,] 
Example 17
Project: pyroomacoustics   Author: LCAV   File: dft.py    License: MIT License 5 votes vote down vote up
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 
Example 18
Project: pyroomacoustics   Author: LCAV   File: deconvolution.py    License: MIT License 5 votes vote down vote up
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] 
Example 19
Project: pyroomacoustics   Author: LCAV   File: sync.py    License: MIT License 5 votes vote down vote up
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)]]) 
Example 20
Project: matrixprofile-ts   Author: target   File: utils.py    License: Apache License 2.0 5 votes vote down vote up
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 :] 
Example 21
Project: gcc-nmf   Author: seanwood   File: gccNMFProcessor.py    License: MIT License 5 votes vote down vote up
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 
Example 22
Project: pyroomacoustics   Author: LCAV   File: util.py    License: MIT License 4 votes vote down vote up
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,] 
Example 23
Project: alphacsc   Author: alphacsc   File: custom_distances.py    License: BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up
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