Python numpy.fft.fft() Examples

The following are code examples for showing how to use numpy.fft.fft(). They are from open source Python projects. You can vote up the examples you like or vote down the ones you don't like.

Example 1
Project: LaserTOF   Author: kyleuckert   File: fir_filter_design.py    MIT License 6 votes vote down vote up
def _dhtm(mag):
    """Compute the modified 1D discrete Hilbert transform

    Parameters
    ----------
    mag : ndarray
        The magnitude spectrum. Should be 1D with an even length, and
        preferably a fast length for FFT/IFFT.
    """
    # Adapted based on code by Niranjan Damera-Venkata,
    # Brian L. Evans and Shawn R. McCaslin (see refs for `minimum_phase`)
    sig = np.zeros(len(mag))
    # Leave Nyquist and DC at 0, knowing np.abs(fftfreq(N)[midpt]) == 0.5
    midpt = len(mag) // 2
    sig[1:midpt] = 1
    sig[midpt+1:] = -1
    # eventually if we want to support complex filters, we will need a
    # np.abs() on the mag inside the log, and should remove the .real
    recon = ifft(mag * np.exp(fft(sig * ifft(np.log(mag))))).real
    return recon 
Example 2
Project: FX-RER-Value-Extraction   Author: tsKenneth   File: fir_filter_design.py    MIT License 6 votes vote down vote up
def _dhtm(mag):
    """Compute the modified 1D discrete Hilbert transform

    Parameters
    ----------
    mag : ndarray
        The magnitude spectrum. Should be 1D with an even length, and
        preferably a fast length for FFT/IFFT.
    """
    # Adapted based on code by Niranjan Damera-Venkata,
    # Brian L. Evans and Shawn R. McCaslin (see refs for `minimum_phase`)
    sig = np.zeros(len(mag))
    # Leave Nyquist and DC at 0, knowing np.abs(fftfreq(N)[midpt]) == 0.5
    midpt = len(mag) // 2
    sig[1:midpt] = 1
    sig[midpt+1:] = -1
    # eventually if we want to support complex filters, we will need a
    # np.abs() on the mag inside the log, and should remove the .real
    recon = ifft(mag * np.exp(fft(sig * ifft(np.log(mag))))).real
    return recon 
Example 3
Project: dockerizeme   Author: dockerizeme   File: snippet.py    Apache License 2.0 6 votes vote down vote up
def fourierExtrapolation(x, n_predict):
    n = x.size
    n_harm = 10                     # number of harmonics in model
    t = np.arange(0, n)
    p = np.polyfit(t, x, 1)         # find linear trend in x
    x_notrend = x - p[0] * t        # detrended x
    x_freqdom = fft.fft(x_notrend)  # detrended x in frequency domain
    f = fft.fftfreq(n)              # frequencies
    indexes = range(n)
    # sort indexes by frequency, lower -> higher
    indexes.sort(key = lambda i: np.absolute(f[i]))
 
    t = np.arange(0, n + n_predict)
    restored_sig = np.zeros(t.size)
    for i in indexes[:1 + n_harm * 2]:
        ampli = np.absolute(x_freqdom[i]) / n   # amplitude
        phase = np.angle(x_freqdom[i])          # phase
        restored_sig += ampli * np.cos(2 * np.pi * f[i] * t + phase)
    return restored_sig + p[0] * t 
Example 4
Project: ocelot   Author: ocelot-collab   File: madx.py    GNU General Public License v3.0 6 votes vote down vote up
def get_tunes(tr_x, tr_y):
    v = np.zeros(len(tr_x))
    for i in range(len(v)): v[i] = float(tr_x[i])

    u = np.zeros(len(tr_y))
    for i in range(len(u)): u[i] = float(tr_y[i])

    sv = fft.fft(v)
    su = fft.fft(u)

    sv = np.abs(sv * np.conj(sv))
    su = np.abs(su * np.conj(su))

    freq = np.fft.fftfreq(len(sv))

    return (freq[np.argmax(sv[1:len(sv) / 2 - 1], axis=0)], freq[np.argmax(su[1:len(su) / 2 - 1], axis=0)]) 
Example 5
Project: vnpy_crypto   Author: birforce   File: fftarma.py    MIT License 6 votes vote down vote up
def fftar(self, n=None):
        '''Fourier transform of AR polynomial, zero-padded at end to n

        Parameters
        ----------
        n : int
            length of array after zero-padding

        Returns
        -------
        fftar : ndarray
            fft of zero-padded ar polynomial
        '''
        if n is None:
            n = len(self.ar)
        return fft.fft(self.padarr(self.ar, n)) 
