DFT в Python занимает значительно больше времени, чем C

Question

В настоящее время я работаю над переводом кода C на Python. Этот код используется, чтобы помочь идентифицировать ошибки, возникающие из алгоритма CLEAN, используемого в радиоастрономии. Чтобы выполнить этот анализ, значения карт преобразования интенсивности Фурье, карты Q Stokes и карты U Stokes должны быть найдены при определенных значениях пикселей (заданных ANT_pix). Эти Карты - только 257*257 массивов.

Приведенный ниже код занимает несколько секунд для запуска с C, но занимает несколько часов для запуска с Python. Я почти уверен, что он ужасно оптимизирован, так как мои знания Python весьма скудны.

Спасибо за любую помощь, которую вы можете оказать.

Обновление Мой вопрос: есть ли лучший способ реализовать циклы в Python, который ускорит процесс. Я прочитал здесь довольно много ответов на другие вопросы по Python, которые рекомендуют избегать вложенных циклов for в Python, если это возможно, и мне просто интересно, знает ли кто-нибудь хороший способ реализации чего-то вроде приведенного ниже кода Python без циклов или с лучшими оптимизированные петли. Я понимаю, что это может быть непростой задачей!

Я использовал FFT до сих пор, но мой супервайзер хочет посмотреть, какую разницу сделает DFT. Это связано с тем, что положение антенны, как правило, не возникает при точных значениях пикселей. Использование FFT требует округления до ближайшего значения пикселя.

Я использую Python в качестве CASA, компьютерная программа, используемая для сокращения наборов данных радиоастрономии, написана на python, и реализация скриптов Python в ней намного проще, чем C.

Оригинальный код

def DFT_Vis(ANT_Pix="",IMap="",QMap="",UMap="", NMap="", Nvis=""):

UV=numpy.zeros([Nvis,6])
Offset=(NMap+1)/2
ANT=ANT_Pix+Offset;

i=0
l=0
k=0
SumI=0
SumRL=0
SumLR=0


z=0

RL=QMap+1j*UMap
LR=QMap-1j*UMap

Factor=[math.e**(-2j*math.pi*z/NMap) for z in range(NMap)]

for i in range(Nvis):
    X=ANT[i,0]
    Y=ANT[i,1]

    for l in range(NMap):

        for k in range(NMap):

            Temp=Factor[int((X*l)%NMap)]*Factor[int((Y*k)%NMap)];

            SumI+=IMap[l,k]*Temp
            SumRL+=RL[l,k]*Temp
            SumLR+=IMap[l,k]*Temp               


        k=1

    UV[i,0]=SumI.real
    UV[i,1]=SumI.imag
    UV[i,2]=SumRL.real
    UV[i,3]=SumRL.imag
    UV[i,4]=SumLR.real
    UV[i,5]=SumLR.imag
    l=1
    k=1
    SumI=0
    SumRL=0
    SumLR=0

return(UV)

Answer 1

Вы, вероятно, должны использовать код преобразования Фурье от numpy, а не писать свой собственный: http://docs.scipy.org/doc/numpy/reference/routines.fft.html

Answer 2

Если вы заинтересованы в повышении производительности вашего скрипта, возможно, вам подойдет Cython.

Answer 3

Я не эксперт по БПФ, но, насколько я понимаю, БПФ - это просто быстрый способ вычисления ДПФ. Поэтому для меня ваш вопрос звучит так, будто вы пытаетесь написать алгоритм сортировки по пузырькам, чтобы увидеть, дает ли он лучший ответ, чем быстрая сортировка. Оба алгоритма сортировки дают одинаковый результат!

Поэтому я ставлю под сомнение вашу основную посылку. Мне интересно, можете ли вы просто изменить округление ваших данных и получить тот же результат из кода SciPy FFT.

Кроме того, согласно моему учебнику по DSP, БПФ может дать более точный ответ, чем вычисление ДПФ на долгом пути, просто потому, что операции с плавающей запятой неточны, и БПФ вызывает меньше операций с плавающей запятой на пути к поиску правильного ответа.

Если у вас есть какой-то рабочий код C, который выполняет необходимые вычисления, вы всегда можете обернуть код C, чтобы он вызывался из Python. Обсуждение здесь: Упаковка библиотеки C в Python: C, Cython или ctypes?

Чтобы ответить на ваш актуальный вопрос: как отметил @ZoZo123, было бы большой победой изменить с range() в xrange(), С range(), Python должен построить список чисел, а затем уничтожить список, когда закончите; с xrange() Python просто делает итератор, который выдает числа по одному. (Но учтите, что в Python 3.x range() делает итератор и нет xrange().)

Кроме того, если этот код не должен интегрироваться с остальной частью вашего кода, вы можете попробовать запустить этот код под PyPy. Это именно тот код, который PyPy может лучше всего оптимизировать. Проблема с PyPy состоит в том, что в настоящее время ваш проект должен быть "чистым" Python, и похоже, что вы используете NumPy. (Существуют проекты для совместной работы NumPy и PyPy, но это еще не сделано.) http://pypy.org/

Если этот код должен интегрироваться с остальной частью вашего кода, то я думаю, что вам нужно взглянуть на Cython (как отметил @Krzysztof Rosiński).