python - 在 python 中使用光谱进行多锥度光谱分析
问题描述
我正在使用 python ( https://pyspectrum.readthedocs.io/en/latest/install.html )上的频谱库进行多锥度分析,但我无法完全理解输出的幅度。
这里有一段代码用于说明:
from spectrum import *
N=500
dt=2*10**-3
# Creating a signal with 2 sinus waves.
x = np.linspace(0.0, N*dt, N)
y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
# classical FFT
yf = fft.fft(y)
xf = np.linspace(0.0, 1.0/(2.0*dt), N//2)
# The multitapered method
NW=2.5
k=4
[tapers, eigen] = dpss(N, NW, k)
Sk_complex, weights, eigenvalues=pmtm(y, e=eigen, v=tapers, NFFT=500, show=False)
Sk = abs(Sk_complex)
Sk = np.mean(Sk * np.transpose(weights), axis=0)
# ploting both the results
plt.plot(xf,abs(yf[0:N//2])*dt*2)
plt.plot(xf,Sk[0:N//2])
两个结果相似,并在 50 和 80 Hz 处找到频率峰值。经典的 FFT 也能找到合适的幅度(1 和 0.5),但多锥度法找不到合适的幅度。在这个例子中,它的重要性大约是 5 倍。有谁知道如何正确显示结果?谢谢
解决方案
据我了解,这里有几个因素在起作用。
首先,要获得功率谱密度的多锥度估计,您应该这样计算:
Sk = abs(Sk_complex)**2
Sk = np.mean(Sk * np.transpose(weights), axis=0) * dt
即,您需要平均功率谱,而不是傅里叶分量。
其次,要获得功率谱,您只需使用 fft 将能量谱除以估计的 N 并乘以 dt (并且您需要 **2 才能从傅立叶分量中获得功率):
plt.plot(xf,abs(yf[0:N//2])**2 / N * dt)
plt.plot(xf,Sk[0:N//2])
最后,应该直接比较的不是功率谱密度中的幅度,而是总功率。你可以看看:
print(np.sum(abs(yf[0:N//2])**2/N * dt), np.sum(Sk[0:N//2]))
非常匹配。
所以你的整个代码变成:
from spectrum import *
N=500
dt=2*10**-3
# Creating a signal with 2 sinus waves.
x = np.linspace(0.0, N*dt, N)
y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
# classical FFT
yf = fft.fft(y)
xf = np.linspace(0.0, 1.0/(2.0*dt), N//2)
# The multitapered method
NW=2.5
k=4
[tapers, eigen] = dpss(N, NW, k)
Sk_complex, weights, eigenvalues=pmtm(y, e=eigen, v=tapers, NFFT=N, show=False)
Sk = abs(Sk_complex)**2
Sk = np.mean(Sk * np.transpose(weights), axis=0) * dt
# ploting both results
plt.figure()
plt.plot(xf,abs(yf[0:N//2])**2 / N * dt)
plt.plot(xf,Sk[0:N//2])
# ploting both results in log scale
plt.semilogy(xf, abs(yf[0:N // 2]) ** 2 / N * dt)
plt.semilogy(xf, Sk[0:N // 2])
# comparing total power
print(np.sum(abs(yf[0:N//2])**2 / N * dt), np.sum(Sk[0:N//2]))
推荐阅读
- jquery - 调用函数完成后执行 AJAX
- c++ - 如何在 C++ 中取消对 grpc::ClientReader 的 Read() 调用?
- javascript - 通过 Android webview 中的 BLOB 链接下载文件
- javascript - 使用 Lodash 获取数组中每个数组的第一个元素
- fuzzing - 达到超时时如何使afl-fuzz不跳过测试用例
- python - 在 python 中调用上下文管理器的多种方法
- excel - 动态条件格式?
- reactjs - 在 React 无状态组件中调用函数
- python - AttributeError:“列表”对象没有属性“解码”
- circleci-2.0 - Circleci 在合并的分支上运行