Python numpy.fft.irfft() Examples

The following are 25 code examples for showing how to use numpy.fft.irfft(). 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: ocelot   Author: ocelot-collab   File: reswake.py    License: GNU General Public License v3.0 5 votes vote down vote up
def impedance2wake(f, y):
    """
    Fourier transform with exp(-iwt)
                f -    Hz
                y -    Om
                s    - Meter
                w -    V/C
    """
    df = f[1] - f[0]
    n = len(f)
    s = 1./df*np.arange(n)/n*speed_of_light
    w = n*df*irfft(y, n)
    return s, w 
Example 4
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 irfft as numpy_irfft

        print()
        print('Inverse 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)
            x1 = zeros(size/2+1,dtype=cdouble)
            x1[0] = x[0]
            for i in range(1,size/2):
                x1[i] = x[2*i-1] + 1j * x[2*i]
            if not size % 2:
                x1[-1] = x[-1]
            y = irfft(x)

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

            assert_array_almost_equal(numpy_irfft(x1,size),y)
            print('|%8.2f' % measure('numpy_irfft(x1,size)',repeat), end=' ')
            sys.stdout.flush()

            print(' (secs for %s calls)' % (repeat))

        sys.stdout.flush() 
Example 5
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 6
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 7
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 8
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 9
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 10
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 11
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 12
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 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 synthesis(self, X=None):
        """
        Perform time synthesis of frequency domain to real signal using the 
        inverse DFT.

        Parameters
        ----------
        X : numpy array, optional
            Complex signal in frequency domain. Default is to use DFT computed
            from ``analysis``.

        Returns
        -------
        numpy array
            IDFT of ``self.X`` or input if given.
        """

        # check for valid input
        if X is not None:
            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]))

            self.X[:,] = X

        # inverse DFT
        if self.transform == 'fftw':
            self.b[:] = self.X
            self.x[:,] = self._backward()
        elif self.transform == 'mkl':
            self.x[:,] = mkl_fft.irfft_numpy(self.X, self.nfft, axis=self.axis)
        else:
            self.x[:,] = irfft(self.X, self.nfft, axis=self.axis)

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

        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: liblsl-Python   Author: labstreaminglayer   File: PerformanceTest.py    License: MIT License 5 votes vote down vote up
def generate(self):
        X = np.random.randn(self.N // 2 + 1 + self.uneven) + 1j * np.random.randn(self.N // 2 + 1 + self.uneven)
        y = (irfft(X / self.S)).real
        if self.uneven:
            y = y[:-1]
        return normalize(y) 
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: nara_wpe   Author: fgnt   File: utils.py    License: MIT License 4 votes vote down vote up
def istft_single_channel(stft_signal, size=1024, shift=256,
          window=signal.blackman, fading=True, window_length=None,
          use_amplitude_for_biorthogonal_window=False,
          disable_sythesis_window=False):
    """
    Calculated the inverse short time Fourier transform to exactly reconstruct
    the time signal.

    Notes:
        Be careful if you make modifications in the frequency domain (e.g.
        beamforming) because the synthesis window is calculated according to the
        unmodified! analysis window.

    Args:
        stft_signal: Single channel complex STFT signal
            with dimensions frames times size/2+1.
        size: Scalar FFT-size.
        shift: Scalar FFT-shift. Typically shift is a fraction of size.
        window: Window function handle.
        fading: Removes the additional padding, if done during STFT.
        window_length: Sometimes one desires to use a shorter window than
            the fft size. In that case, the window is padded with zeros.
            The default is to use the fft-size as a window size.

    Returns:
        Single channel complex STFT signal
        Single channel time signal.
    """
    assert stft_signal.shape[1] == size // 2 + 1, str(stft_signal.shape)

    if window_length is None:
        window = window(size)
    else:
        window = window(window_length)
        window = np.pad(window, (0, size-window_length), mode='constant')
    window = _biorthogonal_window_fastest(window, shift,
                                          use_amplitude_for_biorthogonal_window)
    if disable_sythesis_window:
        window = np.ones_like(window)

    time_signal = scipy.zeros(stft_signal.shape[0] * shift + size - shift)

    for j, i in enumerate(range(0, len(time_signal) - size + shift, shift)):
        time_signal[i:i + size] += window * np.real(irfft(stft_signal[j]))

    # Compensate fade-in and fade-out
    if fading:
        time_signal = time_signal[size-shift:len(time_signal)-(size-shift)]

    return time_signal 
Example 24
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 
Example 25
Project: nolitsa   Author: manu-mannattil   File: data.py    License: BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up
def falpha(length=8192, alpha=1.0, fl=None, fu=None, mean=0.0, var=1.0):
    """Generate (1/f)^alpha noise by inverting the power spectrum.

    Generates (1/f)^alpha noise by inverting the power spectrum.
    Follows the algorithm described by Voss (1988) to generate
    fractional Brownian motion.

    Parameters
    ----------
    length : int, optional (default = 8192)
        Length of the time series to be generated.
    alpha : float, optional (default = 1.0)
        Exponent in (1/f)^alpha.  Pink noise will be generated by
        default.
    fl : float, optional (default = None)
        Lower cutoff frequency.
    fu : float, optional (default = None)
        Upper cutoff frequency.
    mean : float, optional (default = 0.0)
        Mean of the generated noise.
    var : float, optional (default = 1.0)
        Variance of the generated noise.

    Returns
    -------
    x : array
        Array containing the time series.

    Notes
    -----
    As discrete Fourier transforms assume that the input data is
    periodic, the resultant series x_{i} (= x_{i + N}) is also periodic.
    To avoid this periodicity, it is recommended to always generate
    a longer series (two or three times longer) and trim it to the
    desired length.
    """
    freqs = fft.rfftfreq(length)
    power = freqs[1:] ** -alpha
    power = np.insert(power, 0, 0)  # P(0) = 0

    if fl:
        power[freqs < fl] = 0

    if fu:
        power[freqs > fu] = 0

    # Randomize complex phases.
    phase = 2 * np.pi * np.random.random(len(freqs))
    y = np.sqrt(power) * np.exp(1j * phase)

    # The last component (corresponding to the Nyquist frequency) of an
    # RFFT with even number of points is always real.  (We don't have to
    # make the mean real as P(0) = 0.)
    if length % 2 == 0:
        y[-1] = np.abs(y[-1] * np.sqrt(2))

    x = fft.irfft(y, n=length)

    # Rescale to proper variance and mean.
    x = np.sqrt(var) * x / np.std(x)
    return mean + x - np.mean(x)