Example 6
Project: vnpy_crypto   Author: birforce   File: fftarma.py    MIT License 6 votes vote down vote up
def fftma(self, n):
        '''Fourier transform of MA polynomial, zero-padded at end to n

        Parameters
        ----------
        n : int
            length of array after zero-padding

        Returns
        -------
        fftar : ndarray
            fft of zero-padded ar polynomial
        '''
        if n is None:
            n = len(self.ar)
        return fft.fft(self.padarr(self.ma, n))

    #@OneTimeProperty  # not while still debugging things 
Example 7
Project: vnpy_crypto   Author: birforce   File: fftarma.py    MIT License 6 votes vote down vote up
def fftarma(self, n=None):
        '''Fourier transform of ARMA polynomial, zero-padded at end to n

        The Fourier transform of the ARMA process is calculated as the ratio
        of the fft of the MA polynomial divided by the fft of the AR polynomial.

        Parameters
        ----------
        n : int
            length of array after zero-padding

        Returns
        -------
        fftarma : ndarray
            fft of zero-padded arma polynomial
        '''
        if n is None:
            n = self.nobs
        return (self.fftma(n) / self.fftar(n)) 
Example 8
Project: vnpy_crypto   Author: birforce   File: fftarma.py    MIT License 6 votes vote down vote up
def filter(self, x):
        '''
        filter a timeseries with the ARMA filter

        padding with zero is missing, in example I needed the padding to get
        initial conditions identical to direct filter

        Initial filtered observations differ from filter2 and signal.lfilter, but
        at end they are the same.

        See Also
        --------
        tsa.filters.fftconvolve

        '''
        n = x.shape[0]
        if n == self.fftarma:
            fftarma = self.fftarma
        else:
            fftarma = self.fftma(n) / self.fftar(n)
        tmpfft = fftarma * fft.fft(x)
        return fft.ifft(tmpfft) 
Example 9
Project: vnpy_crypto   Author: birforce   File: fftarma.py    MIT License 6 votes vote down vote up
def invpowerspd(self, n):
        '''autocovariance from spectral density

        scaling is correct, but n needs to be large for numerical accuracy
        maybe padding with zero in fft would be faster
        without slicing it returns 2-sided autocovariance with fftshift

        >>> ArmaFft([1, -0.5], [1., 0.4], 40).invpowerspd(2**8)[:10]
        array([ 2.08    ,  1.44    ,  0.72    ,  0.36    ,  0.18    ,  0.09    ,
                0.045   ,  0.0225  ,  0.01125 ,  0.005625])
        >>> ArmaFft([1, -0.5], [1., 0.4], 40).acovf(10)
        array([ 2.08    ,  1.44    ,  0.72    ,  0.36    ,  0.18    ,  0.09    ,
                0.045   ,  0.0225  ,  0.01125 ,  0.005625])
        '''
        hw = self.fftarma(n)
        return np.real_if_close(fft.ifft(hw*hw.conj()), tol=200)[:n] 
Example 10
Project: ble5-nrf52-mac   Author: tomasero   File: fir_filter_design.py    MIT License 6 votes vote down vote up
def _dhtm(mag):
    """Compute the modified 1D discrete Hilbert transform

    Parameters
    ----------
    mag : ndarray
        The magnitude spectrum. Should be 1D with an even length, and
        preferably a fast length for FFT/IFFT.
    """
    # Adapted based on code by Niranjan Damera-Venkata,
    # Brian L. Evans and Shawn R. McCaslin (see refs for `minimum_phase`)
    sig = np.zeros(len(mag))
    # Leave Nyquist and DC at 0, knowing np.abs(fftfreq(N)[midpt]) == 0.5
    midpt = len(mag) // 2
    sig[1:midpt] = 1
    sig[midpt+1:] = -1
    # eventually if we want to support complex filters, we will need a
    # np.abs() on the mag inside the log, and should remove the .real
    recon = ifft(mag * np.exp(fft(sig * ifft(np.log(mag))))).real
    return recon 
Example 11
Project: Test-stock-prediction-algorithms   Author: timestocome   File: FFT_Stock_Prediction.py    MIT License 6 votes vote down vote up
def fourierEx(x, n_predict, harmonics):

    n = len(x)                  # number of input samples
    n_harmonics = harmonics            # f, 2*f, 3*f, .... n_harmonics  ( 1,2, )
    t = np.arange(0, n)         # place to store data
    p = np.polyfit(t, x, 1)     # find trend
    x_no_trend = x - p[0] * t 
    x_frequency_domains = fft.fft(x_no_trend)
    f = np.fft.fftfreq(n)       # frequencies
    indexes = list(range(n))
    indexes.sort(key=lambda i: np.absolute(f[i]))

    t = np.arange(0, n + n_predict)
    restored_signal = np.zeros(t.size)
    for i in indexes[:1 + n_harmonics * 2]:
        amplitude = np.absolute(x_frequency_domains[i] / n)
        phase = np.angle(x_frequency_domains[i])
        restored_signal += amplitude * np.cos(2 * np.pi * f[i] * t + phase)
    
    return restored_signal + p[0] * t 
