On 1 September 2015 at 11:38, Joseph Codadeen <jdm...@hotmail.com> wrote: > >> And while you zero-pad, you can zero-pad to a sequence that is a power of >> two, thus preventing awkward factorizations. > > Does numpy have an easy way to do this, i.e. for a given number, find the > next highest number (within a range) that could be factored into small, > prime numbers as Phil explained? It would help if it gave a list, > prioritised by number of factors.
Just use the next power of 2. Pure powers of 2 are the most efficient for FFT algorithms so it potentially works out better than finding a smaller but similarly composite size to pad to. Finding the next power of 2 is easy to code and never a bad choice. To avoid the problems mentioned about zero-padding distorting the FFT you can use the Bluestein transform as below. This pads up to a power of two (greater than 2N-1) but does it in a special way that ensures that it is still calculating the exact (up to floating point error) same DFT: from numpy import array, exp, pi, arange, concatenate from numpy.fft import fft, ifft def ceilpow2(N): ''' >>> ceilpow2(15) 16 >>> ceilpow2(16) 16 ''' p = 1 while p < N: p *= 2 return p def fftbs(x): ''' >>> data = [1, 2, 5, 2, 5, 2, 3] >>> from numpy.fft import fft >>> from numpy import allclose >>> from numpy.random import randn >>> for n in range(1, 1000): ... data = randn(n) ... assert allclose(fft(data), fftbs(data)) ''' N = len(x) x = array(x) n = arange(N) b = exp((1j*pi*n**2)/N) a = x * b.conjugate() M = ceilpow2(N) * 2 A = concatenate((a, [0] * (M - N))) B = concatenate((b, [0] * (M - 2*N + 1), b[:0:-1])) C = ifft(fft(A) * fft(B)) c = C[:N] return b.conjugate() * c if __name__ == "__main__": import doctest doctest.testmod() -- Oscar _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion