FFT с использованием рекурсивной функции python

Я пытаюсь использовать следующий код для поиска FFT данного списка.

После многих попыток я обнаружил, что этот код работает только для входного списка, содержащего элементы 2^m или 2^m+1.

Не могли бы вы пояснить, почему это так и можно ли изменить его, чтобы использовать входной список, содержащий некоторое другое число элементов. (P.S. У меня есть входной список из 16001 элементов) 

    from cmath import exp, pi

    def fft(x):
        N = len(x)
        if N <= 1: return x
        even = fft(x[0::2])
        odd =  fft(x[1::2])
        T= [exp(-2j*pi*k/N)*odd[k] for k in xrange(N/2)]
        return [even[k] + T[k] for k in xrange(N/2)] + 
        [even[k] - T[k] for k in xrange(N/2)]

    print( ' '.join("%5.3f" % abs(f) 
            for f in fft([1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0])) )

Правка 1 Не могли бы вы, пожалуйста, сказать разницу между предыдущим и следующее определение функции: 

def fft(x):
    N = len(x)
    T = exp(-2*pi*1j/N)
    if N > 1:
        x = fft(x[::2]) + fft(x[1::2])
        for k in xrange(N/2):
            xk = x[k]
            x[k] = xk + T**k*x[k+N/2]
            x[k+N/2] = xk - T**k*x[k+N/2]
    return x

Edit 2: на самом деле этот код(в разделе Edit 1) действительно работает (извините за ошибку в отступе и именовании переменных ранее), поэтому я хочу понять разницу между ними.(это работает и для 16001 элемента!)

3 2
python fft

3 ответа:

Ответ на отредактированную версию

В то время как эта версия: 

from __future__ import print_function

import sys

if sys.version_info.major < 3:
    range = xrange

from cmath import exp, pi

def fft2(x):
    N = len(x)
    T = exp(-2*pi*1j/N)
    if N > 1:
        x = fft2(x[::2]) + fft2(x[1::2])
        for k in range(N//2):
            xk = x[k]
            x[k] = xk + T**k*x[k+N//2]
            x[k+N//2] = xk - T**k*x[k+N//2]
    return x

, по-видимому, работает для числа входов в степени двух: 

import numpy as np
from numpy.fft import fft as np_fft

data = [1, 2, 3, 4]
np.allclose(fft2(data), np_fft(data))

- это True.

Он не дает правильных результатов для различного числа входных данных.

data2 = [1, 2, 3, 4, 5]
np.allclose(fft2(data2), np_fft(data2))

- это False.

Он по-прежнему основан на предположении, что число входов равно степени двух, хотя и не создает исключения.
0

Алгоритм

Этот алгоритм работает только для степени двух чисел входных данных. Взгляните на теорию .

Это улучшенная версия, которая проверяет это предварительное условие: 

from __future__ import print_function

import sys

if sys.version_info.major < 3:
    range = xrange

from cmath import exp, pi

def fft(x):
    N = len(x)
    if N <= 1: 
        return x
    if N % 2 > 0:
        raise ValueError("size of x must be a power of 2")
    even = fft(x[::2])
    odd =  fft(x[1::2])
    r = range(N//2)
    T = [exp(-2j * pi * k / N) * odd[k] for k in r]
    [even[k] for k in r]
    res = ([even[k] + T[k] for k in r] +
           [even[k] - T[k] for k in r])
    return res

input_data = [1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0]
print(' '.join("%5.3f" % abs(f) for f in fft(input_data)))

NumPy поставляется с оптимизированной версией функции FFT. Он работает для размеров, отличных от мощности двух. Но это будет гораздо медленнее, потому что ему нужно применить другой алгоритм.

Например: 

from numpy.fft import fft as np_fft
​
n = 1000
a = np.arange(n ** 2)
b = np.arange(n ** 2 - 10)

Время выполнения для power-of-two дело: 

%timeit np_fft(a)
10 loops, best of 3: 59.3 ms per loop

Гораздо меньше, чем там, где это не так: 

%timeit np_fft(b)
1 loops, best of 3: 511 ms per loop

Предел Рекурсии

Python имеет встроенный предел рекурсии 1000: 

>>> import sys
>>> sys.getrecursionlimit()
1000

Но вы можете увеличить предел рекурсии: 

sys.setrecursion(50000)

The docs расскажут вам, почему: 

getrecursionlimit()

Возвращает текущее значение предела рекурсии, максимальной глубины стека интерпретатора Python. Это ограничение предотвращает бесконечную рекурсию от переполнения стека C и грохот питона. Его можно задать с помощью setrecursionlimit ().

setrecursionlimit()

Установите максимальную глубину стека интерпретатора Python в значение limit. Это ограничение предотвращает бесконечную рекурсию от переполнения стека C и сбоя Python.

Максимально возможный предел зависит от платформы. Пользователю может потребоваться установить более высокий предел, если у него есть программа, требующая глубокой рекурсии, и платформа, поддерживающая более высокий предел. Это следует делать с осторожностью, потому что слишком высокий предел может привести к краху.

0

В приведенном выше коде в этой строке: 

T = [exp(-2j * pi * k / N) * odd[k] for k in r]

Если вы внимательно понаблюдаете за этой частью: exp(-2j * pi * k / N) приведет все ваши выходные векторы по часовой стрелке, для правильного ответа против часовой стрелки замените exp(-2j * pi * k / N) на exp(2j * pi * k / N) как omega = exp(2j * pi * k / N)

0