我正在尝试使用 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 编写过代码,所以有点难以理解。
所有帮助都非常准确:)
最佳答案
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
关于python - Cython 的功率谱,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/10772513/