1.功率谱
from scipy.fftpack import fft, fftshift, ifft
from scipy.fftpack import fftfreq
import numpy as np
import matplotlib.pyplot as plt
fs = 1000
num_fft = 1024;
# generate original signal
t = np.arange(0, 1, 1/fs)
f0 = 100
f1 = 200
x = np.cos(2*np.pi*f0*t) + 3*np.cos(2*np.pi*f1*t) + np.random.randn(t.size)
# FFT , 由于采样点数比fs大,多余的用0填充,在fft源码里能看到
Y = fft(x, num_fft)
Y = np.abs(Y)
# 直接法,power spectrum 关键的一步
ps = Y**2 / num_fft
# 间接法 power spectrum using correlate
cor_x = np.correlate(x, x, 'same')
cor_X = fft(cor_x, num_fft)
ps_cor = np.abs(cor_X)
ps_cor = ps_cor / np.max(ps_cor)
plt.figure(figsize=(15, 12),dpi=100)
plt.subplot(411)
plt.plot(x,label='original wave')
plt.legend(loc='best')
# 这里取20log 只是为了转化成db单位
plt.subplot(412)
plt.plot(20*np.log10(Y[:num_fft//2]),label='fft result')
plt.legend(loc='best')
plt.subplot(413)
plt.plot(20*np.log10(ps[:num_fft//2]),label='功率谱')
plt.legend(loc='best')
plt.subplot(414)
plt.plot(20*np.log10(ps_cor[:num_fft//2]),label='间接法')
plt.legend(loc='best')
plt.show()
2.参考链接
添加链接描述