Python numpy.fft.fft() Examples

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 

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 

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 

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)]) 

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)) 

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 

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)) 

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) 

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] 

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 

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 

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 

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) 

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) 

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) 

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) 

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) 

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)) 

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)) 

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) 

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) 

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) 

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) 

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() 

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) 

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) 

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) 

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) 

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) 

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)) 

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)) 

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) 

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) 

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]] 

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 

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 

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 

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 

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:] 

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] 

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 

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] 

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) 

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) 

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) 

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) 

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) 

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) 

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) 

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)