python实现时间序列信号的频谱、倒频谱以及功率谱

''' 功率谱 ''' from scipy.fftpack import fft, fftshift, ifft from scipy.fftpack import fftfreq import numpy as np import matplotlib.pyplot as plt import xlrd workbook = xlrd.open_workbook(r'E:\speed.xlsx') #读取文件 sheet1 = workbook.sheet_by_index(0) #根据工作表的索引获取sheet工作表页的对象 idx = sheet1.col_values(0) #获取的是列表，需要转换为数组形式 idx = idx[1:] idx = np.array(idx).transpose() fs = 25600 #采样频率 Sampling_points = 32768 #采样点数，fft后的点数就是这个数 df=1/fs #采样间隔时间 # generate original signal t = np.arange(0.0, 1.28, 1/fs) y = idx def get_fft_power_spectrum(y_values, N, f_s, fs_n): f_values = np.linspace(0.0, f_s/fs_n, N//fs_n) fft_values_ = np.abs(fft(y_values)) fft_values = 2.0/N * (fft_values_[0:N//fs_n]) #频率真实幅值分布 # power spectrum 直接法 ps_values = fft_values**2 / N # power spectrum using correlate cor_x = np.correlate(y_values, y_values, 'same') #自相关 cor_X = fft(cor_x, Sampling_points) ps_cor = np.abs(cor_X) ps_cor_values = 10*np.log10(ps_cor[0:N//fs_n] / np.max(ps_cor)) return f_values, fft_values, ps_values, ps_cor_values f_v, fft_v, ps_v, ps_cor_v = get_fft_power_spectrum(y, Sampling_points, fs, 20) fig = plt.figure(figsize=(15, 12)) #原始加噪信号 pltref = plt.subplot(411) plt.plot(t, y) pltref.set_xlabel('Time (seconds)') pltref.set_ylabel('speed') pltref.set_title('Template') #FFT pltref = plt.subplot(412) plt.plot(f_v, fft_v) plt.plot(f_v, fft_v, color='red') pltref.set_xlabel('frequency (Hz)') pltref.set_ylabel('Amplitude') pltref.set_title('FFT') #直接法 (傅立叶变换的平方)/(区间长度) pltref = plt.subplot(413) plt.plot(f_v, ps_v) pltref.set_xlabel('frequency (Hz)') pltref.set_ylabel('power ') pltref.set_title('power spectrum') #自相关法 pltref = plt.subplot(414) plt.plot(f_v, ps_cor_v) pltref.set_xlabel('frequency (Hz)') pltref.set_ylabel('power (dB)') pltref.set_title('power spectrum using correlate') plt.show()