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 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 mode(self, on):
        if on:
            self.changeState=1
            self.current=self.pwspec
            self.btnMode.setText("scope")
            self.btnMode.setIcon(Qt.QIcon(Qt.QPixmap(icons.scope)))
            self.btnMode.setChecked(True)
        else:
            self.changeState=0
            self.current=self.scope
            self.btnMode.setText("fft")
            self.btnMode.setIcon(Qt.QIcon(Qt.QPixmap(icons.pwspec)))
            self.btnMode.setChecked(False)
        if self.changeState==1:
            self.stack.setCurrentIndex(self.changeState)
            self.scope.plot.setDatastream(None)
            self.pwspec.plot.setDatastream(stream)
        else:
            self.stack.setCurrentIndex(self.changeState)
            self.pwspec.plot.setDatastream(None)
            self.scope.plot.setDatastream(stream) 

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 to_powerspec(x, sr, win_len, time_step):
    nperseg = int(win_len*sr)
    nperstep = int(time_step*sr)
    nfft = int(2**(ceil(log(nperseg)/log(2))))
    window = hanning(nperseg+2)[1:nperseg+1]

    indices = arange(int(nperseg/2), x.shape[0] - int(nperseg/2) + 1, nperstep)
    num_frames = len(indices)

    pspec = zeros((num_frames,int(nfft/2)+1))
    for i in range(num_frames):
        X = x[indices[i]-int(nperseg/2):indices[i]+int(nperseg/2)]
        X = X * window
        fx = fft(X, n = nfft)
        power = abs(fx[:int(nfft/2)+1])**2
        pspec[i,:] = power
    return pspec 

def cconv(a, b):
    """
    Circular convolution of vectors

    Computes the circular convolution of two vectors a and b via their
    fast fourier transforms

    a \ast b = \mathcal{F}^{-1}(\mathcal{F}(a) \odot \mathcal{F}(b))

    Parameter
    ---------
    a: real valued array (shape N)
    b: real valued array (shape N)

    Returns
    -------
    c: real valued array (shape N), representing the circular
       convolution of a and b
    """
    return ifft(fft(a) * fft(b)).real 

def ccorr(a, b):
    """
    Circular correlation of vectors

    Computes the circular correlation of two vectors a and b via their
    fast fourier transforms

    a \ast b = \mathcal{F}^{-1}(\overline{\mathcal{F}(a)} \odot \mathcal{F}(b))

    Parameter
    ---------
    a: real valued array (shape N)
    b: real valued array (shape N)

    Returns
    -------
    c: real valued array (shape N), representing the circular
       correlation of a and b
    """

    return ifft(np.conj(fft(a)) * fft(b)).real 

def cconv(a, b):
	"""
	Circular convolution of vectors

	Computes the circular convolution of two vectors a and b via their
	fast fourier transforms

	a \ast b = \mathcal{F}^{-1}(\mathcal{F}(a) \odot \mathcal{F}(b))

	Parameter
	---------
	a: real valued array (shape N)
	b: real valued array (shape N)

	Returns
	-------
	c: real valued array (shape N), representing the circular
	   convolution of a and b
	"""
	return ifft(fft(a) * fft(b)).real 

def ccorr(a, b):
	"""
	Circular correlation of vectors

	Computes the circular correlation of two vectors a and b via their
	fast fourier transforms

	a \ast b = \mathcal{F}^{-1}(\overline{\mathcal{F}(a)} \odot \mathcal{F}(b))

	Parameter
	---------
	a: real valued array (shape N)
	b: real valued array (shape N)

	Returns
	-------
	c: real valued array (shape N), representing the circular
	   correlation of a and b
	"""

	return ifft(np.conj(fft(a)) * fft(b)).real 

def _frft2(x, alpha):
    assert x.ndim == 2, "x must be a 2 dimensional array"
    m, n = x.shape
    # TODO please remove this confusing comment. Is it 'm' or 'm-1' ?
    # TODO If 'p = m', more code cleaning is easy to do.
    p = m  # m-1 # deveria incrementarse el sigiente pow
    y = zeros((2 * p, n), dtype=complex)
    z = zeros((2 * p, n), dtype=complex)

    j = indices(z.shape)[0]
    y[(p - m) // 2 : (p + m) // 2, :] = x * exp(
        -1.0j * pi * (j[0:m] ** 2) * float(alpha) / m
    )

    z[0:m, :] = exp(1.0j * pi * (j[0:m] ** 2) * float(alpha) / m)
    z[-m:, :] = exp(1.0j * pi * ((j[-m:] - 2 * p) ** 2) * float(alpha) / m)

    d = exp(-1.0j * pi * j ** 2 ** float(alpha) / m) * ifft(
        fft(y, axis=0) * fft(z, axis=0), axis=0
    )

    return d[0:m]


# TODO better docstring 

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 _ncc_c(x, y):
    """
    >>> _ncc_c([1,2,3,4], [1,2,3,4])
    array([ 0.13333333,  0.36666667,  0.66666667,  1.        ,  0.66666667,
            0.36666667,  0.13333333])
    >>> _ncc_c([1,1,1], [1,1,1])
    array([ 0.33333333,  0.66666667,  1.        ,  0.66666667,  0.33333333])
    >>> _ncc_c([1,2,3], [-1,-1,-1])
    array([-0.15430335, -0.46291005, -0.9258201 , -0.77151675, -0.46291005])
    """
    den = np.array(norm(x) * norm(y))
    den[den == 0] = np.Inf

    x_len = len(x)
    fft_size = 1 << (2*x_len-1).bit_length()
    cc = ifft(fft(x, fft_size) * np.conj(fft(y, fft_size)))
    cc = np.concatenate((cc[-(x_len-1):], cc[:x_len]))
    return np.real(cc) / den 

def find_frequencies(data, freq=44100, bits=16):
    """Convert audio data into a frequency-amplitude table using fast fourier
    transformation.

    Return a list of tuples (frequency, amplitude).

    Data should only contain one channel of audio.
    """
    # Fast fourier transform
    n = len(data)
    p = _fft(data)
    uniquePts = math.ceil((n + 1) / 2.0)

    # Scale by the length (n) and square the value to get the amplitude
    p = [(abs(x) / float(n)) ** 2 * 2 for x in p[0:uniquePts]]
    p[0] = p[0] / 2
    if n % 2 == 0:
        p[-1] = p[-1] / 2

    # Generate the frequencies and zip with the amplitudes
    s = freq / float(n)
    freqArray = numpy.arange(0, uniquePts * s, s)
    return list(zip(freqArray, p)) 

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