Python numpy.fft() 使用实例

The following are code examples for showing how to use . They are extracted from open source Python projects. You can vote up the examples you like or vote down the exmaples you don’t like. You can also save this page to your account.

Example 1

def make_node(self, frames, n, axis):
        """
        Compute an n-point fft of frames along given axis.

        """
        _frames = tensor.as_tensor(frames, ndim=2)
        _n = tensor.as_tensor(n, ndim=0)
        _axis = tensor.as_tensor(axis, ndim=0)
        if self.half and _frames.type.dtype.startswith('complex'):
            raise TypeError('Argument to HalfFFT must not be complex', frames)
        spectrogram = tensor.zmatrix()
        buf = generic()
        # The `buf` output is present for future work
        # when we call FFTW directly and re-use the 'plan' that FFTW creates.
        # In that case, buf would store a CObject encapsulating the plan.
        rval = Apply(self, [_frames, _n, _axis], [spectrogram, buf])
        return rval 

Example 2

def perform(self, node, inp, out):
        frames, n, axis = inp
        spectrogram, buf = out
        if self.inverse:
            fft_fn = numpy.fft.ifft
        else:
            fft_fn = numpy.fft.fft

        fft = fft_fn(frames, int(n), int(axis))
        if self.half:
            M, N = fft.shape
            if axis == 0:
                if (M % 2):
                    raise ValueError(
                        'halfFFT on odd-length vectors is undefined')
                spectrogram[0] = fft[0:M / 2, :]
            elif axis == 1:
                if (N % 2):
                    raise ValueError(
                        'halfFFT on odd-length vectors is undefined')
                spectrogram[0] = fft[:, 0:N / 2]
            else:
                raise NotImplementedError()
        else:
            spectrogram[0] = fft 

Example 3

def spike_lfp_filters(spikes,lfp):
    '''
    TODO: documentation
    
    Cross-spectral densities between spikes and LFP
    NTrials,NNeurons,NSamples = np.shape(spikes)
    
    Parameters
    ----------
    
    Returns
    -------
    '''
    NTrials,NNeurons,NSamples = np.shape(spikes)
    # precomute lfp fft
    window = hanning(NSamples)
    lfpmean   = np.mean(lfp)
    spikemean = np.mean(spikes,(0,2))
    fftlfp  = [fft((lfp[t]-np.mean(lfp[t]))*window) for t in np.arange(NTrials)]
    result = np.zeros((NSamples,NNeurons),dtype=np.complex128)
    for i in np.arange(NNeurons):
        cspectra = [np.conj(fft((spikes[t,i,:]-np.mean(spikes[t,i,:]))*window))*fftlfp[t] for t in np.arange(NTrials)]
        result[:,i] = np.mean(cspectra,0)
    return result 

Example 4

def spectreconstruct(k,B,spikes=None,fftspikes=None):
    '''
    TODO: documentation
    
    Reconstructs LFP from spikes using cross-spectral matrix.
    Can optionally pass the fts if they are already available
    NTrials,NNeurons,NSamples = np.shape(spikes)
    
    Parameters
    ----------
    
    Returns
    -------
    '''
    if spikes!=None:
        NTrials,NNeurons,NSamples = np.shape(spikes)
    else:
        NTrials,NNeurons,NSamples = np.shape(fftspikes)
    if ffts==None:
        assert spikes!=None
        fftspikes = np.array([[fft(trial[i]-np.mean(trial[i])) for i in np.arange(NNeurons)] for trial in spikes])
    result = [ifft(sum(fftspikes[t]*B.T,0)) for t in np.arange(NTrials)]
    return result 

Example 5

def diffusion_basis(N=range(1,6),t=np.arange(100)):
    '''
    Note: conceptually similar to other basis functions in this file
    with base=2 and offset=1
    repeatly convolves exponential with itself to generate basis
    
    Parameters
    ----------
    
    Returns
    -------
    '''
    print('THIS IS BAD')
    assert 0
    normalize = lambda x:x/np.sum(x)
    first = np.fft(np.exp(-t))
    kernels = [normalize(np.real(np.ifft(first**(2**(1+(n-1)*0.5))))) for n in N]
    return np.array(kernels) 

Example 6