Example 12
Project: P3_image_processing   Author: latedude2   File: fir_filter_design.py    MIT License 6 votes vote down vote up
def _dhtm(mag):
    """Compute the modified 1D discrete Hilbert transform

    Parameters
    ----------
    mag : ndarray
        The magnitude spectrum. Should be 1D with an even length, and
        preferably a fast length for FFT/IFFT.
    """
    # Adapted based on code by Niranjan Damera-Venkata,
    # Brian L. Evans and Shawn R. McCaslin (see refs for `minimum_phase`)
    sig = np.zeros(len(mag))
    # Leave Nyquist and DC at 0, knowing np.abs(fftfreq(N)[midpt]) == 0.5
    midpt = len(mag) // 2
    sig[1:midpt] = 1
    sig[midpt+1:] = -1
    # eventually if we want to support complex filters, we will need a
    # np.abs() on the mag inside the log, and should remove the .real
    recon = ifft(mag * np.exp(fft(sig * ifft(np.log(mag))))).real
    return recon 
Example 13
Project: LaserTOF   Author: kyleuckert   File: test_basic.py    MIT License 5 votes vote down vote up
def test_djbfft(self):
        from numpy.fft import fft as numpy_fft
        for i in range(2,14):
            n = 2**i
            x = list(range(n))
            y2 = numpy_fft(x)
            y1 = zeros((n,),dtype=double)
            y1[0] = y2[0].real
            y1[-1] = y2[n//2].real
            for k in range(1, n//2):
                y1[2*k-1] = y2[k].real
                y1[2*k] = y2[k].imag
            y = fftpack.drfft(x)
            assert_array_almost_equal(y,y1) 
Example 14
Project: LaserTOF   Author: kyleuckert   File: test_max_len_seq.py    MIT License 5 votes vote down vote up
def test_mls_output(self):
        # define some alternate working taps
        alt_taps = {2: [1], 3: [2], 4: [3], 5: [4, 3, 2], 6: [5, 4, 1], 7: [4],
                    8: [7, 5, 3]}
        # assume the other bit levels work, too slow to test higher orders...
        for nbits in range(2, 8):
            for state in [None, np.round(np.random.rand(nbits))]:
                for taps in [None, alt_taps[nbits]]:
                    if state is not None and np.all(state == 0):
                        state[0] = 1  # they can't all be zero
                    orig_m = max_len_seq(nbits, state=state,
                                         taps=taps)[0]
                    m = 2. * orig_m - 1.  # convert to +/- 1 representation
                    # First, make sure we got all 1's or -1
                    err_msg = "mls had non binary terms"
                    assert_array_equal(np.abs(m), np.ones_like(m),
                                       err_msg=err_msg)
                    # Test via circular cross-correlation, which is just mult.
                    # in the frequency domain with one signal conjugated
                    tester = np.real(ifft(fft(m) * np.conj(fft(m))))
                    out_len = 2**nbits - 1
                    # impulse amplitude == test_len
                    err_msg = "mls impulse has incorrect value"
                    assert_allclose(tester[0], out_len, err_msg=err_msg)
                    # steady-state is -1
                    err_msg = "mls steady-state has incorrect value"
                    assert_allclose(tester[1:], -1 * np.ones(out_len - 1),
                                    err_msg=err_msg)
                    # let's do the split thing using a couple options
                    for n in (1, 2**(nbits - 1)):
                        m1, s1 = max_len_seq(nbits, state=state, taps=taps,
                                             length=n)
                        m2, s2 = max_len_seq(nbits, state=s1, taps=taps,
                                             length=1)
                        m3, s3 = max_len_seq(nbits, state=s2, taps=taps,
                                             length=out_len - n - 1)
                        new_m = np.concatenate((m1, m2, m3))
                        assert_array_equal(orig_m, new_m) 
Example 15
Project: LaserTOF   Author: kyleuckert   File: test_ndimage.py    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 16
Project: LaserTOF   Author: kyleuckert   File: test_ndimage.py    MIT License 5 votes vote down vote up
def test_fourier_gaussian_complex01(self):
        for shape in [(32, 16), (31, 15)]:
            for type in [numpy.complex64, numpy.complex128]:
                a = numpy.zeros(shape, type)
                a[0, 0] = 1.0
                a = fft.fft(a, shape[0], 0)
                a = fft.fft(a, shape[1], 1)
                a = ndimage.fourier_gaussian(a, [5.0, 2.5], -1,
                                                       0)
                a = fft.ifft(a, shape[1], 1)
                a = fft.ifft(a, shape[0], 0)
                assert_almost_equal(ndimage.sum(a.real), 1.0) 
Example 17
Project: LaserTOF   Author: kyleuckert   File: test_ndimage.py    MIT License 5 votes vote down vote up
def test_fourier_uniform_complex01(self):
        for shape in [(32, 16), (31, 15)]:
            for type in [numpy.complex64, numpy.complex128]:
                a = numpy.zeros(shape, type)
                a[0, 0] = 1.0
                a = fft.fft(a, shape[0], 0)
                a = fft.fft(a, shape[1], 1)
                a = ndimage.fourier_uniform(a, [5.0, 2.5], -1, 0)
                a = fft.ifft(a, shape[1], 1)
                a = fft.ifft(a, shape[0], 0)
                assert_almost_equal(ndimage.sum(a.real), 1.0) 
Example 18
Project: LaserTOF   Author: kyleuckert   File: test_ndimage.py    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 19
Project: LaserTOF   Author: kyleuckert   File: test_ndimage.py    MIT License 5 votes vote down vote up
def test_fourier_shift_complex01(self):
        for shape in [(32, 16), (31, 15)]:
            for type in [numpy.complex64, numpy.complex128]:
                expected = numpy.arange(shape[0] * shape[1],
                                       dtype=type)
                expected.shape = shape
                a = fft.fft(expected, shape[0], 0)
                a = fft.fft(a, shape[1], 1)
                a = ndimage.fourier_shift(a, [1, 1], -1, 0)
                a = fft.ifft(a, shape[1], 1)
                a = fft.ifft(a, shape[0], 0)
                assert_array_almost_equal(a.real[1:, 1:], expected[:-1, :-1])
                assert_array_almost_equal(a.imag, numpy.zeros(shape)) 
Example 20
Project: LaserTOF   Author: kyleuckert   File: test_ndimage.py    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 21
Project: LaserTOF   Author: kyleuckert   File: test_ndimage.py    MIT License 5 votes vote down vote up
def test_fourier_ellipsoid_complex01(self):
        for shape in [(32, 16), (31, 15)]:
            for type in [numpy.complex64, numpy.complex128]:
                a = numpy.zeros(shape, type)
                a[0, 0] = 1.0
                a = fft.fft(a, shape[0], 0)
                a = fft.fft(a, shape[1], 1)
                a = ndimage.fourier_ellipsoid(a, [5.0, 2.5], -1,
                                                        0)
                a = fft.ifft(a, shape[1], 1)
                a = fft.ifft(a, shape[0], 0)
                assert_almost_equal(ndimage.sum(a.real), 1.0) 
Example 22
Project: FX-RER-Value-Extraction   Author: tsKenneth   File: test_basic.py    MIT License 5 votes vote down vote up
def test_djbfft(self):
        from numpy.fft import fft as numpy_fft
        for i in range(2,14):
            n = 2**i
            x = list(range(n))
            y2 = numpy_fft(x)
            y1 = zeros((n,),dtype=double)
            y1[0] = y2[0].real
            y1[-1] = y2[n//2].real
            for k in range(1, n//2):
                y1[2*k-1] = y2[k].real
                y1[2*k] = y2[k].imag
            y = fftpack.drfft(x)
            assert_array_almost_equal(y,y1) 
Example 23
Project: FX-RER-Value-Extraction   Author: tsKenneth   File: test_max_len_seq.py    MIT License 5 votes vote down vote up
def test_mls_output(self):
        # define some alternate working taps
        alt_taps = {2: [1], 3: [2], 4: [3], 5: [4, 3, 2], 6: [5, 4, 1], 7: [4],
                    8: [7, 5, 3]}
        # assume the other bit levels work, too slow to test higher orders...
        for nbits in range(2, 8):
            for state in [None, np.round(np.random.rand(nbits))]:
                for taps in [None, alt_taps[nbits]]:
                    if state is not None and np.all(state == 0):
                        state[0] = 1  # they can't all be zero
                    orig_m = max_len_seq(nbits, state=state,
                                         taps=taps)[0]
                    m = 2. * orig_m - 1.  # convert to +/- 1 representation
                    # First, make sure we got all 1's or -1
                    err_msg = "mls had non binary terms"
                    assert_array_equal(np.abs(m), np.ones_like(m),
                                       err_msg=err_msg)
                    # Test via circular cross-correlation, which is just mult.
                    # in the frequency domain with one signal conjugated
                    tester = np.real(ifft(fft(m) * np.conj(fft(m))))
                    out_len = 2**nbits - 1
                    # impulse amplitude == test_len
                    err_msg = "mls impulse has incorrect value"
                    assert_allclose(tester[0], out_len, err_msg=err_msg)
                    # steady-state is -1
                    err_msg = "mls steady-state has incorrect value"
                    assert_allclose(tester[1:], -1 * np.ones(out_len - 1),
                                    err_msg=err_msg)
                    # let's do the split thing using a couple options
                    for n in (1, 2**(nbits - 1)):
                        m1, s1 = max_len_seq(nbits, state=state, taps=taps,
                                             length=n)
                        m2, s2 = max_len_seq(nbits, state=s1, taps=taps,
                                             length=1)
                        m3, s3 = max_len_seq(nbits, state=s2, taps=taps,
                                             length=out_len - n - 1)
                        new_m = np.concatenate((m1, m2, m3))
                        assert_array_equal(orig_m, new_m) 
Example 24
Project: att   Author: Centre-Alt-Rendiment-Esportiu   File: bench_basic.py    GNU General Public License v3.0 5 votes vote down vote up
def bench_random(self):
        from numpy.fft import fft as numpy_fft
        print()
        print('                 Fast Fourier Transform')
        print('=================================================')
        print('      |    real input     |   complex input    ')
        print('-------------------------------------------------')
        print(' size |  scipy  |  numpy  |  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()

            for x in [random([size]).astype(double),
                      random([size]).astype(cdouble)+random([size]).astype(cdouble)*1j
                      ]:
                if size > 500:
                    y = fft(x)
                else:
                    y = direct_dft(x)
                assert_array_almost_equal(fft(x),y)
                print('|%8.2f' % measure('fft(x)',repeat), end=' ')
                sys.stdout.flush()

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

            print(' (secs for %s calls)' % (repeat))
        sys.stdout.flush() 
Example 25
Project: att   Author: Centre-Alt-Rendiment-Esportiu   File: test_basic.py    GNU General Public License v3.0 5 votes vote down vote up
def test_djbfft(self):
        from numpy.fft import fft as numpy_fft
        for i in range(2,14):
            n = 2**i
            x = list(range(n))
            y2 = numpy_fft(x)
            y1 = zeros((n,),dtype=double)
            y1[0] = y2[0].real
            y1[-1] = y2[n//2].real
            for k in range(1, n//2):
                y1[2*k-1] = y2[k].real
                y1[2*k] = y2[k].imag
            y = fftpack.drfft(x)
            assert_array_almost_equal(y,y1) 
Example 26
Project: att   Author: Centre-Alt-Rendiment-Esportiu   File: test_max_len_seq.py    GNU General Public License v3.0 5 votes vote down vote up
def test_mls_output(self):
        # define some alternate working taps
        alt_taps = {2: [1], 3: [2], 4: [3], 5: [4, 3, 2], 6: [5, 4, 1], 7: [4],
                    8: [7, 5, 3]}
        # assume the other bit levels work, too slow to test higher orders...
        for nbits in range(2, 8):
            for state in [None, np.round(np.random.rand(nbits))]:
                for taps in [None, alt_taps[nbits]]:
                    if state is not None and np.all(state == 0):
                        state[0] = 1  # they can't all be zero
                    orig_m = max_len_seq(nbits, state=state,
                                         taps=taps)[0]
                    m = 2. * orig_m - 1.  # convert to +/- 1 representation
                    # First, make sure we got all 1's or -1
                    err_msg = "mls had non binary terms"
                    assert_array_equal(np.abs(m), np.ones_like(m),
                                       err_msg=err_msg)
                    # Test via circular cross-correlation, which is just mult.
                    # in the frequency domain with one signal conjugated
                    tester = np.real(ifft(fft(m) * np.conj(fft(m))))
                    out_len = 2**nbits - 1
                    # impulse amplitude == test_len
                    err_msg = "mls impulse has incorrect value"
                    assert_allclose(tester[0], out_len, err_msg=err_msg)
                    # steady-state is -1
                    err_msg = "mls steady-state has incorrect value"
                    assert_allclose(tester[1:], -1 * np.ones(out_len - 1),
                                    err_msg=err_msg)
                    # let's do the split thing using a couple options
                    for n in (1, 2**(nbits - 1)):
                        m1, s1 = max_len_seq(nbits, state=state, taps=taps,
                                             length=n)
                        m2, s2 = max_len_seq(nbits, state=s1, taps=taps,
                                             length=1)
                        m3, s3 = max_len_seq(nbits, state=s2, taps=taps,
                                             length=out_len - n - 1)
                        new_m = np.concatenate((m1, m2, m3))
                        assert_array_equal(orig_m, new_m) 
Example 27
Project: att   Author: Centre-Alt-Rendiment-Esportiu   File: test_ndimage.py    GNU General Public License v3.0 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 28
Project: att   Author: Centre-Alt-Rendiment-Esportiu   File: test_ndimage.py    GNU General Public License v3.0 5 votes vote down vote up
def test_fourier_gaussian_complex01(self):
        for shape in [(32, 16), (31, 15)]:
            for type in [numpy.complex64, numpy.complex128]:
                a = numpy.zeros(shape, type)
                a[0, 0] = 1.0
                a = fft.fft(a, shape[0], 0)
                a = fft.fft(a, shape[1], 1)
                a = ndimage.fourier_gaussian(a, [5.0, 2.5], -1,
                                                       0)
                a = fft.ifft(a, shape[1], 1)
                a = fft.ifft(a, shape[0], 0)
                assert_almost_equal(ndimage.sum(a.real), 1.0) 
Example 29
Project: att   Author: Centre-Alt-Rendiment-Esportiu   File: test_ndimage.py    GNU General Public License v3.0 5 votes vote down vote up
def test_fourier_uniform_complex01(self):
        for shape in [(32, 16), (31, 15)]:
            for type in [numpy.complex64, numpy.complex128]:
                a = numpy.zeros(shape, type)
                a[0, 0] = 1.0
                a = fft.fft(a, shape[0], 0)
                a = fft.fft(a, shape[1], 1)
                a = ndimage.fourier_uniform(a, [5.0, 2.5], -1, 0)
                a = fft.ifft(a, shape[1], 1)
                a = fft.ifft(a, shape[0], 0)
                assert_almost_equal(ndimage.sum(a.real), 1.0) 
Example 30
Project: att   Author: Centre-Alt-Rendiment-Esportiu   File: test_ndimage.py    GNU General Public License v3.0 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 31
Project: att   Author: Centre-Alt-Rendiment-Esportiu   File: test_ndimage.py    GNU General Public License v3.0 5 votes vote down vote up
def test_fourier_shift_complex01(self):
        for shape in [(32, 16), (31, 15)]:
            for type in [numpy.complex64, numpy.complex128]:
                expected = numpy.arange(shape[0] * shape[1],
                                       dtype=type)
                expected.shape = shape
                a = fft.fft(expected, shape[0], 0)
                a = fft.fft(a, shape[1], 1)
                a = ndimage.fourier_shift(a, [1, 1], -1, 0)
                a = fft.ifft(a, shape[1], 1)
                a = fft.ifft(a, shape[0], 0)
                assert_array_almost_equal(a.real[1:, 1:], expected[:-1, :-1])
                assert_array_almost_equal(a.imag, numpy.zeros(shape)) 
Example 32
Project: att   Author: Centre-Alt-Rendiment-Esportiu   File: test_ndimage.py    GNU General Public License v3.0 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 33
Project: att   Author: Centre-Alt-Rendiment-Esportiu   File: test_ndimage.py    GNU General Public License v3.0 5 votes vote down vote up
def test_fourier_ellipsoid_complex01(self):
        for shape in [(32, 16), (31, 15)]:
            for type in [numpy.complex64, numpy.complex128]:
                a = numpy.zeros(shape, type)
                a[0, 0] = 1.0
                a = fft.fft(a, shape[0], 0)
                a = fft.fft(a, shape[1], 1)
                a = ndimage.fourier_ellipsoid(a, [5.0, 2.5], -1,
                                                        0)
                a = fft.ifft(a, shape[1], 1)
                a = fft.ifft(a, shape[0], 0)
                assert_almost_equal(ndimage.sum(a.real), 1.0) 
Example 34
Project: vawt-wake-model   Author: byuflowlab   File: Poisson_Solver.py    MIT License 5 votes vote down vote up
def dst(x,axis=-1):
    """Discrete Sine Transform (DST-I)

    Implemented using 2(N+1)-point FFT
    xsym = r_[0,x,0,-x[::-1]]
    DST = (-imag(fft(xsym))/2)[1:(N+1)]

    adjusted to work over an arbitrary axis for entire n-dim array
    """
    n = len(x.shape)
    N = x.shape[axis]
    slices = [None]*3
    for k in range(3):
        slices[k] = []
        for j in range(n):
            slices[k].append(slice(None))
    newshape = list(x.shape)
    newshape[axis] = 2*(N+1)
    xsym = np.zeros(newshape,np.float)
    slices[0][axis] = slice(1,N+1)
    slices[1][axis] = slice(N+2,None)
    slices[2][axis] = slice(None,None,-1)
    for k in range(3):
        slices[k] = tuple(slices[k])
    xsym[slices[0]] = x
    xsym[slices[1]] = -x[slices[2]]
    DST = fft(xsym,axis=axis)
    #print xtilde
    return (-(DST.imag)/2)[slices[0]] 
Example 35
Project: ocelot   Author: ocelot-collab   File: csr.py    GNU General Public License v3.0 5 votes vote down vote up
def csr_convolution(a, b):
    P = len(a)
    Q = len(b)
    L = P + Q - 1
    K = 2 ** nextpow2(L)
    a_pad = np.pad(a, (0, K - P), 'constant', constant_values=(0))
    b_pad = np.pad(b, (0, K - Q), 'constant', constant_values=(0))
    c = ifft(fft(a_pad)*fft(b_pad))
    c = c[0:L-1].real
    return c 
Example 36
Project: ocelot   Author: ocelot-collab   File: reswake.py    GNU General Public License v3.0 5 votes vote down vote up
def wake2impedance(s, w):
    """
    Fourier transform with exp(iwt)
                 s    - Meter
                 w -    V/C
                 f -    Hz
                 y -    Om
    """
    ds = s[1] - s[0]
    dt = ds/speed_of_light
    n = len(s)
    shift = 1.
    y = dt*fft(w, n)*shift
    return y 
Example 37
Project: vnpy_crypto   Author: birforce   File: fftarma.py    MIT License 5 votes vote down vote up
def spd(self, npos):
        '''raw spectral density, returns Fourier transform

        n is number of points in positive spectrum, the actual number of points
        is twice as large. different from other spd methods with fft
        '''
        n = npos
        w = fft.fftfreq(2*n) * 2 * np.pi
        hw = self.fftarma(2*n)  #not sure, need to check normalization
        #return (hw*hw.conj()).real[n//2-1:]  * 0.5 / np.pi #doesn't show in plot
        return (hw*hw.conj()).real * 0.5 / np.pi, w 
Example 38
Project: vnpy_crypto   Author: birforce   File: fftarma.py    MIT License 5 votes vote down vote up
def spdshift(self, n):
        '''power spectral density using fftshift

        currently returns two-sided according to fft frequencies, use first half
        '''
        #size = s1+s2-1
        mapadded = self.padarr(self.ma, n)
        arpadded = self.padarr(self.ar, n)
        hw = fft.fft(fft.fftshift(mapadded)) / fft.fft(fft.fftshift(arpadded))
        #return np.abs(spd)[n//2-1:]
        w = fft.fftfreq(n) * 2 * np.pi
        wslice = slice(n//2-1, None, None)
        #return (hw*hw.conj()).real[wslice], w[wslice]
        return (hw*hw.conj()).real, w 
Example 39
Project: vnpy_crypto   Author: birforce   File: fftarma.py    MIT License 5 votes vote down vote up
def _spddirect2(self, n):
        '''this looks bad, maybe with an fftshift
        '''
        #size = s1+s2-1
        hw = (fft.fft(np.r_[self.ma[::-1],self.ma], n)
                / fft.fft(np.r_[self.ar[::-1],self.ar], n))
        return (hw*hw.conj()) #.real[n//2-1:] 
Example 40
Project: SEELablet   Author: jithinbp   File: analyticsClass.py    GNU General Public License v3.0 5 votes vote down vote up
def find_frequency(self, v, si):  # voltages, samplimg interval is seconds
		from numpy import fft
		NP = len(v)
		v = v -v.mean()			# remove DC component
		frq = fft.fftfreq(NP, si)[:NP/2]	# take only the +ive half of the frequncy array
		amp = abs(fft.fft(v)[:NP/2])/NP		# and the fft result
		index =  amp.argmax()				# search for the tallest peak, the fundamental
		return frq[index] 
Example 41
Project: SEELablet   Author: jithinbp   File: analyticsClass.py    GNU General Public License v3.0 5 votes vote down vote up
def amp_spectrum(self, v, si, nhar=8):  
		# voltages, samplimg interval is seconds, number of harmonics to retain
		from numpy import fft
		NP = len(v)
		frq = fft.fftfreq(NP, si)[:NP/2]	# take only the +ive half of the frequncy array
		amp = abs(fft.fft(v)[:NP/2])/NP		# and the fft result
		index =  amp.argmax()				# search for the tallest peak, the fundamental
		if index == 0:						# DC component is dominating
			index =  amp[4:].argmax()		# skip frequencies close to zero
		return frq[:index*nhar], amp[:index*nhar]	# restrict to 'nhar' harmonics 
Example 42
Project: SEELablet   Author: jithinbp   File: analyticsClass.py    GNU General Public License v3.0 5 votes vote down vote up
def fft(self,ya, si):
		'''
		Returns positive half of the Fourier transform of the signal ya. 
		Sampling interval 'si', in milliseconds
		'''
		ns = len(ya)
		if ns %2 == 1:  # odd values of np give exceptions
			ns=ns-1 # make it even
			ya=ya[:-1]
		v = np.array(ya)
		tr = abs(np.fft.fft(v))/ns
		frq = np.fft.fftfreq(ns, si)
		x = frq.reshape(2,ns/2)
		y = tr.reshape(2,ns/2)
		return x[0], y[0] 
Example 43
Project: ble5-nrf52-mac   Author: tomasero   File: test_basic.py    MIT License 5 votes vote down vote up
def test_djbfft(self):
        from numpy.fft import fft as numpy_fft
        for i in range(2,14):
            n = 2**i
            x = list(range(n))
            y2 = numpy_fft(x)
            y1 = zeros((n,),dtype=double)
            y1[0] = y2[0].real
            y1[-1] = y2[n//2].real
            for k in range(1, n//2):
                y1[2*k-1] = y2[k].real
                y1[2*k] = y2[k].imag
            y = fftpack.drfft(x)
            assert_array_almost_equal(y,y1) 
Example 44
Project: ble5-nrf52-mac   Author: tomasero   File: test_max_len_seq.py    MIT License 5 votes vote down vote up
def test_mls_output(self):
        # define some alternate working taps
        alt_taps = {2: [1], 3: [2], 4: [3], 5: [4, 3, 2], 6: [5, 4, 1], 7: [4],
                    8: [7, 5, 3]}
        # assume the other bit levels work, too slow to test higher orders...
        for nbits in range(2, 8):
            for state in [None, np.round(np.random.rand(nbits))]:
                for taps in [None, alt_taps[nbits]]:
                    if state is not None and np.all(state == 0):
                        state[0] = 1  # they can't all be zero
                    orig_m = max_len_seq(nbits, state=state,
                                         taps=taps)[0]
                    m = 2. * orig_m - 1.  # convert to +/- 1 representation
                    # First, make sure we got all 1's or -1
                    err_msg = "mls had non binary terms"
                    assert_array_equal(np.abs(m), np.ones_like(m),
                                       err_msg=err_msg)
                    # Test via circular cross-correlation, which is just mult.
                    # in the frequency domain with one signal conjugated
                    tester = np.real(ifft(fft(m) * np.conj(fft(m))))
                    out_len = 2**nbits - 1
                    # impulse amplitude == test_len
                    err_msg = "mls impulse has incorrect value"
                    assert_allclose(tester[0], out_len, err_msg=err_msg)
                    # steady-state is -1
                    err_msg = "mls steady-state has incorrect value"
                    assert_allclose(tester[1:], -1 * np.ones(out_len - 1),
                                    err_msg=err_msg)
                    # let's do the split thing using a couple options
                    for n in (1, 2**(nbits - 1)):
                        m1, s1 = max_len_seq(nbits, state=state, taps=taps,
                                             length=n)
                        m2, s2 = max_len_seq(nbits, state=s1, taps=taps,
                                             length=1)
                        m3, s3 = max_len_seq(nbits, state=s2, taps=taps,
                                             length=out_len - n - 1)
                        new_m = np.concatenate((m1, m2, m3))
                        assert_array_equal(orig_m, new_m) 
Example 45
Project: ble5-nrf52-mac   Author: tomasero   File: test_ndimage.py    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 46
Project: ble5-nrf52-mac   Author: tomasero   File: test_ndimage.py    MIT License 5 votes vote down vote up
def test_fourier_gaussian_complex01(self):
        for shape in [(32, 16), (31, 15)]:
            for type_, dec in zip([numpy.complex64, numpy.complex128], [6, 14]):
                a = numpy.zeros(shape, type_)
                a[0, 0] = 1.0
                a = fft.fft(a, shape[0], 0)
                a = fft.fft(a, shape[1], 1)
                a = ndimage.fourier_gaussian(a, [5.0, 2.5], -1, 0)
                a = fft.ifft(a, shape[1], 1)
                a = fft.ifft(a, shape[0], 0)
                assert_almost_equal(ndimage.sum(a.real), 1.0, decimal=dec) 
Example 47
Project: ble5-nrf52-mac   Author: tomasero   File: test_ndimage.py    MIT License 5 votes vote down vote up
def test_fourier_uniform_complex01(self):
        for shape in [(32, 16), (31, 15)]:
            for type_, dec in zip([numpy.complex64, numpy.complex128], [6, 14]):
                a = numpy.zeros(shape, type_)
                a[0, 0] = 1.0
                a = fft.fft(a, shape[0], 0)
                a = fft.fft(a, shape[1], 1)
                a = ndimage.fourier_uniform(a, [5.0, 2.5], -1, 0)
                a = fft.ifft(a, shape[1], 1)
                a = fft.ifft(a, shape[0], 0)
                assert_almost_equal(ndimage.sum(a.real), 1.0, decimal=dec) 
Example 48
Project: ble5-nrf52-mac   Author: tomasero   File: test_ndimage.py    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 49
Project: ble5-nrf52-mac   Author: tomasero   File: test_ndimage.py    MIT License 5 votes vote down vote up
def test_fourier_shift_complex01(self):
        for shape in [(32, 16), (31, 15)]:
            for type_, dec in zip([numpy.complex64, numpy.complex128], [4, 11]):
                expected = numpy.arange(shape[0] * shape[1], dtype=type_)
                expected.shape = shape
                a = fft.fft(expected, shape[0], 0)
                a = fft.fft(a, shape[1], 1)
                a = ndimage.fourier_shift(a, [1, 1], -1, 0)
                a = fft.ifft(a, shape[1], 1)
                a = fft.ifft(a, shape[0], 0)
                assert_array_almost_equal(a.real[1:, 1:], expected[:-1, :-1],
                                          decimal=dec)
                assert_array_almost_equal(a.imag, numpy.zeros(shape),
                                          decimal=dec) 
Example 50
Project: ble5-nrf52-mac   Author: tomasero   File: test_ndimage.py    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)