python程序作功率谱、倒频谱

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.参考链接

添加链接描述