def delta(N):
    '''
    Parameters
    ----------
    
    Returns
    -------
    '''
    # Get discrete delta but make it even
    delta = np.zeros((N,))
    if (N%2==1):
        delta[N//2]=1
    else:
        delta[N//2+1]=delta[N//2]=0.5
    x = np.fft(delta)
    return x 

Example 7

def fft(a, b=None, axis=0, threads=1, **kw):
    if b is None:
        return numpy.fft.fft(a, axis=axis)
    else:
        b[:] = numpy.fft.fft(a, axis=axis)
        return b 

Example 8

def ifft(a, b=None, axis=0, threads=1, **kw):
    if b is None:
        return numpy.fft.ifft(a, axis=axis)
    else:
        b[:] = numpy.fft.ifft(a, axis=axis)
        return b 

Example 9

def rfft(a, b=None, axis=0, threads=1, **kw):
    if b is None:
        return numpy.fft.rfft(a, axis=axis)
    else:
        b[:] = numpy.fft.rfft(a, axis=axis)
        return b 

Example 10

def irfft(a, b=None, axis=0, threads=1, **kw):
    if b is None:
        return numpy.fft.irfft(a, axis=axis)
    else:
        b[:] = numpy.fft.irfft(a, axis=axis)
        return b 

Example 11

def fft2(a, b=None, axes=(0, 1), threads=1, **kw):
    if b is None:
        return numpy.fft.fft2(a, axes=axes)
    else:
        b[:] = numpy.fft.fft2(a, axes=axes)
        return b 

Example 12

def rfft2(a, b=None, axes=(0, 1), threads=1, **kw):
    if b is None:
        return numpy.fft.rfft2(a, axes=axes)
    else:
        b[:] = numpy.fft.rfft2(a, axes=axes)
        return b 

Example 13

def irfft2(a, b=None, axes=(0, 1), threads=1, **kw):
    if b is None:
        return numpy.fft.irfft2(a, axes=axes)
    else:
        b[:] = numpy.fft.irfft2(a, axes=axes)
        return b 

Example 14

def fftn(a, b=None, axes=(0, 1, 2), threads=1, **kw):
    if b is None:
        return numpy.fft.fftn(a, axes=axes)
    else:
        b[:] = numpy.fft.fftn(a, axes=axes)
        return b 

Example 15

def ifftn(a, b=None, axes=(0, 1, 2), threads=1, **kw):
    if b is None:
        return numpy.fft.ifftn(a, axes=axes)
    else:
        b[:] = numpy.fft.ifftn(a, axes=axes)
        return b 

Example 16

def rfftn(a, b=None, axes=(0, 1, 2), threads=1, **kw):
    if b is None:
        return numpy.fft.rfftn(a, axes=axes)
    else:
        b[:] = numpy.fft.rfftn(a, axes=axes)
        return b 

Example 17

def _xx2k(self ):
        
        """
        Private: oversampled FFT on the heterogeneous device
        
        Firstly, zeroing the self.k_Kd array
        Second, copy self.x_Nd array to self.k_Kd array by cSelect
        Third: inplace FFT
        """
        
        self.cMultiplyScalar(self.zero_scalar, self.k_Kd, local_size=None, global_size=int(self.Kdprod))
        self.cSelect(self.NdGPUorder,      self.KdGPUorder,  self.x_Nd, self.k_Kd, local_size=None, global_size=int(self.Ndprod))
        self.fft( self.k_Kd,self.k_Kd,inverse=False)
        self.thr.synchronize() 

Example 18

def _k2xx(self):
        """
        Private: the inverse FFT and image cropping (which is the reverse of _xx2k() method)
        """
        self.fft( self.k_Kd2, self.k_Kd2,inverse=True)
#         self.x_Nd._zero_fill()
        self.cMultiplyScalar(self.zero_scalar, self.x_Nd,  local_size=None, global_size=int(self.Ndprod ))
#         self.cSelect(self.queue, (self.Ndprod,), None,   self.KdGPUorder.data,  self.NdGPUorder.data,     self.k_Kd2.data, self.x_Nd.data )
        self.cSelect(  self.KdGPUorder,  self.NdGPUorder,     self.k_Kd2, self.x_Nd, local_size=None, global_size=int(self.Ndprod ))
        
        self.thr.synchronize() 

Example 19

def xx2k(self, xx):
        '''
        fft of the image
        input:
            xx:    scaled 2D image
        output:
            k:    k-space grid
        '''
#         dd = numpy.size(self.st['Kd'])      
        output_x = numpy.zeros(self.st['Kd'], dtype=self.dtype,order='C')
#         output_x[crop_slice_ind(xx.shape)] = xx
        output_x.flat[self.KdCPUorder]=xx.flat[self.NdCPUorder]
        k = numpy.fft.fftn(output_x, self.st['Kd'], range(0, self.ndims))

        return k 

Example 20

def k2xx(self, k):
        '''
        Transform regular k-space (Kd array) to scaled image xx (Nd array)
        '''
#         dd = numpy.size(self.st['Kd'])
        
        k = numpy.fft.ifftn(k, self.st['Kd'], range(0, self.ndims))
        xx= numpy.zeros(self.st['Nd'],dtype=dtype, order='C')
        xx.flat[self.NdCPUorder]=k.flat[self.KdCPUorder]
#         xx = xx[crop_slice_ind(self.st['Nd'])]
        return xx 

Example 21

def xx2k(self, xx):
        """
        Private: oversampled FFT on CPU
        
        Firstly, zeroing the self.k_Kd array
        Second, copy self.x_Nd array to self.k_Kd array by cSelect
        Third: inplace FFT
        """
#         dd = numpy.size(self.Kd)      
        output_x = numpy.zeros(self.Kd, dtype=self.dtype,order='C')
#         output_x[crop_slice_ind(xx.shape)] = xx
        output_x.flat[self.KdCPUorder]=xx.flat[self.NdCPUorder]
        k = numpy.fft.fftn(output_x, self.Kd, range(0, self.ndims))

        return k 

Example 22

def _xx2k(self ):
        
        """
        Private: oversampled FFT on the heterogeneous device
        
        Firstly, zeroing the self.k_Kd array
        Second, copy self.x_Nd array to self.k_Kd array by cSelect
        Third: inplace FFT
        """
        
        self.cMultiplyScalar(self.zero_scalar, self.k_Kd, local_size=None, global_size=int(self.Kdprod))
        self.cSelect(self.NdGPUorder,      self.KdGPUorder,  self.x_Nd, self.k_Kd, local_size=None, global_size=int(self.Ndprod))
        self.fft( self.k_Kd,self.k_Kd,inverse=False)
        self.thr.synchronize() 

Example 23

def _k2xx(self):
        """
        Private: the inverse FFT and image cropping (which is the reverse of _xx2k() method)
        """
        self.fft( self.k_Kd2, self.k_Kd2,inverse=True)
#         self.x_Nd._zero_fill()
        self.cMultiplyScalar(self.zero_scalar, self.x_Nd,  local_size=None, global_size=int(self.Ndprod ))
#         self.cSelect(self.queue, (self.Ndprod,), None,   self.KdGPUorder.data,  self.NdGPUorder.data,     self.k_Kd2.data, self.x_Nd.data )
        self.cSelect(  self.KdGPUorder,  self.NdGPUorder,     self.k_Kd2, self.x_Nd, local_size=None, global_size=int(self.Ndprod ))
        
        self.thr.synchronize() 

Example 24

def fftpack_lite_rfftb(buf, s, temp_buf=[()]):
	n = len(buf)
	m = (n - 1) * 2
	temp = temp_buf[0]
	if m >= len(temp): temp_buf[0] = temp = numpy.empty(m * 2, buf.dtype)
	numpy.divide(buf, m, temp[0:n])
	temp[n:m] = 0
	result = numpy.fft.fftpack_lite.rfftb(temp[0:m], s)
	return result 

Example 25

def psf(aperture, wavefront, overfill=1):
  """
  Transform an aperture and the wavefront to a PSF
  
  Args:
      aperture (TYPE): aperture
      wavefront (TYPE): wavefront
      overfill (int, optional): number of extra pixels
  
  Returns:
      real: PSF
  """
  npix = len(wavefront)
  nbig = npix*overfill

  wfbig = numpy.zeros((nbig,nbig),dtype='d')

  half = (nbig - npix)/2
  wfbig[half:half+npix,half:half+npix] = wavefront

  illum = numpy.zeros((nbig,nbig),dtype='d')
  illum[half:half+npix,half:half+npix] = aperture

  phase = numpy.exp(wfbig*(0.+1.j))
  input = illum*phase

  ft = numpy.fft.fft2(input)
  powft = numpy.real(numpy.conj(ft)*ft)

  sorted = numpy.zeros((nbig,nbig),dtype='d')
  sorted[:nbig/2,:nbig/2] = powft[nbig/2:,nbig/2:]
  sorted[:nbig/2,nbig/2:] = powft[nbig/2:,:nbig/2]
  sorted[nbig/2:,:nbig/2] = powft[:nbig/2,nbig/2:]
  sorted[nbig/2:,nbig/2:] = powft[:nbig/2,:nbig/2]

  crop =  sorted[half:half+npix,half:half+npix]

  fluxrat = numpy.sum(crop)/numpy.sum(sorted)
  # print("Cropped PSF has %.2f%% of the flux" % (100*fluxrat))

  return crop 

Example 26

def test_lib(self):
        from pyculib.fft.binding import libcufft
        cufft = libcufft()
        self.assertNotEqual(libcufft().version, 0) 

Example 27

def test_plan1d(self):
        from pyculib.fft.binding import Plan, CUFFT_C2C
        n = 10
        data = np.arange(n, dtype=np.complex64)
        orig = data.copy()
        d_data = cuda.to_device(data)
        fftplan = Plan.one(CUFFT_C2C, n)
        fftplan.forward(d_data, d_data)
        fftplan.inverse(d_data, d_data)
        d_data.copy_to_host(data)
        result = data / n
        self.assertTrue(np.allclose(orig, result.real)) 

Example 28

def test_plan2d(self):
        from pyculib.fft.binding import Plan, CUFFT_C2C
        n = 2**4
        data = np.arange(n, dtype=np.complex64).reshape(2, n//2)
        orig = data.copy()
        d_data = cuda.to_device(data)
        fftplan = Plan.two(CUFFT_C2C, *data.shape)
        fftplan.forward(d_data, d_data)
        fftplan.inverse(d_data, d_data)
        d_data.copy_to_host(data)
        result = data / n
        self.assertTrue(np.allclose(orig, result.real)) 

Example 29

def test_plan3d(self):
        from pyculib.fft.binding import Plan, CUFFT_C2C
        n = 32
        data = np.arange(n, dtype=np.complex64).reshape(2, 2, 8)

        orig = data.copy()
        d_data = cuda.to_device(data)
        fftplan = Plan.three(CUFFT_C2C, *data.shape)
        fftplan.forward(d_data, d_data)
        fftplan.inverse(d_data, d_data)
        d_data.copy_to_host(data)
        result = data / n
        self.assertTrue(np.allclose(orig, result.real)) 

Example 30

def test_against_fft_1d(self):
        from pyculib.fft.binding import Plan, CUFFT_R2C
        N = 128
        x = np.asarray(np.arange(N), dtype=np.float32)
        xf = np.fft.fft(x)
        d_x_gpu = cuda.to_device(x)
        xf_gpu = np.zeros(N//2+1, np.complex64)
        d_xf_gpu = cuda.to_device(xf_gpu)
        plan = Plan.many(x.shape, CUFFT_R2C)
        plan.forward(d_x_gpu, d_xf_gpu)
        d_xf_gpu.copy_to_host(xf_gpu)
        self.assertTrue( np.allclose(xf[0:N//2+1], xf_gpu,
                                      atol=1e-6) ) 

Example 31

def test_against_fft_3d(self):
        from pyculib.fft.binding import Plan, CUFFT_R2C
        depth = 2
        colsize = 2
        rowsize = 64
        N = depth * colsize * rowsize
        x = np.arange(N, dtype=np.float32).reshape(depth, colsize, rowsize)
        xf = np.fft.fftn(x)
        d_x_gpu = cuda.to_device(x)
        xf_gpu = np.zeros(shape=(depth, colsize, rowsize//2 + 1), dtype=np.complex64)
        d_xf_gpu = cuda.to_device(xf_gpu)
        plan = Plan.many(x.shape, CUFFT_R2C)
        plan.forward(d_x_gpu, d_xf_gpu)
        d_xf_gpu.copy_to_host(xf_gpu)
        self.assertTrue(np.allclose(xf[:, :, 0:rowsize//2+1], xf_gpu, atol=1e-6)) 

Example 32

def test_fft_1d_single(self):
        from pyculib.fft import fft
        N = 32
        x = np.asarray(np.arange(N), dtype=np.float32)
        xf = np.fft.fft(x)

        xf_gpu = np.empty(shape=N//2 + 1, dtype=np.complex64)
        fft(x, xf_gpu)

        self.assertTrue( np.allclose(xf[0:N//2+1], xf_gpu, atol=1e-6) ) 

Example 33

def test_fft_1d_double(self):
        from pyculib.fft import fft
        N = 32
        x = np.asarray(np.arange(N), dtype=np.float64)
        xf = np.fft.fft(x)

        xf_gpu = np.zeros(shape=N//2 + 1, dtype=np.complex128)
        fft(x, xf_gpu)

        self.assertTrue( np.allclose(xf[0:N//2+1], xf_gpu, atol=1e-6) ) 

Example 34

def test_fft_2d_single(self):
        from pyculib.fft import fft
        N2 = 2
        N1 = 32
        N = N1 * N2
        x = np.asarray(np.arange(N), dtype=np.float32).reshape(N2, N1)
        xf = np.fft.fft2(x)

        xf_gpu = np.empty(shape=(N2, N1//2 + 1), dtype=np.complex64)
        fft(x, xf_gpu)

        self.assertTrue( np.allclose(xf[:, 0:N1//2+1], xf_gpu, atol=1e-6) ) 

Example 35

def test_fft_2d_single_col_major(self):
        from pyculib.fft import fft
        N2 = 2
        N1 = 8
        N = N1 * N2
        x = np.asarray(np.arange(N), dtype=np.float32).reshape((N2, N1), order='F')
        xf_ref = np.fft.rfft2(x)

        xf = np.empty(shape=(N2, N1//2 + 1), dtype=np.complex64, order='F')
        fft(x, xf)
        self.assertTrue( np.allclose(xf_ref, xf, atol=1e-6) ) 

Example 36

def test_fft_2d_double(self):
        from pyculib.fft import fft
        N2 = 2
        N1 = 32
        N = N1 * N2
        x = np.asarray(np.arange(N), dtype=np.float64).reshape(N2, N1)
        xf = np.fft.fft2(x)

        xf_gpu = np.empty(shape=(N2, N1//2 + 1), dtype=np.complex128)
        fft(x, xf_gpu)

        self.assertTrue( np.allclose(xf[:, 0:N1//2+1], xf_gpu, atol=1e-6) ) 

Example 37

def test_fft_3d_single(self):
        from pyculib.fft import fft
        N3 = 2
        N2 = 2
        N1 = 32
        N = N1 * N2 * N3
        x = np.asarray(np.arange(N), dtype=np.float32).reshape(N3, N2, N1)
        xf = np.fft.fftn(x)

        xf_gpu = np.empty(shape=(N3, N2, N1//2 + 1), dtype=np.complex64)
        fft(x, xf_gpu)

        self.assertTrue( np.allclose(xf[:, :, 0:N1//2+1], xf_gpu, atol=1e-6) ) 

Example 38

def test_fft_3d_double(self):
        from pyculib.fft import fft
        N3 = 2
        N2 = 2
        N1 = 32
        N = N1 * N2 * N3
        x = np.asarray(np.arange(N), dtype=np.float64).reshape(N3, N2, N1)
        xf = np.fft.fftn(x)

        xf_gpu = np.empty(shape=(N3, N2, N1//2 + 1), dtype=np.complex128)
        fft(x, xf_gpu)

        self.assertTrue( np.allclose(xf[:, :, 0:N1//2+1], xf_gpu, atol=1e-6) ) 

Example 39

def test_fft_1d_roundtrip_single(self):
        from pyculib.fft import fft, ifft
        N = 32
        x = np.asarray(np.arange(N), dtype=np.float32)
        x0 = x.copy()
        xf_gpu = np.empty(shape=N//2 + 1, dtype=np.complex64)
        fft(x, xf_gpu)
        ifft(xf_gpu, x)
        self.assertTrue( np.allclose(x / N, x0, atol=1e-6) ) 

Example 40

def test_fft_1d_roundtrip_double(self):
        from pyculib.fft import fft, ifft
        N = 32
        x = np.asarray(np.arange(N), dtype=np.float64)
        x0 = x.copy()
        xf_gpu = np.empty(shape=N//2 + 1, dtype=np.complex128)
        fft(x, xf_gpu)
        ifft(xf_gpu, x)
        self.assertTrue( np.allclose(x / N, x0, atol=1e-6) ) 

Example 41

def test_fft_3d_roundtrip_single(self):
        from pyculib.fft import fft, ifft
        N3 = 2
        N2 = 2
        N1 = 32
        N = N3 * N2 * N1
        x = np.asarray(np.arange(N), dtype=np.float32).reshape(N3, N2, N1)
        x0 = x.copy()
        xf_gpu = np.empty(shape=(N3, N2, N1//2 + 1), dtype=np.complex64)
        fft(x, xf_gpu)
        ifft(xf_gpu, x)
        self.assertTrue( np.allclose(x / N, x0, atol=1e-6) ) 

Example 42

def test_fft_inplace_1d_single(self):
        from pyculib.fft import fft_inplace
        N = 32
        x = np.asarray(np.arange(N), dtype=np.complex64)
        xf = np.fft.fft(x)

        fft_inplace(x)

        self.assertTrue( np.allclose(xf, x, atol=1e-6) ) 

Example 43

def test_fft_inplace_1d_double(self):
        from pyculib.fft import fft_inplace
        N = 32
        x = np.asarray(np.arange(N), dtype=np.complex128)
        xf = np.fft.fft(x)

        fft_inplace(x)

        self.assertTrue( np.allclose(xf, x, atol=1e-6) ) 

Example 44

def test_fft_inplace_2d_single(self):
        from pyculib.fft import fft_inplace
        N1 = 32
        N2 = 2
        N = N1 * N2
        x = np.asarray(np.arange(N), dtype=np.complex64).reshape(N2, N1)
        xf = np.fft.fft2(x)

        fft_inplace(x)

        self.assertTrue( np.allclose(xf, x, atol=1e-6) ) 

Example 45

def test_fft_inplace_2d_double(self):
        from pyculib.fft import fft_inplace
        N1 = 32
        N2 = 2
        N = N1 * N2
        x = np.asarray(np.arange(N), dtype=np.complex128).reshape(N2, N1)
        xf = np.fft.fft2(x)

        fft_inplace(x)

        self.assertTrue( np.allclose(xf, x, atol=1e-6) ) 

Example 46

def test_fft_1d_roundtrip_double_2(self):
        from pyculib.fft import fft_inplace, ifft_inplace
        N = 32
        x = np.asarray(np.arange(N), dtype=np.complex128)
        x0 = x.copy()

        fft_inplace(x)
        ifft_inplace(x)

        self.assertTrue( np.allclose(x / N, x0, atol=1e-6) ) 

Example 47

def test_fft_2d_roundtrip_single_2(self):
        from pyculib.fft import fft_inplace, ifft_inplace
        N2 = 2
        N1 = 32
        N = N1 * N2
        x = np.asarray(np.arange(N), dtype=np.complex64).reshape(N2, N1)
        x0 = x.copy()

        fft_inplace(x)
        ifft_inplace(x)

        self.assertTrue( np.allclose(x / N, x0, atol=1e-6) ) 

Example 48

def test_fft_3d_roundtrip_double(self):
        from pyculib.fft import fft_inplace, ifft_inplace
        N3 = 2
        N2 = 2
        N1 = 8
        N = N3 * N2 * N1
        x = np.asarray(np.arange(N), dtype=np.complex128).reshape(N3, N2, N1)
        x0 = x.copy()

        fft_inplace(x)
        ifft_inplace(x)

        self.assertTrue( np.allclose(x / N, x0, atol=1e-6) ) 

Example 49

def _fft_module(da):
    if da.chunks:
        return dsar.fft
    else:
        return np.fft 

Example 50

def fast_autocorrelate(x):
    """
    Compute the autocorrelation of the signal, based on the properties of the
    power spectral density of the signal.

    Note that the input's length may be reduced before the correlation is
    performed due to a pathalogical case in numpy.fft:
    http://stackoverflow.com/a/23531074/679081

    > The FFT algorithm used in np.fft performs very well (meaning O(n log n))
    > when the input length has many small prime factors, and very bad
    > (meaning a naive DFT requiring O(n^2)) when the input size is a prime
    > number.
    """

    # This is one simple way to ensure that the input array
    # has a length with many small prime factors, although it
    # doesn't guarantee that (also hopefully we don't chop too much)
    optimal_input_length = int(numpy.sqrt(len(x))) ** 2
    x = x[:optimal_input_length]
    xp = x - numpy.mean(x)
    f = numpy.fft.fft(xp)
    p = numpy.absolute(numpy.power(f, 2))
    pi = numpy.fft.ifft(p)
    result = numpy.real(pi)[:x.size / 2] / numpy.sum(numpy.power(xp, 2))
    return result 
点赞

发表评论

电子邮件地址不会被公开。 必填项已用*标注