2

我目前正在将一些 C 代码翻译成 Python。此代码用于帮助识别射电天文学中使用的 CLEAN 算法引起的错误。为了进行这种分析，强度图、Q Stokes 图和 U Stokes 图的傅里叶变换值必须在特定像素值处找到（由 ANT_pix 给出）。这些地图只是 257*257 的数组。

下面的代码用 C 运行需要几秒钟，但用 Python 运行需要几个小时。我很确定它经过了极大的优化，因为我对 Python 的了解很差。

谢谢你提供的所有帮助。

更新我的问题是是否有更好的方法来实现 Python 中的循环，这将加快速度。我已经在这里阅读了很多关于 Python 的其他问题的答案，建议尽可能避免在 Python 中嵌套 for 循环，我只是想知道是否有人知道一种很好的方法来实现下面的 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)
4

3 回答 3

3

您可能应该使用 numpy 的傅立叶变换代码，而不是自己编写： http: //docs.scipy.org/doc/numpy/reference/routines.fft.html

于 2012-04-08T15:03:55.713 回答
1

如果您有兴趣提高脚本的性能，cython可能是一种选择。

于 2012-04-08T17:50:40.080 回答
1

我不是 FFT 专家，但我的理解是 FFT 只是一种计算 DFT 的快速方法。所以对我来说，你的问题听起来像是你正在尝试编写一个冒泡排序算法，看看它是否能给出比快速排序更好的答案。它们都是会给出相同结果的排序算法！

所以我质疑你的基本前提。我想知道您是否可以仅更改数据的舍入并从 SciPy FFT 代码中获得相同的结果。

此外，根据我的 DSP 教科书，FFT 可以产生比长期计算 DFT 更准确的答案，这仅仅是因为浮点运算不精确，而且 FFT 在找到正确答案的过程中调用的浮点运算更少。

如果你有一些工作的 C 代码来执行你想要的计算，你总是可以包装 C 代码让你从 Python 调用它。此处讨论：在 Python 中包装 C 库：C、Cython 还是 ctypes？

回答您的实际问题：正如@ZoZo123 指出的那样，从 更改为 将是一个巨大的range()胜利xrange()。使用range()，Python 必须构建一个数字列表，然后在完成后销毁该列表；Python 只是创建了xrange()一个迭代器，一次生成一个数字。（但请注意，在 Python 3.x 中，range()创建了一个迭代器并且没有xrange().）

此外，如果此代码不必与您的其余代码集成，您可以尝试在 PyPy 下运行此代码。这正是 PyPy 最能优化的那种代码。PyPy 的问题是目前你的项目必须是“纯”Python，而且看起来你正在使用 NumPy。（有一些项目可以让 NumPy 和 PyPy 一起工作，但这还没有完成。） http://pypy.org/

如果此代码确实需要与您的其余代码集成，那么我认为您需要查看 Cython（如 @Krzysztof Rosiński 所述）。

于 2012-04-08T19:32:47.800 回答