Cython 的功率谱

Posted 2023-02-25

【中文标题】Cython 的功率谱

【英文标题】：Power spectrum with Cython

【发布时间】：2012-05-27 07:59:13

【问题描述】：

我正在尝试使用 Cython 优化我的代码。它正在做一个功率谱，而不是使用 FFT，因为这是我们在课堂上被告知要做的。我尝试在 Cython 中编写代码，但没有发现任何区别。 这是我的代码

#! /usr/bin/env python
# -*- coding: utf8 -*-

from __future__ import division
cimport numpy as np
import numpy as np
cimport cython

@cython.boundscheck(False)
def power_spectrum(time, data, double f_min, double f_max, double df,w=1 ):

    cdef double com,f
    cdef double s,c,sc,cc,ss
    cdef np.ndarray[double, ndim=1] power
    cdef np.ndarray[double, ndim=1] freq

    alfa, beta = [],[] 
    m = np.mean(data)
    data -= m       

    freq = np.arange( f_min,f_max,df )
    for f in freq:
            sft = np.sin(2*np.pi*f*time)
            cft = np.cos(2*np.pi*f*time)
            s   = np.sum( w*data*sft )
            c   = np.sum( w*data*cft )
            ss  = np.sum( w*sft**2  )
            cc  = np.sum( w*cft**2  )
            sc  = np.sum( w*sft*cft )

            alfa.append( ( s*cc-c*sc )/( ss*cc-sc**2 ))
            beta.append( ( c*ss-s*sc )/( ss*cc-sc**2 ))
            com = -(f-f_min)/(f_min-f_max)*100
            print "%0.3f%% complete" %com

    power = np.array(alfa)**2 + np.array(beta)**2
    return freq,power,alfa,beta

时间和数据通过 numpy.loadtxt 加载并发送到此函数。 当我这样做时

cython -a power_spectrum.pyx

.html 文件很黄，所以效率不高。尤其是整个 for 循环以及计算幂和返回所有内容。

我曾尝试阅读 Cython 的官方指南，但由于我从未用 C 编写过代码，因此有点难以理解。

非常感谢所有帮助:)

【问题讨论】：

为什么不用 C 来编写代码？

我以前从未用 C 编写过代码，所以工作量很大。我喜欢 Python 中的代码，但我希望它更快。

如果没有 print 语句，它的性能会更好吗？您也可能希望将alfa 和beta 实例化为np.zeroes(len(freq)) 而不是列表，将转换保存到最后np.array 并避免附加到列表。

是的，没有 print 语句会好一点。当然，我不应该将它粘贴到我的代码中。制作 alfa 和 beta np.zeros(len(freq)) 的问题是，for 循环中的 f 不是整数。

run a profiler

【参考方案1】：

Cython 可以读取 numpy 数组 according to this，但它不会神奇地编译像 np.sum 这样的东西 - 你仍然只是调用 numpy 方法。

你需要做的是用纯 cython 重写你的内部循环，然后它可以为你编译它。所以你需要重新实现np.sum、np.sin等。预分配aplfa和beta是一个好主意，所以你不要使用append并尝试cdef尽可能多的变量可能。

编辑

这是一个完整的示例，显示了完全 C 编译的内部循环（没有黄色）。我不知道代码是否正确，但它应该是一个很好的起点！特别注意在任何地方使用cdef，打开cdivision 并从标准库中导入sin 和cos。

from __future__ import division
cimport numpy as np
import numpy as np
cimport cython
from math import pi

cdef extern from "math.h":
    double cos(double theta)
    double sin(double theta)

@cython.boundscheck(False)
@cython.cdivision(True)
def power_spectrum(np.ndarray[double, ndim=1] time, np.ndarray[double, ndim=1] data, double f_min, double f_max, double df, double w=1 ):

    cdef double com,f
    cdef double s,c,sc,cc,ss,t,d
    cdef double twopi = 6.283185307179586
    cdef np.ndarray[double, ndim=1] power
    cdef np.ndarray[double, ndim=1] freq = np.arange( f_min,f_max,df )
    cdef int n = len(freq)
    cdef np.ndarray[double, ndim=1] alfa = np.zeros(n)
    cdef np.ndarray[double, ndim=1] beta = np.zeros(n)
    cdef int ndata = len(data)
    cdef int i, j

    m = np.mean(data)
    data -= m       

    for i in range(ndata):
        f = freq[i]

        s = 0.0
        c = 0.0
        ss = 0.0
        cc = 0.0
        sc = 0.0
        for j in range(n):
            t = time[j]
            d = data[j]
            sf = sin(twopi*f*t)
            cf = cos(twopi*f*t)
            s += w*d*sf
            c += w*d*cf
            ss += w*sf**2
            cc += w*cf**2
            sc += w*sf*cf

        alfa[i] = ( s*cc-c*sc )/( ss*cc-sc**2 )
        beta[i] = ( c*ss-s*sc )/( ss*cc-sc**2 )

    power = np.array(alfa)**2 + np.array(beta)**2
    return freq,power,alfa,beta

【讨论】：

感谢您的回复。您能否举例说明如何重新实施例如np.sum 在纯 cython 代码中？当f 不是整数时，我不确定如何不使用append。

alfa 和 beta 与 freq 具有相同的长度，因此您可以以相同的方式定义它们，然后在循环外用零填充它们。重新安排你的循环来循环一个索引i 说它来自0..len(freq)-1 并使用它来读取freq 并存储alfa 和beta。至于np.sum - 它是sft 和data 的循环，逐点相乘然后相加。如果您在 cython 中编写这一切，那么您可以将内部循环的前 8 行组合成一个针对 time 数组索引的 for 循环。

这里是 python 中的一个示例 dft：numericalrecipes.wordpress.com/2009/04/30/… 注意内部循环中嵌套的 fors。事实上，你对 cython 的代码非常想要，给出或接受一些 cdef 语句

虽然很容易让 for 循环在 i 上运行，但我仍然看不出如何以一种智能且快速的方式求和。 alfa、beta 和 power 现在在 .html 文件中的黄色减少了，但仍然与纯 python 一样慢。

给你 :) 祝你有美好的一天。